Learning the galaxy-environment connection with graph neural networks
Galaxies co-evolve with their host dark matter halos. Models of the
galaxy-halo connection, calibrated using cosmological hydrodynamic simulations,
can be used to populate dark matter halo catalogs with galaxies. We present a
new method for inferring baryonic properties from dark matter subhalo
properties using message-passing graph neural networks (GNNs). After training
on subhalo catalog data from the Illustris TNG300-1 hydrodynamic simulation,
our GNN can infer stellar mass from the host and neighboring subhalo positions,
kinematics, masses, and maximum circular velocities. We find that GNNs can also
robustly estimate stellar mass from subhalo properties in 2d projection. While
other methods typically model the galaxy-halo connection in isolation, our GNN
incorporates information from galaxy environments, leading to more accurate
stellar mass inference.
John F. Wu, Christian Kragh Jespersen
# Learning the galaxy-environment connection with graph neural networks

###### Abstract

Galaxies co-evolve with their host dark matter halos. Models of the galaxy-halo connection, calibrated using cosmological hydrodynamic simulations, can be used to populate dark matter halo catalogs with galaxies. We present a new method for inferring baryonic properties from dark matter subhalo properties using message-passing graph neural networks (GNNs). After training on subhalo catalog data from the Illustris TNG300-1 hydrodynamic simulation, our GNN can infer stellar mass from the host and neighboring subhalo positions, kinematics, masses, and maximum circular velocities. We find that GNNs can also robustly estimate stellar mass from subhalo properties in \(2d\) projection. While other methods typically model the galaxy-halo connection in isolation, our GNN incorporates information from galaxy environments, leading to more accurate stellar mass inference.

Machine Learning, Galaxy-halo connection, Galaxy-halo connection, Galaxy-halo connection, Galaxy-halo connection, Galaxy-halo connection, Galaxy-halo

## 1 Introduction

In the current \(\Lambda\)CDM paradigm of hierarchical galaxy formation, the galaxy-halo connection is crucial for understanding how galaxies form and evolve, and for constraining the small-scale clustering of matter (Somerville and Dave, 2015; Wechsler and Tinker, 2018; Vogelsberger et al., 2020). Techniques for modeling the co-evolution of galaxies and dark matter range from simple, non-parametric approaches to full-physics magnetohydrodynamic simulations which require \(>10^{8}\) CPU hours of computation (e.g., Vale and Ostriker, 2004; Pillepich et al., 2018). Detailed simulations contribute important insights into galaxy formation, but due to their complexity and heavy computational costs, they are hard to analyze and cannot be performed for cosmologically significant volumes. Machine learning (ML) is a natural option for making progress on both of these problems.

We present an equivariant Graph Neural Network (GNN), which takes as its input a graph composed of halos linked on a linking scale of 5 Mpc, and predicts baryonic properties. The GNN incorporates the effects of a galaxy's environment, thereby improving the prediction of its baryonic properties compared to traditional methods. We are also able to train a network on the Illustris TNG300-1 box in 10 minutes on a single NVIDIA A10G GPU; inference takes one second. In this work, we focus on estimating stellar mass from a catalog of subhalo positions, velocities, \(M_{\rm halo}\), and \(V_{\rm max}\).

## 2 Related work

The connection between galaxies and their dark matter halos has been characterized via abundance matching or halo occupation distribution (HOD) models of central halos (Berlind and Weinberg, 2002; Wechsler et al., 2002), conditional luminosity or mass functions (Yang et al., 2003; Moster et al., 2010), subhalo abundance matching (Kravtsov et al., 2004; Vale and Ostriker, 2004; Conroy et al., 2006), and empirical models of the galaxy-halo connection (e.g., Reddick et al., 2013; Behroozi et al., 2019). Several works have also attempted to perform abundance matching or paint baryons (i.e., stars) onto dark matter maps by using classical machine learning algorithms (e.g., Kamdar et al., 2016; Agarwal et al., 2018; Calderon and Berlind, 2019) and/or neural networks (e.g., Zhang et al., 2019; Moster et al., 2021; Mohammad et al., 2022).

