Physics-informed neural networks in the recreation of hydrodynamic
  simulations from dark matter
Physics-informed neural networks have emerged as a coherent framework for
building predictive models that combine statistical patterns with domain
knowledge. The underlying notion is to enrich the optimization loss function
with known relationships to constrain the space of possible solutions.
Hydrodynamic simulations are a core constituent of modern cosmology, while the
required computations are both expensive and time-consuming. At the same time,
the comparatively fast simulation of dark matter requires fewer resources,
which has led to the emergence of machine learning algorithms for baryon
inpainting as an active area of research; here, recreating the scatter found in
hydrodynamic simulations is an ongoing challenge. This paper presents the first
application of physics-informed neural networks to baryon inpainting by
combining advances in neural network architectures with physical constraints,
injecting theory on baryon conversion efficiency into the model loss function.
We also introduce a punitive prediction comparison based on the
Kullback-Leibler divergence, which enforces scatter reproduction. By
simultaneously extracting the complete set of baryonic properties for the Simba
suite of cosmological simulations, our results demonstrate improved accuracy of
baryonic predictions based on dark matter halo properties, successful recovery
of the fundamental metallicity relation, and retrieve scatter that traces the
target simulation's distribution.
Zhenyu Dai, Ben Moews, Ricardo Vilalta, Romeel Dave
# Physics-informed neural networks in the recreation of hydrodynamic simulations from dark matter

###### Abstract

Physics-informed neural networks have emerged as a coherent framework for building predictive models that combine statistical patterns with domain knowledge. The underlying notion is to enrich the optimization loss function with known relationships to constrain the space of possible solutions. Hydrodynamic simulations are a core constituent of modern cosmology, while the required computations are both expensive and time-consuming. At the same time, the comparatively fast simulation of dark matter requires fewer resources, which has led to the emergence of machine learning algorithms for baryon inpainting as an active area of research; here, recreating the scatter found in hydrodynamic simulations is an ongoing challenge. This paper presents the first application of physics-informed neural networks to baryon inpainting by combining advances in neural network architectures with physical constraints, injecting theory on baryon conversion efficiency into the model loss function. We also introduce a punitive prediction comparison based on the Kullback-Leibler divergence, which enforces scatter reproduction. By simultaneously extracting the complete set of baryonic properties for the Simba suite of cosmological simulations, our results demonstrate improved accuracy of baryonic predictions based on dark matter halo properties, successful recovery of the fundamental metallicity relation, and retrieve scatter that traces the target simulation's distribution.

keywords: galaxies: evolution - galaxies: haloes - methods: analytical - methods: statistical

## 1 Introduction

The \(\Lambda\)CDM model, coined the standard model of cosmology due to its widespread adoption and explanatory power, plays a crucial role in modern cosmology and astrophysics. Galaxy formation and evolution occur within virialized structures resulting from density perturbations in the early Universe, subjected to gravitational collapse (Frenk & White, 2012). This leads to large-scale structure in the form of a cosmic web evolving from a somewhat smooth starting point, with dark matter halos as gravitationally bound overdensities of the postulated main contributor to the matter content of the Universe.

As galaxies form through the condensation of cooling gas within these halos, their resulting baryonic properties share a natural relationship with the dark matter accumulations they live in (see, for example, Rees & Ostriker, 1977; Blumenthal et al., 1984). As dark matter distributions are often sufficient to extract cosmological parameters of interest, \(N\)-body simulations restricted to simulating dark matter particles or grids through gravitational interactions are common in cosmology (Springel et al., 2005; Boylan-Kolchin et al., 2009; Klypin et al., 2011; Riebe et al., 2013; Potter et al., 2017). Options in this regard include direct numerical integration and the inclusion of a scale factor to model the expansion of the Universe via general-relativistic effects. These simulations are computationally cheap and fast to run compared to more complex alternatives (Efstathiou et al., 1985).

In contrast to the success for large-scale structure, modeling galaxies through \(N\)-body simulations has historically been more difficult, as baryonic physics plays a vital role through nonlinear and dissipative processes at this level of granularity (Somerville & Dave, 2015). At the same time, the relevance of the baryonic properties of galaxies is evident by the reason dark matter bears its name; direct observations in the electromagnetic spectrum rely on baryonic processes emitting such radiation. Methods to include the luminous baryonic matter in our simulations are thus required to enable comparisons to these observations in the first place.

Multiple avenues to circumvent these limitations have been developed. These include abundance matching, which connects the halo mass to a range of baryonic properties through the stellar mass. However, this requires the assumption of the preservation of rank
ordering in mass (Reddick et al., 2013). As an alternative to the rank ordering of all galaxies without separating central and satellite galaxies, halo occupancy distribution modeling makes assumptions about the satellite distributions and uses halo mass functions as an input (Berlind & Weinberg, 2002).

The penultimate step in this evolution of complexity are semi-analytic models, which incur a heavier computational burden as a trade-off for using a complete physical framework. The issue these models face is the number of free parameters for which constraints have to be found (see, for example, Somerville & Primack, 1999; Lu et al., 2014). These include, but are not limited to, GALFORM as presented in Cole et al. (2000) and Baugh et al. (2018), GalICS and GalICS 2.0 (see Hatton et al., 2003; Cattaneo et al., 2017), and SAGE (Croton et al., 2016). Simpler analytic formalisms for the evolution of baryonic galaxy properties exist, for example, the bathtub model by Bouche et al. (2010) and the reservoir model by Krumholz & Dekel (2012), as well as the equilibrium model (Dave et al., 2012; Saintonge et al., 2013; Mitra et al., 2015; Mitra et al., 2017).

Hydrodynamic simulations, which require the most computational resources and time, are generally considered the gold standard of an ab initio approach to galaxy formation and evolution (Vogelsberger et al., 2020). Particle-based and grid-based methods exist based on the foundation that baryonic matter can be modeled by treating gas as an ideal fluid. The former rely on computations of discrete masses, or particles, while the latter splits the simulation volume into discrete spaces (Dolag et al., 2008; Somerville & Dave, 2015). Influential recent hydrodynamic simulations include, for example, Illustris and IllustrisTNG (see Genel et al., 2014; Pillepich et al., 2018), EAGLE as introduced in Schaye et al. (2015), and HorizonAGN by Dubois et al. (2016), as well as Mufasa and its successor Simha, the latter of which is used for the lion's share of our experiments and is described in more detail in Section 2.1 (Dave et al., 2016, 2019).

In recent years, the physical sciences have experienced an exponentially rising interest in applying machine learning algorithms (Carleo et al., 2019). These developments include the acceleration of hydrodynamic simulations of galaxy formation and evolution, with Kamdar et al. (2016) providing one of the first examples by predicting baryonic properties in Illustris using extremely randomized trees, commonly abbreviated as 'extra trees', an ensemble model based on decision trees. This heralded the dominance of tree-based ensembles in related research, including Agarwal et al. (2018) for Mufasa and Lovell et al. (2021) for EAGLE, as well as Jo & Kim (2019), McGibbon & Khochfar (2022), and de Santi et al. (2022), for IllustrisTNG.

