Characterizing the Conditional Galaxy Property Distribution using Gaussian Mixture Models

Understanding the physical processes that govern the formation and evolution of galaxies over cosmic time is one of the most fundamental and important problems in astrophysics and cosmology. Numerous astronomical surveys across a wide range of wavelengths have been carried out to observe galaxies, from nearby targets to those probing the earliest epochs (e.g. York et al., 2000; Koekemoer et al., 2011). However, galaxy surveys are limited by the trade-off between area and depth since deeper and fainter targets require more time to achieve a detectable signal-to-noise ratio. Instead of resolving and observing individual targets, a relatively recent technique called line-intensity mapping (LIM) (e.g. Kovetz et al., 2017) has been proposed to measure the integrated emission lines from all sources, bright and faint. Similar to traditional galaxy surveys, LIM not only provides a promising tracer of the large-scale structure (LSS) of the total matter distribution, but also contains a wealth of information about astrophysical properties of galaxies (Schaan & White, 2021a).

One of the potentially powerful aspects of LIM is that there are now numerous experiments planned or online that probe different lines, which trace different gas phases, at different redshifts (see Bernal & Kovetz, 2022, for a recent summary). For example, some of the most commonly considered lines include CO, which is a tracer of dense molecular gas in the interstellar medium (ISM) of galaxies, atomic and ionic fine structure lines such as [CII], which is primarily emitted by photo-dissociation regions in the outer layers of molecular clouds, 21-cm emission from atomic gas in the ISM and intergalactic medium (IGM), and lines arising from ionized gas such as Lyman-$\alpha$ and H-$\alpha$. These are the brightest and most easily detected lines, but there are many others that can provide interesting information. It is clear that combining multi-tracer LIM surveys with traditional galaxy surveys will be the key to extracting the maximum amount of science from these experiments.

In order to do this, it is essential to have a comprehensive modeling framework. The most commonly used framework is halo-based phenomenological models (see Bernal & Kovetz, 2022, for a recent review). In this approach, scaling relations between halo mass and redshift and “intrinsic” physical properties such as stellar mass and star formation rate are derived using abundance matching. In abundance matching, observational estimates of stellar mass functions and SFR are used along with predictions for the number density of dark matter halos as a function of mass to derive constraints on the mapping between these quantities (see Wechsler & Tinker, 2018, for a review). Then, scaling relations between e.g. SFR and a line luminosity of interest are derived using observations of individually detected sources (e.g. Silva et al., 2015; Li et al., 2016; Padmanabhan, 2018, 2019; Sun et al., 2019). This approach has the advantage of being highly computationally efficient, and of being tied to observations.

However, there are a number of drawbacks and limitations to this approach. First, these empirical fitting functions usually lack physical motivation, so they are able to provide only limited insights about the physical processes underlying galaxy formation. Second, LIM detections are sensitive to galaxies that are too faint to detect individually, so the scaling relations have to be extrapolated beyond the observational range. Third, these empirical approaches are not able to self-consistently capture the covariances between different tracers, or between the intrinsic properties (stellar mass, SFR) and various emission lines. It is important to keep in mind that the estimates of stellar mass and SFR that are used in abundance matching models are not directly observable quantities, but are derived from UV-optical-NIR observations and carry significant uncertainties. The prior used to obtain these estimates can have a strong impact on the results (Leja et al., 2022). It is critical to be able to properly fold in the prior, covariances, and uncertainties on these quantities, which lie at the foundation of abundance matching approaches.

Numerical cosmological simulations that incorporate baryonic astrophysical processes are a potentially powerful tool for exploring the physics of galaxy formation and evolution, as well as the connection with observable line emissions. The most fundamental approach is based on explicitly simulating the hydrodynamics and thermodynamics of baryonic material evolving within a $\Lambda$CDM cosmological framework for structure formation (Somerville & Davé, 2015a; Naab & Ostriker, 2017). These numerical hydrodynamic simulations provide direct predictions for the distribution and state of nearly all of the relevant components of galaxies and the IGM (stars, gas at different densities and temperatures, various chemical elements). They can be coupled with modeling tools to provide predictions of many of the emission lines of interest for LIM experiments (Leung et al., 2020; Olsen et al., 2021). However, due to its computational expense, this technique suffers from a fairly severe trade-off between volume and resolution. Simulations that resolve the multi-phase ISM can typically be run for only a few galaxies (e.g. Pallottini et al., 2022). Simulations of cosmological volumes must adopt “sub-grid” recipes to treat most of the critical processes of galaxy formation, including star formation, stellar feedback, and black hole formation and feedback. And even the largest computationally feasible cosmological hydrodynamic simulations fall short of the volumes that are already routinely probed by galaxy and LIM surveys. For example, MillenniumTNG has a volume of (500 Mpc/$h$)³, compared to (1760 Mpc/$h$)³ for the Baryon Oscillations Spectroscopic Survey of $z<0.7$ luminous red galaxies alone (Alam et al., 2017). The currently adopted sub-grid models are highly uncertain, and approaches adopted by different groups can lead to dramatically different predictions for gas properties. It is infeasible to thoroughly explore the parameter space of the sub-grid models for significant volumes using a direct numerical approach.