In general, these previous methods treat halo/galaxy systems as unrelated entities with no formation history. To rectify this, Villanueva-Domingo et al. (2022) construct mathematical graphs to represent group halos, and train a GNN to learn the central halo mass, which was later applied to estimate the halo masses of local Group galaxies (Villanueva-Domingo et al., 2021). GNNs have also been successfully used to model the dependence of galaxy properties on merger history (e.g., Jespersen et al., 2022; Tang and Ting, 2022), and generate synthetic galaxy catalogs (Jagvaral et al., 2022).

In cosmology, several works have already demonstrated the representational power of GNNs, and have used it for simulation-based inference (likelihood-free inference).
Villanueva-Domingo & Villaescusa-Navarro (2022) employ GNNs to infer the cosmological parameters \(\Omega_{m}\) and \(\sigma_{8}\), using \(3d\) galaxy positions and stellar properties from the CAMELS simulation suite (Villaescusa-Navarro et al., 2021). Makinen et al. (2022) show that GNNs can optimally extract and compress catalog data for cosmological parameter inference. Shao et al. (2023) and de Santi et al. (2023) train GNNs to infer cosmological parameters from dark matter-only simulations, and then validate their robustness on other \(N\)-body and hydrodynamic simulations.

## 3 Cosmic graphs

### Simulation data

We use SUBFIND\(z=0\) subhalo catalogs (Springel et al., 2001) derived from the Illustris TNG300-1 hydrodynamic simulation (Nelson et al., 2019; Pillepich et al., 2019). We split the full cosmological box into \(6^{3}=216\) subvolumes in order to fit into 16 GB of memory, such that each subvolume is about \((50~{}\mathrm{Mpc})^{3}\). For consistency with the TNG simulations, we adopt the Planck Collaboration et al. (2016) cosmology and set \(H_{0}=67.74~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}\).

We select unflagged subhalos that have more than 50 star particles, \(\log(M_{\star}/M_{\odot})>9\), and \(\log(M_{\mathrm{halo}}/M_{\odot})>10\). Due to cosmic variance, some subvolumes only have a few hundred subhalos, while others have thousands. In Figure 1, we show an example of a typical subvolume.

### Equivariant graph neural networks

We construct a mathematical graph for each TNG300 subvolume, such as the one depicted in Figure 1. We designate \(\mathcal{V}_{i}=(\mathbf{x}_{i},\mathbf{v}_{i},M_{\mathrm{halo},i},V_{\mathrm{max},i})\) as the eight node features. Subhalos within a linking length of \(L=5\) Mpc are connected with edges. Subvolumes are padded by \(2.5\) Mpc on each side, such that subvolumes do not share connections that would be relevant for the linking length. We allow nodes to be connected to themselves (i.e., self-loops). On each edge \(\mathcal{E}_{ij}\), we compute three features: the squared Euclidean distance \(d_{ij}\equiv||\mathbf{x}_{i}-\mathbf{x}_{j}||\), the inner product between unit vectors \(\mathbf{e}_{i}\cdot\mathbf{e}_{j}\), and the inner product between unit vectors \(\mathbf{e}_{i}\cdot\mathbf{e}_{i-j}\), where unit vectors \(\mathbf{e}_{i}\equiv(\mathbf{x}_{i}-\bar{\mathbf{x}})/||\mathbf{x}_{i}-\bar{\mathbf{x}}||\) are defined using positions \(\mathbf{x}_{i}\) relative to the centroid of the point cloud distribution \(\bar{\mathbf{x}}\), and \(\mathbf{e}_{i-j}\) is the unit vector in the direction of \(\mathbf{x}_{i}-\mathbf{x}_{j}\).