While these works apply tree-based ensembles, de Santi et al. (2022) also compare other algorithms such as \(k\)-nearest neighbors, light gradient-boosting machines, and feed-forward neural networks, and combine them with a linear regressor for improved predictions. Similarly, de Andres et al. (2022) use random forests, feed-forward neural networks, and natural and extreme gradient boosting on The Three Hundred data described by Cui et al. (2018), and provide an example of favoring boosting over tree-based methods. In this case, extreme gradient boosting is reported to reach the best accuracy, while natural boosting retrieves the most scatter as a probabilistic regressor.

Jepersen et al. (2022) apply graph neural networks to IllustrisTNG. However, they only use the dark matter version and compute baryonic properties with a semi-analytic model, effectively making the predictor emulate such a model instead of a hydrodynamic simulation. Other works exist that focus on alternative algorithms or hybrid approaches. Moster et al. (2021), for example, apply wide and deep neural networks (WDNN) and reinforcement learning. However, the baryonic properties are calculated not through a hydrodynamic simulation but emerge, an empirical model that statistically links properties from surveys instead of simulating baryonic physics (Moster et al., 2018). The prediction problem is flipped by von Marttens et al. (2022), who retrieve dark matter halos from their baryonic properties in IllustrisTNG instead.

Moews et al. (2021), on the other hand, extend the equilibrium model by incorporating largest-progenitor merger trees1, and combine it with extra trees into a hybrid approach. One disadvantage of such analytic models is that they predict a limited set of properties. They can, however, be looped into both the training and prediction steps of machine learning algorithms by first predicting a limited set of properties with the analytic model and then using these predictions together with halo properties to predict the complete set of barbyonics. Related work on hybrid models is done by Hearin et al. (2020), who generate mock catalogs by combining empirical and semi-analytic models. This leads to the weighted Monte Carlo sampling of baseline catalogs to improve statistical realism.

Footnote 1: This might be a constant source of frustration for some readers of papers on these topics, as decision trees and merger trees are, apart from technically both having a tree-like structure, two entirely different concepts.

Stiskalek et al. (2022), similarly to Moster et al. (2021), also use a WDNN approach for IllustrisTNG and HorizonAGN, and include a comparison with extra trees, but with a particular focus on reproducing the intrinsic scatter of the galaxy-halo connection. The reason is a well-known failure of commonly used machine learning algorithms when applied to hydrodynamic simulations to retrieve the expected scatter. The prediction of a Gaussian probability distribution for estimating targets does, of course, impose a constraint on the scatter. Still, probabilistic approaches to scatter reproduction are commonly encountered features in the literature (see, for example, Lehmann et al., 2016; Desmond et al., 2017; Mitra et al., 2017; Cao et al., 2020).

Machine learning methods that rely on standard metrics like the mean squared error (MSE) face certain limitations. The dominant drawback for this application area is that these models aim to recreate the provided data without knowledge of the underlying physical theory. Some prior work tries to solve this issue; for example, the mentioned reports on weighted sampling through baseline libraries by Hearin et al. (2020) and pre-prediction of baryonic subsets to aid the machine learning model (Moews et al., 2021). Another more direct way is injecting theory directly into the learning process of neural network architectures, with such models coined physics-informed neural networks (PINN) (Raissi et al., 2019).

Originally developed for the finding of solutions for partial differential equations, PINNs are rooted in older research on neural networks for ordinary and partial differential equations (see, for example, Dissanayake & Phan-Thien, 1994; Lagaris et al., 1998). In such models, knowledge represented as suitable equations offers domain constraints. Various approaches exist that include physical equations into the loss function, while replacing a lengthy simulation with a machine-learning model has become common practice in many scientific endeavors (Deiana et al., 2022). The line of thinking behind this approach is simple; instead of evaluating the model solely on prediction accuracy through metrics such as the mean squared error, additional loss components can enforce compliance with additional parameter relationships. For a more in-depth overview of this rapidly expanding area of physics-driven deep learning, we refer interested
readers to the reviews by Karniadakis et al. (2021) and Cuomo et al. (2022).

PINNs are not entirely alien to the field of astrophysics in general. Recent work by Mishra and Molinaro (2021) targets the simulation of radiative transfer by minimizing the residual of the underlying transfer equations, while Martin and Schaub (2022) bypasses inefficiencies of gravity models to learn representations of small-body gravity fields directly. Other recent examples include Branca and Pallottini (2022) on solutions for interstellar medium chemistry and the finding of quasinormal modes of nonrotating black holes (Cornell et al., 2022). The development and application of these architectures offer a powerful way to add knowledge to observational evidence to enforce adherence to theoretical models in training the machine learning algorithm.

In this paper, we present the first way to incorporate the concept of PINNs into the active research field of the completion of dark matter-only information in cosmological simulations with baryonic properties. To achieve this goal, we modify the training process of a deep learning model through two different novel extensions of the loss function. The first part is based on the standard approach of injecting physical theory into these models and adds the stellar-to-halo mass relation (SHMR) as described by Moster et al. (2010) into the training process. This double power law is subsequently used by Moster et al. (2018) to parameterize the instantaneous baryon conversion efficiency, meaning the efficiency with which gas is transformed into stars, and provides a theory-oriented constraint.

The second part, which presents another novel addition, forces the model to recreate the scatter of the underlying hydrodynamic simulation by including the Kullback-Leibler divergence (KLD) as an asymmetric distance measure between approximations of the probability distributions for training targets and associated model predictions (Kullback and Leibler, 1951; Ferdosi et al., 2011). In doing so, we show that our extensions of existing machine learning approaches in baryon inpainting from dark matter halo properties are a powerful tool for modern cosmological simulations.

Compared to the baseline model, our results demonstrate improvements in both predictive accuracy and the reproduction of scatter, and show suitable correlations between target variables and model predictions. The contribution pertains to the broader field of machine learning in astrophysics, including adding physical theory into predictive models and using distributional loss components.

The remainder of this paper is structured as follows. In Section 2, we provide an overview of our machine learning approach and data. Section 2.1 describes the Simba suite of cosmological simulations and our dataset. Section 2.2 covers the functionality and justification of our baseline PINN model; Sections 2.3 and 2.4 introduce extensions of the loss function for the SHMR and a distribution comparison between predictions and target values, respectively. Section 3 presents our experiments and their results. Specifically, Section 3.1 explains the fitting of our theory constraint to Simba and the weighting of loss function components. Section 3.2 provides the results for our predictions, and Section 3.3 shows correlations for individual data point accuracy. Section 4 discusses our findings, limitations of our approach, and follow-ups. Lastly, Section 5 provides our conclusions.

## 2 Data and Methodology

### Simulation data from the Simba suite

The Simba simulation models the co-evolution of gas and dark matter within an expanding metric using the Gizmo code (see Hopkins, 2015), which itself is based on the Gadget-2(Springel et al., 2005). It employs the Meshless Finite Mass (MFM) hydrodynamics method, which marries the convenience of a mass-conserving particle-based code with the shock and instability-capturing advantages of a Riemann solver-based scheme.

Many so-called sub-grid processes have been added to Gizmo to model the formation and evolution of galaxies. These include radiative cooling and photoionization heating, chemical enrichment from stellar evolution, the formation of stars and supermassive black holes, the energy release ('feedback') from supernovae and black hole accretion discs, and the growth and destruction of dust. The complete model detailing these sub-grid prescriptions is described in Dave et al. (2016, 2019).