An alternative approach to physics-based models is semi-analytic models (SAM), which combine simplified treatments of physical processes with cosmological halo merger histories. These models produce similar predictions to those of numerical hydrodynamic simulations for global galaxy properties, and include a similar suite of physical processes, but at a fraction of the computational cost. However, unlike numerical hydrodynamics simulations they do not provide detailed information about the internal structure of galaxies or the spatial distribution of gas within or around them. Popping et al. (2014) and Popping et al. (2019) coupled semi-analytic models of galaxy formation with the DESPOTIC code (Krumholz, 2014), which computes line luminosities for optically thick clouds using a one-zone model, to make predictions for a broad suite of emission lines, including multiple transitions in the CO ladder, [CI] and [CII]. They validated these models by comparing the predictions with an extensive suite of observations of individually detected galaxies across a broad range of redshifts. Yang et al. (2021) used lightcones populated with this approach to create mock maps (including interlopers, cosmic infrared background, and Galactic foreground) for exemplar LIM experiments in CO and [CI]. This work showed that SAMs coupled with the appropriate line emission models in post-processing are a promising alternative tool for LIM forecasts and interpretation.

Ideally, we would like to be able to use Bayesian inference to constrain the posteriors of both cosmological and astrophysical parameters from a combined suite of multi-tracer LIM and galaxy survey observables. When using a forward-modeling or simulation-based inference approach on a large population of galaxies or a cosmological volume (Henriques et al., 2009, 2013; Lu et al., 2011, 2012, 2014; Forbes et al., 2019), depending on the sampling approach adopted, this requires thousands to tens of thousands of executions of the simulation. Most currently available SAM codes are not fast enough for this to be practical, particularly if large volumes need to be simulated.

An intermediate approach is to develop models that link halo properties to galaxy intrinsic and observable properties using physics-based simulations, such as hydrodynamic simulations and SAMs, instead of empirical methods such as abundance matching. For example, Yang et al. (2022) used the SAM-based predictions described above to derive and present fitting functions between halo mass and all major predicted lines over a broad range of redshifts, and showed that when these fitting functions were combined with dark matter halo properties as in the standard empirical halo model approach, agreement of 5-10% was achieved for common LIM summary statistics such as the power spectrum and intensity.

In this work, we introduce a new approach that we believe can have very general applicability to forecasting and interpreting a wide variety of galaxy and LIM survey data. We introduce the conditional galaxy property distribution (CGPD) $\Phi\left(\bm{p}|M,z\right)$, which gives the number density of galaxies in the property space $\bm{p}$, conditioned on the halo mass $M$ and redshift $z$. Whereas traditional halo models are typically mappings from a scalar quantity (halo mass) to another scalar (e.g. stellar mass or SFR), we can think of this as an extension of the halo model to map to a higher dimensional suite of galaxy properties. In this work, $\bm{p}$ has dimension five and includes the galaxy stellar mass, molecular gas mass, galaxy radius, gas phase metallicity and star formation rate (SFR), which are important for predicting the line luminosities of interest (here, we have in mind mainly CO and [CII]). However, in principle $\bm{p}$ could contain any galaxy property of interest and have an arbitrary number of dimensions. With this set of galaxy properties, CGPD also encapsulates many popular and important galaxy distribution functions and scaling relations such as the stellar mass function and SFR main sequence.

In this paper, we explore the form and effective dimensionality of the CGPD for this chosen set of key galaxy properties, using two different physics based cosmological simulations of galaxy formation: the “Santa Cruz” semi-analytic model (Somerville & Primack, 1999; Somerville et al., 2008, 2012; Porter et al., 2014; Popping et al., 2014; Somerville et al., 2015), and the IllustrisTNG hydrodynamic simulation suite (Springel et al., 2018; Pillepich et al., 2018; Marinacci et al., 2018; Nelson et al., 2018; Naiman et al., 2018). A key point is that it is well known that galaxy properties are highly covariant, giving rise to the many familiar galaxy scaling relations such as the stellar mass vs. SFR relation (SF sequence), mass metallicity relation, etc. Exploring higher dimensional manifolds within galaxy property space (as hinted at by, for example, “fundamental plane” type relations such as the fundamental mass metallicity relation) will in itself provide interesting insights into galaxy formation, as the physics that shapes galaxy properties is encoded in these relationships. Moreover, it will help us to develop tools to characterize and model these high-dimensional relationships in order to fully capture the covariances between galaxy properties. In this work, we adopt a Gaussian mixture model (GMM) based approach as a generative density estimation model to describe the CGPD.