We use a message-passing GNN based on interaction networks (Battaglia et al., 2016, 2018), similar to the model used by Villanueva-Domingo et al. (2022). By design, the GNN is equivariant to permutations and invariant under the \(E(3)\) group action, i.e., invariant to rotations, reflections, and translations. For more details about equivariant GNNs, see the appendices of Garcia Satorras et al. (2021) and Sections 3.1 and 3.2 of Villanueva-Domingo & Villaescusa-Navarro (2022). We aggregate layer inputs at each node by max pooling over information from neighboring nodes.1 Our GNN has one set of fully connected layers with 256 latent channels and 128 hidden channels. We predict two quantities for each node, which correspond to the logarithmic stellar mass \(y_{i}\equiv\log(M_{\star,i}/M_{\odot})\) and the logarithmic variance, \(\log\Sigma_{i}\) (i.e., the logarithm of the squared uncertainty on stellar mass).

Footnote 1: We do not find significant improvements by using a concatenation of sum, max, mean, and variance aggregations, or by using learnable aggregation functions.

### Optimization

Our loss function is composed of two terms: the mean squared error on the logarithmic stellar mass \(||\hat{\mathbf{y}}-\mathbf{y}||^{2}\), and the squared difference between the predicted and measured variance \(||\hat{\mathbf{\Sigma}}-(\hat{\mathbf{y}}-\mathbf{y})^{2}||^{2}\). The latter term ensures that the variance is appropriately estimated (see Moment Networks, described in Section 2 of Jeffrey & Wandelt, 2020). We stabilize training by taking the logarithm of each loss term before summing them. We monitor the loss as well as the root mean squared error (RMSE) on \(\log(M_{\star}/M_{\odot})\).

We perform \(k=6\)-fold cross-validation. For each fold, we train on 180 subvolumes and validate on 36 subvolumes, such that the validation set forms a \(\sim 50\times 300\times 300\) Mpc\({}^{3}\) subbox. We augment the training data set by adding random noise, sampled from a normal distribution with \(10^{-5}\) times

Figure 1: Cosmic graph of a TNG300 subvolume spanning approximately \(45\) Mpc. Subhalos live on nodes and are colored by their logarithmic subhalo mass. Edges are formed between pairs of subhalos separated by less than the linking length of 5 Mpc.
the standard deviation, for each node variable. Based on a preliminary hyperparameter search, we implement a simple optimization schedule over a total of 1000 epochs using the AdamW optimizer (Kingma & Ba, 2014; Loshchilov & Hutter, 2017) and a batch size of 36. We begin with a learning rate of \(10^{-2}\) and weight decay of \(10^{-4}\), and then decrease both by a factor of \(5\) at 500 epochs, and again decrease both by a factor of \(5\) at 750 epochs. We inspect the training and validation losses to ensure that the optimization is converged and does not overfit the training data.

## 4 Results

Overall, we find that the GNN can infer the stellar mass from subhalo properties with remarkable accuracy. We recover the galaxy stellar mass to within RMSE \(=0.129\) dex of its simulated value by using a GNN. The predictions are largely unbiased as a function of mass.

### Comparisons against baseline models

In Figure 2, we compare the performance of different models trained and cross-validated on the same TNG300 data set. The panels show, from left to right: (a) a subhalo abundance matching (SHAM) model, (b) a random forest (RF) trained using \(M_{\rm halo}\) as input, (c) a RF trained using \(V_{\rm max}\), (d) a RF trained using both \(M_{\rm halo}\) and \(V_{\rm max}\), and (e) a GNN trained using \(3d\) positions, \(3d\) velocities, \(M_{\rm halo}\), and \(V_{\rm max}\). In Table 1, we list performance metrics for various RF and GNN models, including the RMSE, mean average error (MAE), normalized median absolute deviation (NMAD),2 Pearson correlation coefficient (\(\rho\)), correlation of determination (\(R^{2}\)), bias, and outlier fraction (\(>3\times\) NMAD).

Footnote 2: We define NMAD(\(\mathbf{x}\)) \(\equiv\)\(k\cdot{\rm median}(|x-{\rm median}(\mathbf{x})|)\), where \(k\)\(\equiv\)\(1.4826\) ensures that the NMAD and standard deviation are equal for a normally distributed \(\mathbf{x}\).

