Explaining dark matter halo density profiles with neural networks

We use explainable neural networks to connect the evolutionary history of dark matter halos with their density profiles. The network captures independent factors of variation in the density profiles within a low-dimensional representation, which we physically interpret using mutual information. Without any prior knowledge of the halos' evolution, the network recovers the known relation between the early time assembly and the inner profile, and discovers that the profile beyond the virial radius is described by a single parameter capturing the most recent mass accretion rate. The results illustrate the potential for machine-assisted scientific discovery in complicated astrophysical datasets.

Introduction -In the modern picture of structure formation, galaxies form at the center of extended, overdense 'halos' of dark matter, which originate from small fluctuations in the density of matter in the early Universe and undergo highly non-linear dynamical processes throughout their evolution [1][2][3][4][5].The history of a halo determines its final structure, commonly parametrized by the spherically-averaged radial density profile.Halo density profiles are not only key ingredients of the galaxy-halo connection in cosmological analyses and of direct and indirect dark matter searches; they are also powerful observational testbeds of fundamental physics.This is because their shape, from the inner core to the outskirts in the proximity of the splashback radius, is sensitive to the nature of dark matter and modifications to gravity [6,7].
Observationally, it has recently become possible to measure weak lensing and 3D density profiles through a combination of multi-wavelength data; upcoming data from Euclid, Rubin, and DESI will provide even more detailed measurements of the density profiles of halos from clusters to dwarf galaxies [8][9][10].Achieving the potential impact of these measurements requires determining the physical effects that control the shape of the density profiles.However, current theoretical models are limited to empirical fitting functions such as the Navarro-Frenk-White (NFW) [11] and Einasto [12] profiles; these do not explain the physical origin of the profiles' universal shape seen in numerical simulations [13,14].Understanding the connection between the formation history and the density profile also offers the possibility of using observational constraints on the latter to estimate the halos' mass accretion rate.This will yield valuable constraints on galaxy formation [15][16][17][18] as well as extensions to the cosmological model.
In this work, we use an explainable AI framework to connect the formation histories of dark matter halos to their density profiles.Our goal differs from typical uses of machine learning in cosmology, such as emulating the output of computationally expensive simulations [19,20] or accelerating the estimation of cosmological parameters from data [21][22][23][24].In Ref. [25], we used an interpretable deep learning framework to build a new model for the spherically-averaged halo density profiles, generalizing over existing empirical fitting functions.The framework, which we denoted an interpretable variational encoder (IVE), was trained to capture all the information used by the neural network to predict the profile, given the 3D density field around the halo center, within a compact, low-dimensional latent representation.We require the representation to be disentangled i.e., each latent component captures different, independent factors of variation in the profiles; the latent representation is equivalent to the profiles' degrees of freedom.We found that three components are required (and sufficient) for modeling the profiles out to the halo outskirts: these three components describe respectively the normalization of the profile, and its shape within and beyond the virial radius.
In this Letter, we turn to the physical interpretation of the learnt IVE latent representation to investigate how halo density profiles are determined from the halos' formation histories.Although the network was trained only on the presentday density field, we explore whether the latent parameters carry memory of the evolution history of the halos.We measure the information encoded within each latent about the halos' evolution history using the information-theoretic measure of mutual information (MI).By this metric, the IVE representation and the NFW parametrization similarly highlight a dependence of the profile on physical accretion history.However the IVE additionally allows us to measure the connection between a halo's recent evolution history and the density in its far outskirts, something that the NFW profile does not capture.
Background -We begin by briefly reviewing the current understanding of the physics of halo density profiles.The NFW profile is the most widely used fitting function for the halo density profile.It is given by where r s and ρ s are the scale radius and characteristic density, respectively.The scale radius is often re-written in terms of a concentration parameter c ≡ r 200m /r s , so that the NFW profile depends on the virial radius r 200m and concentration c.The virial radius r 200m is typically adopted as a proxy for the halo boundary, and defined as the radius which contains a mean density that is 200 times the mean density of the Universe.High-resolution simulations have revealed this functional form to be 'universal': it provides a good fit to stacked profiles of halos for a large range of halo masses [14,26], for several different cosmological models [27][28][29][30][31], and even in the absence of hierarchical growth [32][33][34].This suggests that universal density profiles are a generic feature that arises from collisionless gravitational collapse.
Despite the lack of a first-principles explanation for the selfsimilarity of halo density profiles, some insights have been gained from studying the correlation between the NFW concentration and summary statistics of the halo evolution process.Mass, concentration and halo formation time all correlate: on average, low-mass halos assemble earlier and have higher characteristic densities (or concentration), reflecting the larger background density at earlier times [11,26,30,[35][36][37][38].This description can explain the qualitative trend of the mean concentration as a function of halo mass, but not the large residual scatter in concentration seen in simulations [36,39].It is also limited to the simplest summary statistic of the halo evolution history, i.e. the halo formation time.The scatter in concentration at fixed mass has been shown to be at least in part connected to merger events during the halo assembly process [39][40][41].Further work has suggested that the self-similarity of halos may be related to the self-similarity of the halo mass assembly history [26], although this has only been validated on stacked profiles of well-behaved,'relaxed' halos.
The situation worsens when modeling profiles beyond the virial radius: the halo outskirts strongly deviate from the NFW form due to the presence of the splashback radius, where particles reach the apocenter of their first orbit.Recent work has focused on modeling the location of the splashback radius, finding that it is sensitive to the late-time mass accretion rate [42][43][44][45][46][47].Modeling the full shape of the outer profile remains a difficult task due to its intrinsically non-equilibrium nature, leading to a reliance on multi-parameter fitting functions with little physical explainability [48,49].
Deep learning model -The IVE architecture used in this work has two main components: the encoder, mapping the 3D density field to a low-dimensional latent representation, and the decoder, mapping the latent representation and the query radius log(r) to the output profile log[ρ(r)].By design, all the information used by the model to predict the density profiles is captured within the latent representation.An illustration of the model is shown in the top-half of Fig. 1.The encoder is a 3D convolutional neural network (CNN) with parameters ϕ that maps the inputs x to a multivariate distribution in the latent space p ϕ (z|x).We choose the latent representation to be a set of independent Gaussians, , where L is the dimensionality of the latent space; under this assumption, the encoder maps the inputs x to the vectors µ = µ i , .., µ L and σ = σ i , .., σ L .The decoder of the IVE consists of another neural network model with parameters θ that maps a sampled latent vector z ∼ p ϕ (z|x) and a value of the query log(r) to a single predicted estimate for log[ρ pred (r)].
A crucial aspect of the IVE that makes the latent space interpretable is that it is disentangled: independent factors of variation in the density profiles are captured by different, inde- pendent latents.This is achieved through the design of a loss function that minimizes the mean squared error between predicted and ground-truth profiles, while simultaneously maximizing the degree of independence between the latent variables by encouraging those to be as close as possible to independent Gaussians of mean 0 and variance 1 [50].More details on the encoder and decoder architectures and the loss function are presented in Ref. [25].
Methods -We generated the training data from four darkmatter-only N -body simulations produced with GADGET-4 [51], each containing 512 3 particles in a (50 Mpc h −1 ) 3 box.We trained two IVE models for different tasks: one (IVE virial ) was trained to model the density profile up to the halo virial radius r 200m , and the second (IVE infall ) was trained to model profiles beyond the halo boundary out to 2 r 200m .The former is used for direct comparison with the NFW profile, which is also designed to model the profile out to the virial radius, and the latter is used to investigate the less studied halo outer profile.The innermost radius of the profiles we consider is r min = 3 ϵ, where ϵ is the gravitational softening of the simulation; this choice ensures that we can robustly trust the inner profile.The inputs are given by the 3D density field within a N = 131 3 sub-box of size L sub−box = 0.4 Mpc h −1 for the IVE virial model, and of size L sub−box = 0.6 Mpc h −1 for the IVE infall one.We considered halos with log 10 (M/M ⊙ ) ∈ [11,13], but for the IVE infall model, we further restricted our analysis to halos with r 200m ≤ 150 kpc h −1 .These cuts yielded ∼17,000 (13,000) halos for training the IVE virial (IVE infall ) model.Further discussion on the training data of the IVE virial and IVE infall models is presented in Ref. [25].To compare the IVE virial results with NFW, we fitted the NFW formula in Eq. (1) to each halo's density profile using least-squares minimization, and recovered the best-fitting parameters r s and ρ s .The concentration was then derived using c =r 200m /r s .A description of the simulations used for training and testing the IVE models can be found in the Supplemental Material.The first step of the analysis was to verify that the IVE models learn to predict the density profiles at ∼ 5% accuracy, comparable to the accuracy of NFW fits (see Supplemental Material).Crucially, only the z = 0 snapshots were used for training the IVE models to construct the latent representations mapping the 3D density field to dark matter halo profilesi.e., the model had no access to the merger histories of the halos during training.The resulting disentangled latent space directly corresponds to the underlying degrees of freedom in the halo density profiles.Following recent works [25,[52][53][54][55], we then used the MI to (i) quantify the information captured by the latent space about the halo density profiles and (ii) connect the IVE latents to the halo's evolution history, showing how the latter determines the present-day density profile.
The MI was estimated using GMM-MI [55], which performs density estimation using Gaussian mixtures and provides MI uncertainties through bootstrap.Background details on MI are provided in the Supplemental Material.We first measured the MI between each latent and the density profile ρ(r); this allows us to directly link each latent to a degree of freedom in the profile that affects its shape over a certain radial range.We then measured the MI between each latent and the mass assembly history of each halo.This in turn allowed us to connect each degree of freedom describing the density profile of the halo directly to characteristics of the halos' evolution which determines that component.
Results -Figure 2 quantifies the information contained within the latents of the IVE infall (top panel) and the IVE virial (bottom panel) models about the ground-truth density pro- files1 .We show the MI between each latent parameter and the ground-truth profiles, which we denote as MI ρtrue(r) .The three latents discovered by the IVE infall describe (i) the normalization of the profile, which dominates the variation in the profiles out ∼ r 200m /2, (ii) the shape of the inner profile, which becomes informative on radial scales approaching r 200m , and (iii) the shape of the outer profile beyond r 200m .
The first two are analogous to the two NFW parameters, mass and concentration, respectively.A closer comparison between the inner shape latent of the IVE virial model and concentration (bottom panel of Fig. 2) shows that both parameters carry information about the density in the core and on radial scales close to r 200m2 .This bimodality is due to a compensation effect between the density in the inner region and that close to the virial boundary: at fixed normalization, halos with denser cores become less dense in the outskirts, and vice versa.The MI ρtrue(r) of the inner latent is shifted towards larger radii compared to that of concentration, suggesting that the former is sensitive to variations in the shape of the profile on larger radial scales than the latter; this distinction will become relevant when physically interpreting the latent and comparing it to concentration.
We now move on to a physical interpretation of the latents in relation to characteristics of the halos' evolution histories.Recall that the network did not have access to this information during training.The interpretation of the normalization latent is straightforward: it captures the z = 0 mass of the halo, M 200m .Their MI is ∼ 2.07 ± 0.01 nats, implying a strong correlation between the two.This also matches expectations from the literature [11,12], as halo mass also controls the normalization in the NFW and Einasto fitting functions.To physically interpret the inner and outer shape latents, we measure their MI with two quantities that describe the assembly history of the halos over cosmic time.The first is the mass accretion history, M 200m (z)/M 200m (z = 0), which describes the evolution of the halo mass as a function of time M 200m (z) normalized to the present-day halo mass M 200m (z = 0).The second is the mass accretion rate Γ(t) ≡ ∆ ln M 200m (a)/∆ ln a [46], which describes the rate of change in halo mass with respect to the scale factor a(t).The value of the accretion rate depends on the time interval used to compute the change in mass and scale factor; we compute Γ(t) by taking the finite difference of the halo masses at each consecutive timestep in the simulation.
Figure 3 shows the MI between the latents and mass accretion history (MI M (z) ; top row) and that between the latents and the mass accretion rate (MI dM (z)/dz ; bottom row).We first focus on the inner shape latent, which we compare to the NFW concentration.The MI M (z) of the inner shape latent increases with time during the early formation period, peaks at z ∼ 1, and declines rapidly towards z = 0; recall that this is the MI with the mass assembly history normalized to the present-day halo mass.This result reveals that the inner shape latent is sensitive to the early assembly history of halos.The MI dM (z)/dz of the same latent reveals that the latter is also sensitive to the later time mass accretion rate.This dual dependence explains the bimodal shape of the MI between the inner latent and the profile (Fig. 2, bottom panel): the early assembly phase determines the shape of the profile in the innermost region of the halo, while the later time mass accretion rate determines the shape of the profile close to the virial radius.We further validate this interpretation in the Supplemental Material.
The NFW concentration shows a similar picture to the inner shape latent.However, its MI M (z) peaks at earlier times (z ∼ 0.55) compared to the inner shape latent.This implies that the inner shape latent carries information about the build up of mass onto the halo over a longer period of time than concentration, which therefore affects the inner halo structure (and the profile) out to larger scales.The sensitivity of the inner latent to later times/larger scales in the profile explains why the inner latent MI ρtrue(r) is shifted towards larger radial scales than that of concentration (Fig. 2).Moreover, the absolute magnitude of the concentration MI M (z) is higher than that of the inner shape latent; this is because the closer to the halo core the stronger the correlation with the early assembly history due to halos accreting mass 'inside-out'.As a result, concentration, which is sensitive to the profile on smaller r than the latent, has a higher MI with the early assembly history than the latent.Finally, the NFW concentration is related to the later time mass accretion rate in a similar way to the inner latent.
Figure 3 shows that the outer shape latent shares information about the mass accretion rate over the past ∼ 5 Gyr, since this is the period over which the MI roughly doubles (see bottom panel).This timescale corresponds to the halo dynamical time, t dyn ≡ 2 × r 200m /v 200m , defined as the time it takes for material to cross the halo at a typical virial velocity v 200m = GM 200m /r 200m .This suggests that the outer profile is primarily determined by the infall of dynamically unrelaxed material within the last dynamical time, which has not yet virialized within the halo.
Discussion -Our results show that the IVE framework has extracted a direct connection between the assembly history of cold dark matter halos and their density profiles, without having access to explicit information about the time evolution of the halos during training.This has deep implications for understanding the origin of universality in dark matter halos; the universality in the profiles, captured by three degrees of freedom alone, may originate from a universality in the halo assembly histories themselves, since the latents contain comparable amounts of information about both quantities.
Previous work [26] found a resemblance between the shape of the average mass accretion history, expressed in terms of the critical density of the Universe, and the average enclosed mass profile, expressed in terms of its enclosed density, for a selected set of 'well-behaved' halos of similar mass.In the halo outskirts, the profile has been linked to the dynamical accretion history of the halos primarily through the relation between the splashback radius and the mass accretion rate [43,48]; existing models make use of multi-parameter fitting functions to capture the dynamical impact on the outer profile [49].
By contrast, within the IVE framework, the connection between the density profiles and the entire mass accretion history or mass accretion rate is clearly elucidated through MI.This result was obtained using all halos in the simulations, without requiring a curated sample of well-behaved halos.The IVE rediscovers the known correlation between the inner profile and halo formation time [35,36]; it then additionally demonstrates that the complexity of the dynamical, infalling material is encoded in only a single degree of freedom that captures the recent mass accretion rate.In future work, we will use the connection between assembly history, latents, and density profile captured by the IVE framework to build a model that can determine mass accretion histories from density profiles.
In future work we will explore extensions to this work using hydrodynamical simulations.We expect the same IVE framework to successfully disentangle the relevant factors in the baryonic case.Previous work found that baryons primarily impact the inner profile [56][57][58], and that the results can still be encoded with minimal modifications to the pure dark matter expectations [59,60].Conversely, the splashback radius in the halo outskirts remains unchanged when comparing hydrodynamical and dark-matter-only simulations [61].Thus, an IVE with the same dimensionality or a single additional dimension should suffice to account for the impact of baryonic physics on the halo profiles for the baryonic feedback models included in the training set simulations.
More broadly, our results represent progress towards enabling new machine-assisted scientific discoveries, going beyond artificial rediscovery of known physical laws [62][63][64].Our IVE approach towards this goal consisted of compressing the information within a dataset into a set of minimal ingredients which disentangles the independent factors of variation in the output (interpretability), and can be explained in terms of the physics it represents through MI (explainability).The approach shows promise for gaining insight into other emergent properties of the cosmic large-scale structure (e.g.void density profiles [65] and the halo mass function [66]), building physical explanations that are more accurate and complete than traditional methods have achieved.
Mutual information (MI) is a well-established informationtheoretic measure of the correlation between two random variables.In contrast to linear correlation measures such as the r-correlation, MI captures the full (linear and non-linear) dependence between two variables.In other words, the MI between two variables is zero if and only if these are statistically independent.Additionally, MI is invariant under invertible, non-linear transformations, e.g.translations, rotations and any transformation preserving the order of the original elements.
Mathematically, the MI between two continuous random variables x and y is given by MI(x, y) = p(x, y) log p(x, y) p(x)p(y) dxdy, (S.1) where p(x, y) is the joint probability density distribution between x and y, and p(x) and p(y) are their marginal distributions, respectively.The log is the natural logarithm, so that MI is computed in natural units of information (nats).
The main challenge in estimating MI is obtaining a density estimate for the joint probability distribution p(x, y) in Eq. (S.1).We make use of the publicly-available software GMM-MI [55], which performs density estimation using Gaussian mixtures and additionally provides MI uncertainties through bootstrap.

B. Simulations
We begin with a description of the simulations used for training the IVE models and subsequent interpretation of the learned representation.We ran four dark-matter-only N -body cosmological simulations using the publicly-available code GADGET-4 [51] assuming a Planck ΛCDM cosmological model [68].We evolved N = 512 3 dark matter particles in a (50 Mpc h −1 ) 3 box from z = 99 to z = 0.The four simulations are based on different realizations of the initial Gaussian random field, generated using GENETIC [69].Three simulations were used for training the machine learning model and one was set aside for validation and testing.
Dark matter halos were identified at z = 0 using the SUB-FIND halo finder [51,70], as done in Ref. [25].We restricted our analysis to halos within the mass range log 10 (M/M ⊙ ) ∈ [11,13], in order to fully resolve the inner profile of the lowest-mass halos and not be affected by small-number statistics at the high-mass end.To track the evolution history of the dark matter halos, we saved 91 snapshots of the simulations between z ∼ 7 and z = 0. We used the PYNBODY and TANGOS software packages [71] to construct the merger trees of every dark matter halo.TANGOS matches a halo with its successor in time based on the fraction of common particles between the two objects; the procedure is repeated for every timestep in the simulation, thus yielding halo merger trees from z ∼ 7 to z = 0.The merger trees were then used to track the mass of each halo's main progenitor over time.

C. Predictions for the halo density profiles
Figure S.1 shows the accuracy of the predictions of the IVE virial and IVE infall models.We show the mean and 90% confidence interval of the residuals log 10 [ρ predicted /ρ true ], in every radial bin of the profile used for testing.Since every radial bin contains a different value of r for different halos, we plot the residuals as a function of r eff defined as the median of the distribution of radius values within each bin.The grey band shows the residuals of the NFW fits for comparison.Note that the IVE results include uncertainties in the latent distributions, whereas the NFW fits do not include uncertainties as they were obtained through least-squares optimization.
The performance of the IVE models is consistent with that of the NFW profile, meaning that our model contains sufficient predictive accuracy to yield meaningful latent representations.

D. Further physical interpretation of the inner shape latent
We present a further investigation on the physical interpretation of the inner shape latent.As shown in Fig. 2, the MI between the inner shape latent and the ground-truth density profile has a bimodal shape.The MI first peaks at r 1 ∼ 0.1 r 200m and then again at r 2 ∼ 0.9 r 200m , meaning that the latent contains information about the shape of the profile in the inner region of the halo and close to the halo virial radius.When physically interpreting the latent, we found that the latent carries information about the early formation history of the halo and the later time mass accretion rate (Fig. 3).We now verify if the early formation history is responsible for the shape of the profile in the inner region, while the later time mass accretion rate is what determines that closer to the virial radius.To do so, we compute the MI between the ground-truth density at the first MI peak, ρ(r 1 ), and both the halo mass assembly history and mass accretion rate.We then repeat the calculation for the ground-truth density at the second MI peak, ρ(r 2 ).Fig. S.2 shows the MI between the two ground-truth densities with the mass assembly history, M (z), in the top panel 3 and that with the mass accretion rate, d ln M (a)/d ln a, in the bottom panel.The density in the inner region of the halo carries information about build up of mass from early times, but is largely uncorrelated with the mass accretion rate; the density close to the halo boundary is sensitive to the later time mass accretion rate, but does not depend on the early assembly phase.This confirms our physical interpretation of the bimodal shape of the MI between the latent and the profile in Fig. 2: the way the inner profile responds to the latent depends on the early formation history of the halos, while the way the profile close to r 200m responds to that same latent depends on the later time mass accretion rate.The MI between ρ(r 2 ) and the mass accretion rate shows an additional peak at z ∼ 0 that is not as pronounced in the MI between the inner shape latent and the mass accretion rate.This is because, as one approaches the virial radius, the profile also becomes sensitive to the most recent mass accretion rate, which is captured by the outer latent.
FIG.1.A neural network is trained to discover the underlying degrees of freedom in halo density profiles in the form of a latent representation, when presented with the full 3D density structure of a halo.We physically interpret the discovered representation by measuring the MI between the latent parameters and the assembly history of the halos.
FIG.3.The MI between the latent parameters and the mass accretion histories (denoted MI M (z) ; top row), and that between the latent parameters and the mass accretion rate (denoted MI dM (z)/dz ; bottom row).The inner shape latent and the NFW concentration carry memory of the early-time mass assembly history, as well as the later-time mass accretion rate.The outer shape latent carries information about the halos' most recent mass accretion rate over the past dynamical time (indicated by the arrow).
Mean and 90% confidence interval of the residuals log[ρ predicted /ρtrue] of the IVE virial and IVE infall models, as a function of r eff defined as the median radius in each bin.The grey band shows the NFW residuals.
ρ(r 1 ) ρ(r 2 ) FIG. S.2.The MI between the densities on two fixed radial scales, ρ(r1) and ρ(r2), and (i) the mass accretion histories (top row) and (ii) the mass accretion rate (bottom row).The two radial scales, r1 and r2, are the locations at which the MI between the ground-truth profile and the inner shape latent peaks (Fig 2, bottom panel).