Combined with a model that can predict emission line luminosities from these galaxy properties, this parametric model could potentially be used in a statistical inference pipeline, either using “traditional” methods such as Markov Chain Monte Carlo (MCMC) sampling combined with an assumed form for the likelihood, or more modern implicit likelihood inference methods (Alsing et al., 2018; Cranmer et al., 2020). With the formulation presented here, the posterior of the parameters of the GMM (and hence the CGPD) can be constrained based on a chosen set of LIM observables. The resulting CGPD can then be projected and integrated to obtain estimates for galaxy population summary statistics of common interest, such as the SFR or cold gas distribution functions or cosmic SFR or cold gas density, at a series of redshifts. If desired, additional observational constraints from galaxy surveys, such as the stellar mass function, can simultaneously be used as constraints in the inference procedure. In Fig. 1, we provide a graphical representation of how the CGPD model that we develop in this work fits into this broader framework, which we plan to build up in future works.

The paper is organized as follows. In Section 2, we introduce the CGPD and how it is related to more familiar summary statistics that are commonly used to characterize galaxy populations. We introduce the simulated galaxy catalogs with physical properties of interest in Section 3. Based on these simulations, we then present the application of the GMM in characterizing the CGPD in Section 4. We discuss our results in Section 5, and conclude in Section 6. Throughout this work, we adopt a flat $\Lambda$CDM fiducial cosmology with Planck 2015 parameters: $\Omega_{m}=0.3089$, $\Omega_{\Lambda}=0.6911$, $\Omega_{b}=0.0486$, $h=0.6774$, $\sigma_{8}=0.8159$, $n_{s}=0.9667$ (Planck Collaboration et al., 2016), which are consistent with the cosmology assumed in the simulations.

For a halo of mass $M$ at redshift $z$, the CGPD $\Phi\left(\bm{p}|M,z\right)$ is defined as the number density distribution of galaxies in the property space $\bm{p}$. Integrating $\Phi\left(\bm{p}|M,z\right)$ over $\bm{p}$ yields the expected number of galaxies for such a halo. For simplicity, in this work we write the CGPD as a function of $\bm{p}$, although in practice, we perform the fitting of the CGPD in $\tilde{\bm{p}}=\log\bm{p}$ space ¹¹1In this work, we denote a base 10 logarithm with $\log$ and natural logarithm with $\ln$., i.e. $\tilde{\Phi}(\tilde{\bm{p}}|M,z)$. The CGPD may be represented in either of the two spaces for convenience. The transformation of the CGPD between these two spaces can be straightforwardly given by

where $\left|\partial\,\tilde{\bm{p}}(\bm{p})/\partial\,\bm{p}\right|$ is the determinant of the Jacobian matrix. Multiplying the differential volume in $\bm{p}$ space on both sides gives

which is more convenient for the purpose of integration. The galaxy properties we consider are summarized in Table 1.

These are the properties that are the most important for the line luminosities of interest. Note that $\Phi\left(\bm{p}|M,z\right)$ is related to the conditional luminosity function (CLF) $\Phi(\bm{L}|M,z)$ introduced in Schaan & White (2021b) by the line luminosity relation (LLR) $L_{i}(\bm{p},z)$, which gives the luminosity of line $i$ for a galaxy with properties $\bm{p}$ at redshift $z$. The LLR is assumed to be independent of the host halo mass $M$. The modeling of the LLR using SAM simulated galaxies coupled with spectral synthesis models is being developed in a companion work (Breysse et al., 2023).

The CGPD must be coupled with the dark matter halo model in order to connect the cosmological dark matter field with observable galaxies. For example, we obtain the global distribution function for each element of the property vector $\bm{p}$ by performing the integral:

where $n_{\rm halo}(M,z)$ is the halo mass function (HMF) at redshift $z$. In this work, we use the HMF from Tinker et al. (2008) with the implementation in the hmf ²²2https://github.com/halomod/hmf code (Murray et al., 2013). This general galaxy distribution function $\phi(\bm{p}|z)$ gives the comoving volume number density of galaxies in the property space $\bm{p}$ at redshift $z$. Integrating $\phi(\bm{p}|z)$ over $\bm{p}$ yields the total number of galaxies per unit comoving volume. As shown below, by integrating over certain dimensions, $\phi(\bm{p}|z)$ encapsulates many observables of wide interest in galaxy and LIM surveys, such as the cosmic densities of galaxy properties, the galaxy abundance as a function of a certain property, and the scaling relations between two (or more) galaxy properties.