The SHAM model constructs separate monotonic relationships between \(M_{\rm halo}\) or \(V_{\rm max}\) and \(M_{\star}\) for centrals and satellites. Another difference between the SHAM model and other approaches considered here is the former's explicit treatment of subhalo centrality. In order to facilitate an apples-to-apples comparison, we also train an abundance matching (AM) model that does not distinguish between satellites and centrals; however the AM model performs considerably worse than the SHAM counterpart. We note that the AM and SHAM models are trained and evaluated on the same data set, so their performance metrics may be overinflated.

We also train several RF models, which serve as reasonable proxies for AM or conditional luminosity function models (Calderon & Berlind, 2019). By comparing panels (b) and (c), we observe that \(V_{\rm max}\) is more physically connected to \(M_{\star}\) than \(M_{\rm halo}\), in agreement with previous findings (i.e., Conroy et al. 2006; Reddick et al. 2013; we find this to be true for the RF, AM, and SHAM models). A RF trained on both \(M_{\rm halo}\) and \(V_{\rm max}\) provides an even better reconstruction (RMSE = 0.148 dex).

Ultimately, we find that the GNN strongly outperforms all baseline models. While the GNN does not distinguish between centrals and satellites, it may be able to learn whether a given subhalo is a central based on surrounding subhalo properties (see Section 5.2).

### Centrals versus satellites

Satellite dark matter halos are preferentially stripped relative to stars in a host halo's tidal field (Smith et al., 2016). In Appendix A, we show the stellar mass-halo mass relation for satellite and central galaxies in TNG300 (Figure 3). Indeed, we observe that satellite galaxies exhibit significantly more dispersion than centrals \(M_{\star}\)-\(M_{\rm halo}\) relation. Our \(3d\) GNN is also worse at predicting \(\log(M_{\star}/M_{\odot})\) for satellites than for centrals (see bottom two rows of Table 1), but this is due to the inherently larger scatter in the satellite-halo relation. We find that there is an overall negative bias for satellites and and positive bias for centrals, because the GNN must learn separate offset relations for both centrals and satellites.

Figure 2: Predicted stellar mass versus true stellar mass for the TNG300 data. From left to right, we show results for a subhalo abundance matching (SHAM) model, three random forest (RF) models, and our \(3d\) GNN trained using \(\mathbf{x}\), \(\mathbf{v}\), \(M_{\rm halo}\), and \(V_{\rm max}\). We also report the scatter in the reconstructed \(\log(M_{\star}/M_{\odot})\).
### Cosmic substructure in projection

We also construct cosmic graphs in projection, i.e. projected coordinates \(x_{1}\) and \(x_{2}\), and radial velocity \(v_{3}\), instead of the full phase space information (see Appendix B). This \(2d\) GNN model achieves RMSE = \(0.135\) dex scatter, which still exceeds the performance of the best RF estimator (see Table 1). Because the \(2d\) GNN encode projected large scale structure information, it outperforms the RF models that can only learn isolated subhalo information.

## 5 Discussion

We have presented a novel method for populating dark matter subhalos with galaxy stellar masses. Mathematical graphs combine individual halo properties and environmental parameters in an equivariant representation, resulting in robust predictions for both central and satellite galaxies. As shown in Table 1 and Figure 2, the cosmic graphs outperform random forests trained on \(V_{\rm max}\) and \(M_{\rm halo}\). For galaxies with \(\log(M_{\star}/M_{\odot})\geq 9\) and \(\log(M_{\rm halo}/M_{\odot})\geq 10\), we recover the logarithmic stellar mass to within a root mean squared error (RMSE) of \(0.129\) dex.

### Inductive biases of GNNs

We note that previous works have employed convolutional neural networks (CNNs) for painting stars onto dark matter maps (Zhang et al., 2019; Mohammad et al., 2022). Unlike abundance matching models and RFs, CNNs are able to represent local spatial information. However, CNNs and GNNs have different inductive biases: CNNs are well-suited for representing fields discretized onto a Cartesian grid, while GNNs are well-suited for representing objects and relationships between them. Galaxies have small sizes (\(\sim\)kpc) relative to their typical separations (\(\sim\)Mpc), and they interact with each other (and their surrounding media) through multiple physical mechanisms (e.g., gravitational attraction, tides, ram pressure, etc). Therefore, cosmic structures naturally conform to a graphical representation, motivating our use of GNNs in this work.