Simba simulations begin in the linear regime at redshift \(z=249\), and are evolved to \(z=0\), meaning today. 151 snapshot outputs are stored at various redshifts along the way. The main Simba run models a random cube of 147 Mpc (comoving) on a side, represented by 10243 gas elements and 10243 dark matter particles. The minimum (adaptive) spatial resolution is 0.7 kpc. The simulation assumes a _Planck_-concordant cosmology (see Planck Collaboration et al., 2016) of \(\Omega_{m}=0.3\), \(\Lambda=0.7\), \(H_{0}=68\) km/s/Mpc, \(\sigma_{8}=0.82\), and \(n_{s}=0.97\). This results in a mass resolution of \(1.8\times 10^{7}\) per gas element and \(9.5\times 10^{7}\) per dark matter particle.

Footnote 3: [https://simba.roe.ac.uk](https://simba.roe.ac.uk)

Each snapshot is analyzed using the Caesar galaxy/halo catalog package. For each halo identified within Simba using Gizmo's native 3D Friends-of-Friends (FoF) finder, Caesar identifies galaxies as collections of stars and dense gas via a 6D FoF algorithm. The most massive galaxy within a halo is defined as the central, and the others are satellites. A large range of physical and photometric properties are computed for each halo and galaxy. For this work, the key galaxy quantities are the stellar mass (\(M_{*}\)), star formation rate (SFR), SFR-weighted gas-phase metallicity (\(Z\)), neutral hydrogen mass (\(M_{\rm HI}\)), molecular hydrogen mass (\(M_{\rm H2}\)), and central supermassive black hole mass (\(M_{\rm BH}\)).

These properties are obtained by summing the relevant particles in each galaxy. For halos, the relevant quantities are the total mass (\(M_{h}\)), the dark matter half-mass radius (\(r_{h}\)), and the dark matter velocity dispersion (\(\sigma_{h}\)). The catalogs are stored as HDF5 files, and Caesar provides a simple yet powerful Python-based access interface. The Simba snapshots and catalogs are all publicly available online for use by the scientific community4.

Footnote 4: [https://simba.roe.ac.uk](https://simba.roe.ac.uk)

Our work uses central galaxies from the n100ni1024_151 version of the Simba main runs. Entries for which \(M_{\rm BH}=0\) are dropped to avoid zero-mass black holes distorting the distributions of predictions in our experiments. Similarly, we restrict the range of included halo masses to \(11\leq\log_{10}(M_{h})\leq 14\), as in related research mentioned in Section 1. Apart from these preprocessing steps, we refrain from any further alterations. While suitable data selections could lead to improved predictions, the goal of this paper is to maintain generalizability. For the model training and prediction processes, we apply min-max normalization to input and target variables to scale values within the same interval \([0,1]\) and then revert predictions to their proper scales.

### Physics-informed neural network model

The concept of PINNs rests on the assumption that data-driven statistical learning can be enhanced with domain knowledge. This inte
gration can be implemented using different frameworks, as covered in the overview of Section 1 and related reviews (Karniadakis et al., 2021; Cuomo et al., 2022). One approach is to generate more data specifically crafted to enforce domain knowledge into the model. While simple, this approach requires large amounts of additional data to cover a broad region of the input-variable space, effectively introducing observational biases into the inductive analysis (Kashefi et al., 2021; Yang and Perdikaris, 2019).

Another approach involves designing specific learning algorithms that embed knowledge into the learning architecture, for example, convolutional neural networks for images and speech (see LeCun and Bengio, 1995), graph neural networks (see Zhou et al., 2020; Wu et al., 2021), and networks for Hamiltonian systems, among others (Jin et al., 2020). Implementing this approach is difficult because embedding physical laws within a neural network architecture is limited to simple processes. Complex processes require architectural designs that cannot be easily realized with current learning frameworks; even relatively simple processes need complex and elaborate designs. Regarding practical applications, a second complication of these purpose-designed models is creating a network specific to a given problem, making the transfer to new domain applications time-consuming.

The third approach we adopt here enriches the loss function to explicitly incorporate constraints, usually as partial differential equations (Raissi et al., 2019; Lagaris et al., 1998). The benefit is a decoupling between the machine learning strategy and the embedded knowledge; the corresponding framework has broader applicability since the learning architecture is built independently of the underlying physical laws. Normally, the loss function is extended to include a measure of the accuracy of each prediction and the degree of alignment with a physical constraint.

The ability of PINNs to generalize well goes beyond approximation theorems. The representational power of neural networks is well known; under limited assumptions, any continuous function can be approximated to an arbitrarily close fit using a neural network with a finite number of hidden nodes and one hidden layer (Cybenko, 1989; Yarotsky, 2017). It should be noted that representational power is not equivalent to generalization power. Adding domain knowledge to the loss function improves the generalizability of PINNs by reducing the bias component of error while keeping the variance component under control (Hastie et al., 2009).

In traditional PINNs, the final loss, \(\mathcal{L}_{f}\), is a weighted combination of two losses; one is the data-driven empirical loss, \(\mathcal{L}_{s}\), and the other is the domain-knowledge constraint, \(\mathcal{L}_{k}\), with respective weights \(w_{s}\) and \(w_{k}\), so that

\[\mathcal{L}_{f}=w_{s}\mathcal{L}_{s}+w_{k}\mathcal{L}_{k}. \tag{1}\]

The first term is the conventional loss obtained in traditional neural networks through methods such as the squared difference between predictions and target values. Training a network NN\({}_{\theta}\) that aims to find a near-optimal parameter vector \(\theta\), for example, through gradient descent, yields a surrogate for a solution to processes such as complex simulations. The corresponding output \(\hat{f}(x)\), where \(x\) is a feature vector, serves to adjust the weights during training by looking at a loss function \(\mathcal{L}_{s}(\hat{f}(x),y)\). Here, \(y\) is the true response variable, or the target value from the underlying hydrodynamic simulation in our case, and \(\hat{f}(x)\) is our model-provided estimate.

The additional term, \(\mathcal{L}_{k}\), adds domain constraints, usually as differential equations. Specifically, the term captures the partial differential residuals as covered in Section 2.3. In our study, the inputs to the neural network are \(x:=[M_{h},r_{h},\sigma_{h}]\) corresponding to dark matter halo mass at present, dark matter half-mass radius, and dark matter halo velocity dispersion, respectively.

The outputs are \(y:=[M_{s},\text{SFR},Z,M_{\text{HI}},M_{\text{H2}},M_{\text{BH}}]\) corresponding to stellar mass, star formation rate, metallicity, neutral and molecular hydrogen masses, and black hole mass, respectively. During training, we aim to minimize the residual sum of squares,

\[\mathcal{L}_{s}=\sum_{i=1}^{N}(y_{i}-\hat{f}(\mathbf{x}_{i}))^{2}. \tag{2}\]

This data loss assumes points \(\{x_{i},y_{i}\}\) sampled at the initial-boundary locations. Domain knowledge minimizes a different loss function, \(\mathcal{L}_{k}\), forcing the final model to obey the constraint from physical knowledge. Here, we assume points \(\{x_{j}\}\) sampled across the entire input space. In the next section, we explain the domain loss in detail. With the defined loss functions, the neural network NN\({}_{\theta}\) is trained to obtain parameters \(\theta\) using efficient optimization methods, such as gradient descent. The weights \(w_{s}\) and \(w_{k}\) enable different contributions to the final loss function and can be tuned automatically as part of the optimization process.

### Inclusion of baryon conversion efficiency

As the name suggests, the stellar-to-halo mass relation links the stellar mass of a given galaxy to the dark matter halo mass. Suitable parameterizations have been shown to reflect the galaxy mass function observed in the third data release of the Sloan Digital Sky Survey (SDSS DR3; see Panter et al., 2007) more closely, as it does not assume a constant SHMR. The instantaneous baryon conversion efficiency \(\epsilon\) is the rate at which gas is transformed into stars. This efficiency can be described by an SHMR parameterized through a double power law model as shown by Moster et al. (2010).

The latter introduces a parameterization that follows observations by avoiding a surplus of galaxies at low and high masses. For two slopes \(\beta\) and \(\gamma\) that are used to determine the decrease in efficiency at lower and higher masses, respectively, the parameterization takes the form

\[\epsilon(M,z=0)=2\epsilon_{N}\left(\left(\frac{M_{h}}{M_{1}}\right)^{-\beta}+ \left(\frac{M_{h}}{M_{1}}\right)^{\gamma}\right)^{-1}. \tag{3}\]

Here, \(\epsilon_{N}\) is the normalization, while \(M_{1}\) denotes the characteristic mass at which the respective efficiency is the same as its normalization. Moster et al. (2018) use this to parameterize the instantaneous baryon conversion efficiency, showing that peak conversion efficiency takes place at halo masses similar to the characteristic mass,

\[M_{\text{max}}=M_{1}\left(\frac{\beta}{\gamma}\right)^{(\beta+\gamma)^{-1}}, \tag{4}\]

with the general assumption of \(\beta,\gamma>0\). The integrated baryon conversion efficiency is dependent on redshift (see Moster et al., 2013), and the mentioned work allows for parameters of the instantaneous efficiency to vary, with

\[\begin{split}\log_{10}M_{1}(z)&=M_{0}+M_{z}(1-(z+ 1)^{-1})\\ &=M_{0}+M_{z}\left(\frac{z}{z+1}\right),\end{split} \tag{5}\]
and with the normalization and slopes given by

\[\begin{split}\epsilon_{N}\left(z\right)&=\epsilon_{0}+ \epsilon_{z}(1-(z+1)^{-1})=\epsilon_{0}+\epsilon_{z}\left(\frac{z}{z+1}\right), \\ \beta(z)&=\beta_{0}+\beta_{z}(1-(z+1)^{-1})=\beta_{0} +\beta_{z}\left(\frac{z}{z+1}\right),\\ \gamma(z)&=\gamma_{0}.\end{split} \tag{6}\]

As we operate at \(z=0\), these considerations are simplified and need only be optimized for single values. In Section 3.1, we will cover this optimization for Simba data using maximum likelihood estimation, with ranges provided by prior research on this parameterization. This view on the SHMR can then be included, for a given dataset size of \(N\), into the loss function of Eq. 1 as

\[\mathcal{L}_{k}=\sum_{i=1}^{N}\left(\frac{\hat{M}_{*}}{M_{h}}-2\epsilon_{N} \left(\left(\frac{M_{h}}{M_{1}}\right)^{-\beta}+\left(\frac{M_{h}}{M_{1}} \right)^{\gamma}\right)^{-1}\right)^{2}. \tag{7}\]

This injection of domain knowledge provides an additional constraint for the model, which subsequent experiments show to benefit the model to recover the mean relation better.

### Constraints from predictive distributions

Reproducing the scatter found in hydrodynamic simulations when predicting baryonic properties based on dark matter halos is a known challenge in the literature (Cui et al., 2018; Stisalek et al., 2022). As described in Section 1, probabilistic approaches relying on parametrized distribution families are common. In this work, we target the reproduction based on the training data distribution directly by introducing a second extension to the standard loss function in Eq. 2 (see Section 2.2).

The Kullback-Leibler divergence (KLD) is a statistical distance measure to assess the difference between two distributions. Introduced by Kullback and Leibler (1951), it has found various applications in astrophysics in recent years (see, for example, Ben-David et al., 2015; Hee et al., 2016; Moews et al., 2019; Nicola et al., 2019). For a given reference distribution \(P\) and proposal distribution \(Q\), the KLD can be written as

\[D_{\mathrm{KL}}(P||Q)=\sum_{x\in\mathcal{X}}P(x)\mathrm{log}\frac{P(x)}{Q(x)}, \tag{8}\]

or, correspondingly, with an integral for absolutely continuous probability distributions. One important point is that the KLD is not a distance metric due to its status as an asymmetric difference measure. This means that

\[D_{\mathrm{KL}}(P||Q)\neq D_{\mathrm{KL}}(Q||P), \tag{9}\]

as it calculates a directional information loss when approximating \(P\) via \(Q\). It also does not satisfy the triangle equality,

\[d(a,c)\leq d(a,b)+d(b,c), \tag{10}\]

with points \(\{a,b,c\}\in M\) for a given metric space \(M\). Given the above, the KLD is applicable only when a 'true' reference distribution is used. Fortunately, this is the case here, as we want to calculate the difference between the distributions of model predictions and their respective targets.

Following Fussell and Moews (2019), who propose the incorporation of the KLD into the loss function in the context of generative modeling - although without implementing the proposal due to conflicting success metrics - we extend Eq. 1 to

\[\mathcal{L}_{f}=w_{s}\mathcal{L}_{s}+w_{k}\mathcal{L}_{k}+w_{\mathrm{KL}} \mathcal{L}_{\mathrm{KL}}. \tag{11}\]

Here, the KLD is part of the overall loss function through a normal assumption placed on the target and prediction distributions, as the loss needs to remain easily differentiable for error backpropagation during training. This means that

\[\mathcal{L}_{\mathrm{KL}}=D_{\mathrm{KL}}\big{(}N(\mu(y),\sigma(y))||\mathcal{ N}(\mu(\hat{y}),\sigma(\hat{y}))\big{)}, \tag{12}\]

for target values \(y\) and corresponding predictions \(\hat{y}\). In doing so, and as described in the overview in Section 1, we impose an additional constraint on the learning process by forcing the model to recreate the approximate distribution of the underlying training data. In later sections, we will see how this avoids an underprediction of tails in the scatter of the baryonic properties of interest.

Figure 1 provides a schematic for the complete model incorporating both the SHMR from Section 2.3 and the KLD, which we dub the 'Hybrid' model going forward. The lower part of the figure shows the neural network model, going from input data containing information from a set of dark matter halo properties to baryonic predictions, \(\hat{y}\), as the model output. The components of the loss function are encased by a dashed line, listing the mean squared error, the KLD introduced in this section, and the SHMR. The respective weights of these components can be selected based on the predictive performance of the neural network, which we cover as part of our experiments in Section 3.1.

The normal assumption made for the respective distributions is, of course, not without fault. However, it provides an approximation that is differentiable for the purpose of backpropagation during training and is computationally cheap enough to result in reasonable training times on small numbers of GPUs without requiring access to stacks in large-scale supercomputing solutions. We reserve part of the discussion in Section 4 to list possible alternatives in related follow-up research.

Figure 1: Illustration of the hybrid model. The loss function used in this framework is a combination of the mean squared error, the Kullback-Leibler divergence, and the stellar-to-halo mass relationship.
## 3 Experiments and Results

### The stellar-to-halo mass relation in Simba

We first train a simple feed-forward neural network, a multilayer perceptron (MLP), using only the MSE in the loss function. We then use this baseline model to predict all six parameters. Our MLP model has five hidden layers with 50 artificial neurons each. As the top panels of Figure 2 show, this basic model can recover the basic pattern of the SHMR represented in Eq. 3. However, the pattern exhibits a diminished variance at lower and higher halo masses, and is expectedly flat around the characteristic \(M_{\text{max}}\) value.

As a result, the predictions show a flatter relation than 1:1 relative to the true distribution from Simba, indicating a bias. We note that Simba's galaxy sample is limited in stellar mass, which results in a diagonal completeness limit in the SHMR. This leads to galaxies with a high \(M_{s}/M_{h}\) ratio for their halo mass being preferentially included in the sample; such asymmetric completeness (known as Malmquist bias) is common in astrophysics, but can be difficult for machine learning algorithms to recover.

Our PINN, on the other hand, is constructed by encoding Eq. 3 into the loss function as shown in Section 2.3. The parameters of the resulting loss component shown in Eq. 7 are adopted from Behroozi et al. (2019). With this extra constraint, the predicted SHMR pattern better recovers the curvature around \(M_{\text{max}}\) in Figure 2. Eq. 3, however, as a double power-law function, cannot introduce the scatter found in the cosmological simulation, and still displays a biased recovery of the true SHMR relation.

To retrieve the scatter, we need a way to provide information as part of the loss function, effectively forcing the model to recreate it. For this, we extend the loss function with an additional component using the KLD, which we introduced in Section 2.4. This distributional divergence measurement has a history of being used in various cosmological applications, as discussed in Section 1, including the study of properties and auxiliary observational data on baryonic physics (Yasin et al., 2022). The resulting predictions from this hybrid scheme track the respective target values more faithfully and with less bias.

We calculate the mean and covariance matrix of the training data and construct a multivariate Gaussian distribution in six dimensions, one for each prediction target property. By including the KLD-based loss in Eq. 8 into the overall loss function, we create a hybrid model combining the PINN with a KLD measurement as written in Eq. 11. Table 1, listing the MSE as well as the coefficient of determination (\(R^{2}\)) and the Pearson correlation coefficient (\(\rho\)), demonstrates the benefit on the predictive power of these models.

We plot the SHMR relative to the halo mass in Figure 2 to provide a visual overview. This confirms that the hybrid approach between extraneous physical knowledge and distributional adherence outperforms both the baseline model and the PINN alone. The latter provides information on the mean trends of the SHMR, while the KLD helps the neural network to mimic the substantial scatter, which an exact equation does not provide. Compared to the baseline MLP, the figure shows that the distributional component helps predict the SHMR more accurately at lower halo masses while tracing the downward scatter at higher halo masses.

To better understand the correlation between targets and predictions, the bottom panels of Figure 2 show the SHMR values from Simba on the horizontal axis, with SHMR values from model predictions on the vertical axis. The result would follow the plot's diagonal line for a perfect model. The visuals, which indicates a better correlation for the PINN compared to the MLP, and for the hybrid model compared to both, are confirmed by the corresponding correlation measurements in Table 1. Apart from the SHMR, the KLD also improves the model's prediction ability on other baryonic properties, which we will cover in Section 3.2.

Lastly, we perform a maximum likelihood estimation (MLE) to optimize the parameters of Eq. 3 and, accordingly, of the PINN part of the loss function in Eq. 7. Here, the likelihood function is the density function regarded as a function of given parameters \(\theta\),

\[\mathcal{L}(\theta)=\sum_{i=1}^{n}f(x_{i}|\theta),\text{ with }\theta\in \Theta, \tag{13}\]

where \(\Theta\) is the corresponding parameter space and \(f(\cdot)\) is the probability density function. We can estimate the optimal combination of parameter values of interest, \(\hat{\theta}\), as

\[\hat{\theta}=\arg\max_{\theta\in\Theta}\mathcal{L}(\theta|x), \tag{14}\]

which results in a choice of parameter values that make the observed data the most probable. Different values in a wide range have been reported on Eq. 3(Behroozi et al., 2019; Kravtsov et al., 2018; Shankar et al., 2017). Instead of directly applying these numbers, they become the initial guess of our estimator, and the range of values is considered a constraint. The role of the MLE can be framed as knowledge learned, molding the known physical equation to best fit the underlying data from Simba within accepted ranges.

In addition, we split off 10% of the dataset as a validation set to perform a grid search on the weights for the loss function components, \(w=[w_{s},w_{k},w_{\text{KL}}]\). While an MLE approach is infeasible due to the model retraining for each run, this allows us to gauge a suitable combination. Each combination is tested for 20 cross-validations using the MSE and the mean absolute percentage error (MAPE), with the results listed in Table 2.

\begin{table}
\begin{tabular}{l c c c} \hline \hline Model & MSE & \(R^{2}\) & \(\rho\) \\ \hline MLP & 0.024 & 0.821 & 0.906 \\ PINN & 0.023 & 0.827 & 0.910 \\ Hybrid (PINN+KLD) & 0.020 & 0.847 & 0.920 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Performance comparison between different models. The table lists the coefficient of determination (\(R^{2}\)) and Pearson’s correlation coefficient (\(\rho\)) for different neural network models corresponding to Figure 3. The hybrid model is a combination of PINN with KLD in the loss.

\begin{table}
\begin{tabular}{l c c} \hline \hline \((w_{l},w_{k},w_{\text{KL}})\) & MSE & MAPE \\ \hline \((1,0.1,0.1)\) & 0.02623 (0.0013) & 1.126 (0.0527) \\ \((1,0.1,1)\) & 0.02583 (0.0023) & 1.123 (0.0476) \\ \((1,0.1,10)\) & 0.03021 (0.0031) & 1.523 (0.0529) \\ \((1,1,0.1)\) & 0.02671 (0.0020) & 1.234 (0.0588) \\ \((1,1,1)\) & 0.02629 (0.0021) & 1.212 (0.0526) \\ \((1,1,10)\) & 0.03365 (0.0028) & 1.387 (0.0483) \\ \((1,1,10,0.1)\) & 0.03711 (0.0011) & 1.527 (0.0339) \\ \((1,10,1)\) & 0.03521 (0.0015) & 1.434 (0.0415) \\ \((1,10,10)\) & 0.03840 (0.0019) & 1.714 (0.0374) \\ \hline \hline \end{tabular}
\end{table}
Table 2: Grid-based optimization of loss function component weights. The table lists the mean squared error and mean absolute percentage error for various loss function weight combinations in the hybrid model. The values in brackets list the associated standard deviations.
### Predicting baryonic properties from dark matter halos

After establishing the effect on SHMR prediction, we expand our analysis to the full set of six baryonic properties, \(\{M_{\ast},\)SFR, \(Z,M_{\rm HI},M_{\rm H2},M_{\rm BH}\}\). To better understand the cause of these results, we calculate \(\rho\) values between the entire set of available parameters in Figure 3.

Figure 4, where these variables are plotted against the halo mass, shows that each parameter is reasonably well-predicted upon visual inspection. In particular, the model excels for \(M_{\ast}\) and \(M_{\rm BH}\), while the performance on SFR and \(M_{\rm H2}\) is subject to a scatter taper at higher halo masses. In the values of Figure 3, we can see that the correlation between the dark matter halo properties used as inputs, meaning \(\{M_{\rm{H}},r_{h},\sigma_{h}\}\), for these variables are considerably lower compared to the rest of the investigated properties. At the same time, good results are equally reflected in strong correlations for stellar and black hole masses. As covered in Section 3.1, this is aided by the additional SHMR constraint in the loss function.

Following this, we explore secondary correlations between different variables, analogous to a similar analysis performed by Agarwal et al. (2018). Figure 3 indicates that galaxies with higher SFR exhibit positive correlations with \(M_{\rm HI}\) and \(M_{\rm H2}\), and a negative correlation with \(Z\). To test whether our model correctly learns the split in the specific star formation rate,

\[{\rm sSFR}=\frac{{\rm SFR}}{M_{\ast}}, \tag{15}\]

we plot these properties against the stellar mass and color data points using the distance in SSFR values from the mean \(M_{\ast}-{\rm sSFR}\) relation, \(\Delta\log_{10}{\rm sSFR}\), in Figure 5. The resulting mean scaling relations are drawn for model predictions and the corresponding data from the underlying Simba simulation, with the former tracing the results of Agarwal et al. (2018) and Dave et al. (2019).

Both mean scaling relations are obtained by fourth-order polynomial fitting, and we provide the 6\({}^{\rm th}\) to 93\({}^{\rm th}\) percentile range for Simba data indicated by grey shading. We use a bin size of 0.2 for the \(\log_{10}M_{\ast}\) values in solar masses along the horizontal axis for the latter. The divergence between confidence intervals and the mean scaling relation for metallicity in Simba at lower stellar masses is an artifact of a small number of data points available at this range, but we include this left-hand interval to exemplify potential peculiarities encountered in such analyses.

As our approach is a prediction problem using approximations, it is not perfectly consistent, just like other works on machine learning for baryonic inpainting into dark matter halos. Here, the underprediction of scatter at higher stellar masses is notable for \(M_{\rm H2}\) compared to Simba. At the same time, the relationships of \(Z\), \(M_{\rm HI}\), and \(M_{\rm H2}\) to the sSFR are preserved around the mean scaling relations, although in a cleaner split than is the case in the target values provided by the Simba cosmological simul

Figure 3: Correlation matrix of input and output parameters. The color bar indicates the Pearson correlation coefficient values, calculated for variables in the Simba dataset used in the presented work.

Figure 2: Top panels: Hexagonal joint histograms of the predicted SHMR from different models and the true SHMR calculated from the Simba test set. From left to right, the panels show the results from the baseline MLP (yellow), the PINN (red), and the hybrid PINN + KLD model (purple). The Simba data is shown in green. Bottom panels: Hexagonal joint histograms of the predicted versus target SHMR for different models and the test set, with color coding as described above. The diagonal line indicates \(x=y\) for the resulting correlation plots.
of secondary correlations between SFR (or gas content), metallicity, and stellar mass, known as the fundamental metallicity relation, is an important success of this machine learning framework.

### Correlations for separate prediction targets

In this section, we further analyze the quality of model outputs. Similar to our visualizations in Section 3.2, we wish to compare the model predictions to the underlying target data from the Simba suite of cosmological simulations, but separately for different baryonic properties. Fig 6 shows kernel density estimates for all six target properties. For the bandwidth optimization, we make use of Scott's rule as the default heuristic in a variety of statistical software packages for a dataset length \(|X|\) with a given dimensionality,

\[\beta=|X|^{-(\mathrm{dim}(X)+4)^{-1}}. \tag{16}\]

In addition, we also plot the same for the baseline multilayer perception used in previous comparisons. These density plots reveal a reasonably close agreement regarding scatter between Simba and our model's predictions for all baryonic properties. While the baseline model is already performing well on the stellar mass and, to a degree, on the black hole mass, it visibly struggles with the scatter of the remainder of the target properties.

Our hybrid model shows major improvements in the distributions of SFR, \(Z\), \(M_{\mathrm{HI}}\), and \(M_{\mathrm{H2}}\). For the multilayer perceptron, the predicted distributions form a sharp peak around the mean, demonstrating especially difficult-to-retrieve scatter for the star formation rate and the molecular hydrogen mass. While the injected physical knowledge on the SHMR does not aid with scatter predictions but instead with the accuracy of associated properties, the KLD loss component contains extra information on distributions that allows for the production of the necessary scatter.

To quantify these improvements, we compare the baseline model and our hybrid approach in Table 3, providing MSE, \(R^{2}\), and \(\rho\) values as for previous comparisons. Here, we can see the expectedly high correlations for \(M_{*}\), \(M_{\mathrm{BH}}\), and, to a lesser degree, \(Z\), comparable to prior research on different hybrid approaches in this area (Moews et al., 2021).

While the latter does not use neural network architectures or direct injection of information into a loss function, they build an additional analytic model into the prediction pipeline and use full largest-progenitor merger trees, which demonstrates the capability of our loss function extensions.

Despite the correct prediction of mean relations, results for SFR, \(M_{\mathrm{HI}}\), and \(M_{\mathrm{H2}}\) are not quite as perfect when viewed next to the remaining baryonic properties. One reason is that we use a

Figure 4: Hexagonal joint histograms of target values and model predictions. The panels show, for the set of six baryonic properties of interest, plots against the dark matter halo mass. Model predictions are in purple, while the Simba data is in green.

Figure 5: Secondary correlations at \(z=0\). The panels show, from top to bottom, neutral hydrogen, molecular hydrogen, and metallicity as a function of the stellar mass. Mean scaling relations are drawn with blue dot-dashed lines for Simba and red dashed lines for our hybrid model. Galaxies are colored by the distance from the mean \(M_{*}-\mathrm{sSFR}\) relation.
single neural network model to predict all variables simultaneously instead of using separate models. The downside is that the model cannot focus solely on a one-dimensional prediction, but this is not desirable, as it also bars the model from learning the connections between different target variables. The result, however, is that properties with a smaller correlation coefficient in terms of input variables, as shown in Figure 3, are more difficult to recover.

For the star formation rate and both hydrogen masses, we can see in Figure 5 that there is more scatter at higher dark matter halo masses for these targets in particular, which usually have higher scatter compared to \(M_{*}\), \(Z\), and \(M_{\rm BH}\)(Shankar et al., 2017). Another reason for prediction errors is likely found in the normal assumption. To calculate the KLD loss and implant it into the backpropagation, we assume a 6D Gaussian distribution, which is not completely true for our data. More complex approximations, however, are beyond the scope of this initial study due to often prohibitive computational costs and the challenge of maintaining differentiability but are discussed in Section 4.

There are potential solutions, which include the injection of more information, for example, on scatter relations, into the loss function. Another pathway is applying more data, such as merger tree information, into the input features. The main purpose of this paper is to show that PINNs, which make use of extraneous physical knowledge, together with a second distributional loss component, can improve model performance in this challenging area and demonstrate the utility that current developments in deep learning approaches can provide for the simulation of baryonic properties.

## 4 Discussion and Limitations

The application of modern machine learning methods to the completion of \(N\)-body information has emerged as a growing area of interest in recent years. Here, the sometimes mentioned unreasonable effectiveness of tree-based ensembles, most commonly random forests and extra trees, is that these models are comparatively simple in their functionality, and yet, as covered in Section 1, are frequently found to outperform neural network architectures in this area. For standard feed-forward frameworks, the universal approximation theorem even guarantees the ability to represent arbitrary functions under very limited assumptions (Cybenko, 1989;

\begin{table}
\begin{tabular}{l c c c c c c} \hline \hline  & \multicolumn{3}{c}{MLP} & \multicolumn{3}{c}{Hybrid} \\ \cline{2-7} Variable & MSE & \(R^{2}\) & \(\rho\) & MSE & \(R^{2}\) & \(\rho\) \\ \hline \(M_{*}\) & 0.023 & 0.827 & 0.901 & 0.020 & 0.847 & 0.920 \\ SFR & 0.638 & 0.050 & 0.267 & 0.423 & 0.080 & 0.284 \\ \(Z\) & 0.052 & 0.277 & 0.527 & 0.044 & 0.352 & 0.593 \\ \(M_{\rm HI}\) & 0.221 & 0.215 & 0.450 & 0.162 & 0.235 & 0.480 \\ \(M_{\rm H2}\) & 0.308 & 0.087 & 0.284 & 0.208 & 0.097 & 0.313 \\ \(M_{\rm BH}\) & 0.276 & 0.466 & 0.682 & 0.214 & 0.502 & 0.709 \\ \hline \hline \end{tabular}
\end{table}
Table 3: Statistical validation for the full experimental run of hybrid model. The table lists MSE, the coefficient of determination (\(R^{2}\)), and Pearson’s correlation coefficient (\(\rho\)) for each output variable.

Figure 6: Density plots for predictions on all six baryonic target properties. The panels show separate kernel density estimates for stellar mass (\(M_{*}\)), star formation rate (SFR), metallicity (\(Z\)), neutral hydrogen (\(M_{\rm HI}\)), molecular hydrogen (\(M_{\rm H2}\)), and black hole mass (\(M_{\rm BH}\)). The corresponding values from the underlying Sigma test dataset are plotted in yellow, baseline multilayer perceptron predictions in green, and results for the hybrid model are in red.
Hornik et al., 1989; Maiorov and Pinkus, 1999). These theoretical capabilities do not, however, touch upon the 'learnability' of these functions, as there are a wide variety of hyperparameters to consider.

While tree-based ensembles are simple, they too are not immune to such limitations, and more recent research has discussed expected error rates as non-monotonous functions concerning the number of ensemble constituents to be set by the user (Probst and Boulesteix, 2018). With the first application of machine learning to this specific area having arguably started with Kamdar et al. (2016), recent works such as Moster et al. (2021), Jespersen et al. (2022), and Stislaklek et al. (2022) have demonstrated the performance potential of suitably chosen neural network architectures.

One issue that remains is the limitations set by the data used for training and prediction in supervised machine learning models. A popular adage from the early days of computer science still widely used is 'garbage in, garbage out', expressing the constraints on arbitrarily powerful models by insufficient datasets. While the results of cosmological simulations are, of course, by no means garbage, they are also limited by the physics they implement and the chosen assumptions and simplifications. PINNs provide a way to go beyond these ingredients, forcing the model to consider further explicitly specified physical relationships.

Our additional implementation of distributional compliance as part of the loss function puts a further constraint on the model, which targets scatter fidelity. The direct incorporation of this additional knowledge into the architecture also goes beyond prior work on the inclusions of analytic models, which relies on physical computations outside of and before the application of machine learning models, as previously implemented by Moews et al. (2021).

At the same time, the predictive power of our architecture relies on the suitability of the included loss function components in terms of both the physical domain knowledge provided and the choice of a density approximation method. The natural extension of our work, although beyond the scope of this paper, is the identification and analysis of additional domain knowledge to be injected into the loss function, thus allowing the model to root its learning in a more complete set of physics. Another pathway is replacing our normal assumption with more complex approximations that allow for non-Gaussian features and multimodal distributions.

In Table 4, we list the dataset sizes and different dark matter properties used as inputs in similar works providing the same performance metrics (Kamdar et al., 2016; Agarwal et al., 2018; Moews et al., 2021). In contrast to this paper, all three comparable studies employ tree-based models, with the largest dataset used by Kamdar et al. (2016). The latter also provide their model, in addition to \(M_{h}\) and \(\sigma_{h}\), with the three spin components, number of dark matter particles bound to the subhalo, and maximum circular velocity in the subhalo, but omit \(r_{h}\) as used in our work. This supplies additional information on the dark matter halo to their model.

Similarly, Agarwal et al. (2018) do not include \(r_{h}\), but add the halo local density, with \(\{M_{h}\}_{i}\) denoting not only the current halo mass at \(z=0\) but also the five preceding snapshots at higher redshifts, which provides information on the recent merger history. Conversely, \(\{\rho\}_{j}\) indicates the set of nearby halo mass densities within radii \(r\in\{200,500,1000\}\) (in kpc) focused on the halo's mass center, thus providing additional information in a similar way. That being said, their dataset is the smallest in this line-up.

The work most closely aligned with ours in terms of inputs is Moews et al. (2021), although largest-progenitor merger trees over the entirety of a halo's evolution are provided instead of only the current halo mass. More importantly, their hybrid approach uses an analytic formalism to pre-compute a subset of baryonic properties with a physical model, feeding the merger trees both into the latter and the subsequent machine learning model.

In Table 5, we list the performance metrics, \(R^{2}\) and \(\rho\), for these related works to enable a discussion of differences in data and model approach. With regard to Kamdar et al. (2016), our model yields a close performance on \(M_{*}\) and \(M_{\text{BH}}\), but with a less robust prediction of SFR and \(Z\). This could be due to the size of the Illustris dataset used by the authors, which is about eight times larger than our available Simba data, as more suitable data usually leads to a better performance in machine learning models (Zhou et al., 2014).

Agarwal et al. (2018) use Mufasa, the predecessor of Simba, and benefits from their inclusion of prior halo masses and halo mass densities at different distance radii. Our model's results are comparable for \(M_{*}\) and \(M_{\text{HI}}\), with a slight underperformance for \(Z\), but there is a notable degradation in \(MH2\) and SFR predictions. Based on the correlation matrix in Figure 3, these two properties are highly correlated, meaning that an improvement in one in future research should impact the other. One potential reason for this degradation is the larger scatter in Simba higher halo masses, while Agarwal et al. (2018) pre-select their galaxies to be star-forming which strongly reduces the scatter at high masses.

Lastly, Moews et al. (2021) develop a hybrid approach, combining an extra trees ensemble with the equilibrium model, and incorporating merger trees into the latter. The same data and inputs are used in both work, save for the largest-progenitor merger trees that are fed into both the physical and machine learning models. While our model achieves close results on \(M_{*}\), \(Z\), \(M_{\text{HI}}\), and, to a lesser degree, \(M_{\text{BH}}\), metrics of SFR and \(M_{H2}\) are lower. Here, the difficulty to predict the star formation rate due to the large scatter could be data-driven, while Moews et al. (2021) utilize a physical model for this property. At the same time, our model outperforms

\begin{table}
\begin{tabular}{c c c} \hline \hline  & size & input \\ \hline Kamdar et al. (2016) & 249370 & \(M_{h}\), \(S_{X}\), \(S_{y}\), \(S_{z}\), \(\sigma_{h}\), \(N_{h}\), \(V_{\text{c}}\) \\ Agarwal et al. (2018) & 3400 & \(\{M_{h}\}_{i}\), \(\{\rho_{h}\}_{j}\), \(\lambda_{h}\), \(\sigma_{h}\) \\ Moews et al. (2021) & 13132 & \(\{M_{h}\}_{i}\), \(r_{h}\), \(\sigma_{h}\) \\ Dai et al. (this work) & 14247 & \(M_{h}\), \(r_{h}\), \(\sigma_{h}\) \\ \hline \hline \end{tabular}
\end{table}
Table 4: Dataset size and input variables for related works and this paper. Here, \(S_{X,y,z}\) denotes the different components of the spin, \(V_{\text{c}}\) the maximum circular velocity in the subhalo, \(N_{h}\) the number of dark matter particles bound to the subhalo, \(\lambda_{h}\) the halo spin, and \(\rho_{h}\) the halo local density. Sets for additional information at higher redshifts are indicated using braces.

\begin{table}
\begin{tabular}{l c c c c c} \hline \hline  & Kamdar et al. & & Agarwal et al. & Moews et al. \\ \cline{2-6} Variable & \(R^{2}\) & \(\rho\) & \(R^{2}\) & \(\rho\) & \(R^{2}\) & \(\rho\) \\ \hline \(M_{*}\) & 0.91 & 0.95 & 0.90 & 0.95 & 0.82 & 0.94 \\ SFR & 0.63 & 0.79 & 0.55 & 0.74 & 0.73 & 0.87 \\ \(Z\) & 0.93 & 0.96 & 0.73 & 0.86 & 0.21 & 0.66 \\ \(M_{\text{gas}}\) & 0.67 & 0.85 & - & - & - & - \\ \(M_{\text{HI}}\) & - & - & 0.35 & 0.59 & 0.36 & 0.65 \\ \(M_{\text{H2}}\) & - & - & 0.51 & 0.71 & 0.54 & 0.75 \\ \(M_{\text{BH}}\) & 0.72 & 0.85 & - & - & 0.71 & 0.88 \\ \hline \hline \end{tabular}
\end{table}
Table 5: Performance metrics for related works. The table lists the coefficient of determination (\(R^{2}\)) and Pearson’s correlation coefficient (\(\rho\)) for each output variable. Here, in addition to the variables described throughout this paper, \(M_{\text{gas}}\) denotes the total gas mass of a subhalo in the indicated paper.
these results on \(Z\), which reflects the equilibrium model's difficulty with this quantity.

Overall, unlike our method focusing on the retrieval of accurate scatter, the above models are optimized for accuracy and make use of more information on the dark matter halo. The results of both Agarwal et al. (2018) and Moews et al. (2021) show a considerably less accurate scatter retrieval than can be seen in our experiments, with noticeable underprediction of the tails. The goal of this work is not to replace existing research but to demonstrate modifications that can increase the predictive ability of deep learning models in this area. Since the dataset, model, and purpose are different, this comparison provides a direction for future research aimed at combining these strengths.

For follow-up research not using PINNs, or machine learning models that do not require the loss function to be differentiable for backpropagation as described in Section 2.4, we recommend extensions to the comparison of target and prediction distributions. This work uses the KLD to calculate the difference between those distributions under the normal assumption. While a reasonably close proposition for the data used in our experiments, this limits our model's effectiveness when directly transferred to markedly non-Gaussian datasets. In such cases, the spread of predictions through the variance would still be enforced, but in the case of, for example, a starkly multimodal distribution, the recreation of these distributional features would not be a major component of the optimization.

As the comparison needs to be reasonably fast, approaches such as Bayesian mixture models as well as associated methods that are more complex (see, for example, Moews & Zuntz 2020) are likely to slow the training process down too much. As a compromise, a kernel density estimate (KDE), also known as the Parzen-Rosenblatt window after Rosenblatt (1956) and Parzen (1962), can be used, although this is limited to lower dimensionalities. We propose two ways to circumvent the latter limitation. The first is applying one-dimensional KDEs on a variable-by-variable basis and then averaging the Kullback-Leibler divergences between these estimates. While the advantage is the good fit in \(\mathbb{R}^{1}\), this can lead to a subset of variables not being forced to follow the target distribution as long as the KLD average remains small.

The alternative is to make use of dimensionality reduction methods such as principal component analysis, which collapses the \(n\)-dimensional coordinate space, for \(n\) variables, into orthogonal vectors ranked by their ability to explain the variance (for a recent overview, see Jolliffe & Cadima 2016). Reducing the coordinate space to a subspace in \(\mathbb{R}^{2}\) would, for example, still allow for reasonably good KDE approximations while retaining each variable's contribution to a combined KLD. The same line of thought does, of course, also apply to other density approximation and dimensionality reduction methods.

## 5 Conclusion

In this paper, we transfer the paradigm of physics-informed neural networks to predicting baryonic properties for associated dark matter halo variables. We adapt this approach in two different ways. The first includes the stellar-to-halo mass relation as a double power law previously used to parameterize the instantaneous baryon conversion efficiency. While this is a more established way to include physical theory into PINNs, our second extension is the enforcement of baryonic scatter in simulations under a normal assumption using the Kullback-Leibler divergence between the underlying cosmological simulation and predictions. In doing so, we solve the common problem of scatter reproduction in this area, which is merged directly into the machine learning model.

We first test the improvement for the more traditional approach to PINNs, meaning the injection of physical domain knowledge into the loss function, and demonstrate a positive effect on the model's performance. The hybrid approach combines these strengths by including the measurement of distributional differences and outperforms the standard PINN model. Subsequent tests of the scatter retrieval show more faithful reproductions for baryonic properties. These improvements in scatter are especially notable for molecular hydrogen masses and star formation rates, but can also be seen for neutral hydrogen masses and metallicities. In particular, our model successfully recovers the fundamental metallicity relation.

Our experiments demonstrate that PINNs, a rapidly expanding area of research across various subfields of physics, offer a way to directly bake theoretical constraints and distributional adherence into neural network architectures when painting baryonic properties into galactic dark matter halos. As such, they can be used to complete cosmological \(N\)-body simulations based on full hydrodynamic simulation suites, although this comes with the same caveats as other research in this area.

The inference of physics from simulations operates under the assumption that such simulations are a sufficiently close approximation of the real world. Any machine learning models learning from those simulations are subject to the same assumption. That said, including physical models in the learning process enables these algorithms to include domain information beyond the underlying cosmological simulations.

Potential follow-ups include additional physical models specific to galaxy formation and evolution into the loss function, as well as further constraints based on observational data or other simulations to diversify the data sources. Our presented framework is widely applicable to large-scale cosmological simulations and the study of the utility and effect of physical domain knowledge on galaxy evolution emulators. It provides a further piece in the puzzle of fully using modern machine learning in astrophysics.