In this section we briefly introduce the physics-based simulations, which will be used to study and motivate the functional form and modeling approach for characterizing the CGPD. These include the “Santa Cruz” semi-analytic models and hydrodynamic simulations from the IllustrisTNG project. Note that comparing and discussing the differences between the simulated results from the SAM and TNG is outside of the scope of this work, and a recent discussion can be found in Gabrielpillai et al. (2022). It is important to consider different models in order to test the robustness and flexibility of the GMM in describing the CGPD.

Based on the simulated data described above, we now propose and test analytic models that could potentially properly characterize the CGPD.

There are various parametric and nonparametric methods for describing the distribution of point data, also known as density estimation (Ivezić et al., 2019). Histograms might be the most straightforward nonparametric method, which approximate the function value as constants within discrete bins. However, the total number of bins increases exponentially with the dimensionality of the data. Treating each bin as an independent variable makes it infeasible for inference, and also this method does not give any insights into the structure of the distribution. This remains an issue even for more advanced nonparametric models like Gaussian processes. Another class of methods assumes a specific parametric form of the underlying probabilistic density function that generates the data points, among which Gaussian Mixture Models are one of the most powerful approaches due to its simplicity and flexibility. In this section, we present the modeling of the CGPD using GMMs.

In this section, we discuss the significance of the work presented here, as well as some of the caveats and limitations of this analysis.

Traditional galaxy surveys and the relatively recent LIM techniques are both powerful methods to study the formation and evolution of galaxies, one of the most fundamental and important problems in astrophysics and cosmology. Motivated by the necessity of developing a physically motivated, computationally efficient framework that can describe both galaxy and LIM observations, we propose and characterize a new approach based on the CGPD $\Phi\left(\bm{p}|M,z\right)$ function, which describes the probability of a halo with mass $M$ at redshift $z$ to host a galaxy with properties $\bm{p}$. When combined with the Line Luminosity Relation (LLR) $L_{i}(\bm{p},z)$ function developed in a companion work (Breysse et al., 2023), a highly efficient and flexible pipeline that connects halos with an n-dimensional distribution of galaxy physical properties, and with LIM observables, can be constructed.

As a distribution function in the 5-D galaxy property space that we have chosen, the CGPD encapsulates many important galaxy distribution functions and scaling relations etc. Motivated by the form of the CGPD seen in two separate physics-based cosmological simulations (the Santa Cruz SAM and IllustrisTNG), we discuss a parameterized characterization of CGPD using GMM. We find the following main results:

• •

  With just a few ($\lesssim 3$) components, GMM is capable of reproducing the CGPD distribution from both the SAM and TNG to high accuracy.

• •

  The GMM based CGPD yields estimates of common summary statistics of interest including 2-D property relations (scaling relations), property distribution functions, and integrated cosmic densities, with a typical accuracy of a few percent.

• •

  For each bin in halo mass and redshift, we find that the 5-D property distribution can be reduced to only 3 dimensions using PCA, with a 10% loss of information.

The galaxy properties in this work were selected based on their importance for predicting the luminosity of certain emission lines, but a similar approach could be used for other galaxy properties, including galaxy photometry as measured in galaxy surveys.

The next step is to construct a statistical inference pipeline for LIM observations. For a 5-D Gaussian distribution with a full covariance matrix, we have 20 free parameters including 5 parameters in the mean vector and 15 independent parameters in the symmetric covariance matrix. Then if we use 3 Gaussian components for each halo mass bin and 15 mass bins for each redshift as in our main analysis, we would have 900 parameters for the GMM at each redshift. However, these parameters are expected to be highly degenerate. We have shown that the data dimensionality can be reduced significantly using PCA. Another possible approach for doing implicit likelihood inference in high dimensional parameter space is the marginal posterior density with moment networks as proposed in Jeffrey & Wandelt (2020).

We thank the IllustrisTNG team for granting the JupyterLab access to the simulations, and Dylan Nelson for help accessing the TNG data. ARP was supported by NASA under award numbers 80NSSC18K1014, NNH17ZDA001N, and 80NSSC22K0666, and by the NSF under award number 2108411. ARP and rss are supported by the Simons Foundation. The numerical computations in this work were performed in part using the NYU Greene High Performance Computing cluster.