### Galaxy environments

We note that a GNN with no edges except self-loops would essentially model the galaxy-halo connection in isolation; all environmental information is contained and passed along the edges. However, if we remove self-loops from the GNN, then the GNN is still able to infer \(\log(M_{\star}/M_{\odot})\) to within RMSE \(\sim\) 0.145 dex. A GNN without self-loops must estimate galaxy stellar mass _solely_ from neighboring halo information, which demonstrates that galaxy environments are informative for modeling the galaxy-halo connection.

We find that the GNN with max-pooling aggregation function achieves 0.001 dex lower RMSE than a GNN with sum-pooling. This result suggests that the GNN selects the largest value for some combination of \(M_{\rm halo}\), \(V_{\rm max}\), and distance to neighboring subhalos in order to best make predictions. We can speculatively interpret this as evidence that the largest and most nearby subhalo is most informative to a GNN. The largest subhalo might dominate environmental effects (e.g. tides and ram pressure) and control a given subhalo's stellar mass. Meanwhile, the summed information should capture _all_ of the forces, and we expect it to be more robust or transferable across domains. This interpretation

\begin{table}
\begin{tabular}{l c c c c c c c} \hline \hline
**Model** & **RMSE** & **MAE** & **NMAD** & **Pearson \(\rho\)** & \(R^{2}\) & **Bias** & **Outlier fraction** \\  & (dex) & (dex) & (dex) & & & (\(10^{-3}\) dex) & (\%) \\ \hline AM - \(M_{\rm halo}\) & 0.424 & 0.327 & 0.323 & 0.736 & 0.472 & 0.1 & 3.73 \\ AM - \(V_{\rm max}\) & 0.173 & 0.150 & 0.132 & 0.956 & 0.912 & **0.0** & 1.91 \\ SHAM - \(M_{\rm halo}\) & 0.332 & 0.231 & 0.235 & 0.838 & 0.677 & 0.1 & 6.20 \\ SHAM - \(V_{\rm max}\) & 0.151 & 0.133 & 0.115 & 0.966 & 0.933 & **0.0** & 1.75 \\ RF - \(M_{\rm halo}\) & 0.366 & 0.308 & 0.277 & 0.780 & 0.606 & \(\mathbf{-0.0\pm 0.1}\) & \(2.53\pm 0.01\) \\ RF - \(V_{\rm max}\) & 0.197 & 0.177 & 0.152 & 0.942 & 0.886 & \(-0.3\pm 0.0\) & \(1.44\pm 0.01\) \\ RF - \(M_{\rm halo}+V_{\rm max}\) & 0.148 & 0.135 & 0.114 & 0.967 & 0.936 & \(0.3\pm 0.0\) & \(1.31\pm 0.00\) \\ GNN (\(2d\) projection) & 0.135 & 0.131 & 0.106 & 0.973 & 0.946 & \(-3.9\pm 2.2\) & \(0.68\pm 0.01\) \\
**GNN (\(3d\))** & **0.129** & **0.125** & **0.102** & **0.975** & **0.951** & \(0.8\pm 0.6\) & \(\mathbf{0.68\pm 0.00}\) \\ \hline GNN (\(3d\)) - centrals & 0.123 & 0.119 & 0.097 & 0.979 & 0.959 & \(4.6\pm 0.7\) & \(0.67\pm 0.01\) \\ GNN (\(3d\)) - satellites & 0.138 & 0.136 & 0.109 & 0.968 & 0.936 & \(-5.0\pm 0.6\) & \(0.58\pm 0.01\) \\ \hline \hline \end{tabular}
\end{table}
Table 1: Cross-validation metrics for the AM, SHAM, RF, and GNN models discussed in the text. The GNN trained using positions, velocities, \(M_{\rm halo}\), and \(V_{\rm max}\) achieves the best metrics (shown in bold) in nearly every category. The last two rows report metrics for the \(3d\) GNN model, except that only central and satellite subhalos are selected from the cross-validation set. We note that the AM/SHAM models are trained and evaluated on the same data set. For RF and GNN models, we repeat the entire training and cross-validation experiment three times; the scatter is too small to be shown in the displayed significant figures for all columns except the bias and outlier fraction.
requires addition testing and an exhaustive hyperparameter search over GNN architecture and optimization procedures, which we aim to do in a follow-up work.3

Footnote 3: The linking length is a particularly important hyperparameter. In our preliminary tests, we have found 5 Mpc to give good results.

### Applications to observations

The strong performance of \(2d\) GNNs (SS4.3) is promising for facilitating comparisons to observations beyond the Local Group, where we can only reliably measure projected positions and line-of-sight velocities rather than full phase space information. Our method can be used to quickly estimate galaxy properties of constrained \(N\)-body (McAlpine et al., 2022) and Gpc-scale \(N\)-body simulated volumes (Garrison et al., 2018; Maksimova et al., 2021) for comparison with wide-area galaxy surveys in the low-redshift Universe (Ruiz-Macias et al., 2021; Carlsten et al., 2022; Darragh-Ford et al., 2022; Driver et al., 2022; Wu et al., 2022).

### Limitations and caveats

While we have shown that the GNN outperforms other methods, this demonstration does not definitively prove that GNNs are exploiting environmental information. Indeed, we have used a linking length of 5 Mpc, but this hyperparameter may be suboptimal and should be tuned. It is also possible that intrinsic scatter imposes a RMSE floor (i.e., due to the "butterfly effect" in cosmological simulations Genel et al., 2019), although GNN results using merger trees have shown that galaxy stellar mass can be recovered to even lower scatter (Jespersen et al., 2022). Finally, it may be that merger history is more important than environmental information, and that the clustering information learned by a GNN only incrementally improves performance relative to other approaches.

Our results will depend on choice of halo finder, i.e. if we were to use an alternative to the SUBFIND algorithm (e.g. ROCKSTAR; Behroozi et al. 2013). We have not tested our results using different halo finding tools, and it is unclear whether a GNN trained using one halo finder catalog will properly generalize to another catalog produced by a different halo finder. We also note that our results, while promising, must be tested on dark matter only simulations with halo catalogs matched to the hydrodynamic simulation catalogs before we can rely on GNNs to paint galaxies onto dark matter subhalos.

Additionally, domain adaptation will likely be needed to ensure simulated results can transfer to other simulations (e.g., while varying cosmological parameters; Villaescusa-Navarro et al. 2021) or to observations (e.g, Ciprijanovic et al., 2023). As a preliminary test, we repeat our experiment by training on TNG300 and validating on TNG50 data, and vice versa; in both cases the results are poor (\(>0.2\) dex). However, by training on a subset both simulations, we can recover \(\log(M_{\star}/M_{\odot})\) to \(\sim 0.13\) dex for TNG300 and \(\sim 0.14\) dex for TNG50 (Nelson et al., 2019; 2). This test suggests that cross-domain applications, such as transferring GNN results from simulations to observations, will necessitate some form of domain adaptation.

## Software and Data

Our code is completely public on Github: [https://github.com/jwupphysics/halo-gnns/tree/halos-to-stars](https://github.com/jwupphysics/halo-gnns/tree/halos-to-stars). We have used the following software and tools: numpy (Harris et al., 2020), matplotlib(Hunter, 2007), pandas(Wes McKinney, 2010), pytorch(Paszke et al., 2019), and pytorch-geometric (Fey and Lenssen, 2019).

We only use public simulation data from Illustris, which can be downloaded from [https://www.tng-project.org/data/](https://www.tng-project.org/data/).

## Acknowledgments

JFW and CKJ thank Peter Behroozi, Haley Bowden, Francisco Villaescua-Navarro, Tjitske Starkenburg, and Risa Wechsler for valuable discussions that sharpened this work. We also thank the two anonymous reviewers who provided excellent comments and suggestions that improved this manuscript. This research has made use of NASA's Astrophysics Data System Bibliographic Services. The authors are grateful to the Kavli Institute for Theoretical Physics "Building a Physical Understanding of Galaxy Evolution with Data-driven Astronomy" program, where this work began. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958.
