The JWST and ALMA have detected emission lines from the ionized interstellar medium ISM including OII OIII and hydrogen Balmer series lines in some of the first galaxies at zgt6. These measurements present an opportunity to better understand galaxy assembly histories and may allow important tests of stateoftheart galaxy formation simulations. It is challenging however to model these lines in their proper cosmological context given the huge dynamic range in spatial scales involved. In order to meet this challenge we introduce a novel subgrid line emission modeling framework. The framework uses the highz zoomin simulation suite from the Feedback in Realistic Environments FIRE collaboration. The line emission signals from HII regions within each simulated FIRE galaxy are modeled using the semianalytic HIILines code. A machine learning Mixture Density Network approach is then used to determine the conditional probability distribution for the line luminosity to stellarmass ratio from the HII regions around each simulated stellar particle given its age metallicity and its galaxys total stellar mass. This conditional probability distribution can then be applied to predict the line luminosities around stellar particles in lower resolution yet larger volume cosmological simulations. As an example we apply this approach to the IllustrisTNG simulations at z6. The resulting predictions for the OII OIII and Balmer line luminosities as a function of starformation rate agree well with current observations. Our predictions differ however from related work in the literature which lack detailed subgrid ISM models. This highlights the importance of our multiscale simulation modeling framework. Finally we provide forecasts for future line luminosity function measurements from the JWST and quantify the cosmic variance in such surveys.
['Astrophysics - Astrophysics of Galaxies']
A New Framework for ISM Emission Line Models Connecting MultiScale Simulations Across Cosmological Volumes
{'A New Framework for ISM Emission Line Models: Connecting Multi-Scale Simulations Across Cosmological Volumes': 'Shengqi Yang, 1, 2 Adam Lidz, 3 Andrew Benson, 1 Yizhou Zhao, 4 Hui Li, 2 Amelia Zhao, 5 Aaron Smith, 6 Yucheng Zhang, 7 Rachel Somerville, 8 Anthony Pullen, 9 and Hui Li 10  \n1 Carnegie Observatories, 813 Santa Barbara Street, Pasadena, CA 91101, U.S.A  \n2 Los Alamos National Laboratory, Los Alamos, NM 87545, USA  \nDepartment of Physics and Astronomy, University of Pennsylvania, 209 South 33rd Street, Philadelphia, PA 19104, USA  \n4 UCLA Center for Vision, Cognition, Learning, and Autonomy  \n5 Pomona Colledge  \n6 Department of Physics, The University of Texas at Dallas, Richardson, Texas 75080, USA  \n7 Peng Cheng Laboratory, Shenzhen, Guangdong 518066, China  \n8 Center for Computational Astrophysics, Flatiron Institute, New York, NY 10010, USA  \n9 Center for Cosmology and Particle Physics, Department of Physics, New York University, 726 Broadway, New York, NY 10003, USA 10 Department of Astronomy, Tsinghua University, Beijing 100084, China', 'ABSTRACT': "The JWST and ALMA have detected emission lines from the ionized interstellar medium (ISM), including [OII], [OIII], and hydrogen Balmer series lines, in some of the first galaxies at z ≳ 6. These measurements present an opportunity to better understand galaxy assembly histories and may allow important tests of state-of-the-art galaxy formation simulations. It is challenging, however, to model these lines in their proper cosmological context given the huge dynamic range in spatial scales involved. In order to meet this challenge, we introduce a novel sub-grid line emission modeling framework. The framework uses the highz zoom-in simulation suite from the Feedback in Realistic Environments (FIRE) collaboration. The line emission signals from HII regions within each simulated FIRE galaxy are modeled using the semi-analytic HIILines code. A machine learning, Mixture Density Network (MDN), approach is then used to determine the conditional probability distribution for the lineluminosity to stellar-mass ratio from the HII regions around each simulated stellar particle given its age, metallicity, and its galaxy's total stellar mass. This conditional probability distribution can then be applied to predict the line luminosities around stellar particles in lower resolution, yet larger volume cosmological simulations. As an example, we apply this approach to the Illustris-TNG simulations at z = 6. The resulting predictions for the [OII], [OIII], and Balmer line luminosities as a function of starformation rate (SFR) agree well with current observations. Our predictions differ, however, from related work in the literature which lack detailed sub-grid ISM models. This highlights the importance of our multi-scale simulation modeling framework. Finally, we provide forecasts for future line luminosity function measurements from the JWST and quantify the cosmic variance in such surveys.  \nKeywords: galaxies: evolution - galaxies: high-redshift - submillimetre: ISM - (ISM:) H II regions", '1. INTRODUCTION': "The Atacama Large Millimeter/Submillimeter Array (ALMA) and the James Webb Space Telescopes (JWST) have detected multiple emission lines from some of the first galaxies, which formed less than a billion years af-  \nCorresponding author: Shengqi Yang  \nsyang@carnegiescience.edu  \nter the Big Bang (e.g. Inoue et al. (2016); Laporte et al. (2017); Hashimoto et al. (2018); Laporte et al. (2019); Hashimoto et al. (2019); Harikane et al. (2020); Witstok et al. (2022); Heintz et al. (2022); Curti et al. (2023); Sanders et al. (2023); Laseter et al. (2023)). Excitingly, the pace of discovery here is expected to accelerate over the next few years. In addition, SPHEREx will soon be launched and it will, among other measurements, carry out line-intensity mapping (LIM) surveys. This will uniquely complement the bright galaxy detec-  \n3  \ntions from the JWST and ALMA (Dor'e et al. 2014). Among the brightest nebular lines currently detected are collisionally-excited [OII] and [OIII] lines, and hydrogen Balmer series recombination lines. These lines are sensitive to the intensity of the local ionizing radiation and the properties of the surrounding interstellar gas, including its metallicity, density, and temperature. The ALMA, JWST, and LIM measurements can then provide new insights into galaxy formation at high redshift, and can, in-principle, provide strong tests of stateof-the-art galaxy formation models as implemented in numerical simulations.  \nIn order to best interpret the current and upcoming data, reliable emission line models across cosmological volumes are required. For example, one important goal is to use the observations to quantify correlations between the gas-phase metallicities and stellar masses of galaxies (the mass-metallicity relation, or 'MZR'). The amplitude, shape, scatter, and redshift evolution of this correlation are partly shaped by poorly understood feedback processes and the MZR hence provides an important empirical anchor for galaxy formation models (Tremonti et al. 2004; Lee et al. 2006; Maiolino et al. 2008; Mannucci et al. 2009; Zahid et al. 2013; Henry et al. 2021). However, to robustly extract this from current data, one requires models to infer gas-phase metallicities from the observable [OIII], [OII], and Balmer line luminosities, (and unbiased stellar mass estimates, although this is not our focus here). Ideally, this necessitates modeling variations in the properties of the emitting HII regions across each galaxy, as well as accounting for galaxy-to-galaxy fluctuations in the ISM properties. Furthermore, the model must span the observed range in stellar mass. Hence, to best extract the MZR from current data, among other goals, one needs to capture HII-region scale physics across cosmological volumes.  \nAlthough multiple strategies have been proposed for post-processing ISM emission lines on top of cosmological simulations, they are not ideal for making detailed comparisons between simulations and observations. The problem is two-fold. First, large volume cosmological simulations currently lack the resolution to capture the multi-phase ISM. Instead, a simplified two-phase description is used to characterize the sub-grid ISM, with 'hot' and 'cold' phases evolving according to effective equations of state (Springel & Hernquist 2003). This approach is inadequate for modeling the line emission from HII regions, which is partly shaped by the radiation field, metallicity, density, and temperature on sub-grid length scales. Second, it is expensive to perform line emission calculations across many millions of emitting regions using the Cloudy code (Ferland et al. 2017),  \nespecially if one properly accounts for variations in the stellar radiation field and in the local ISM properties. As a compromise, approximations are often adopted such as assuming a constant stellar radiation field and/or gas density across all emitting regions (e.g. Hirschmann et al. (2017); Moriwaki et al. (2018); Ceverino et al. (2021); Kannan et al. (2022); Hirschmann et al. (2023); Nakazato et al. (2023)). Although such simplifications help make the problem tractable, it is unclear how accurate they are. The simple approach also discards information from the simulations regarding the simulated star-formation histories, which should ideally be incorporated into the line emission modeling and lead to region-to-region variations in the ionizing spectral amplitude and shape. On the other hand, there are expensive small-volume cosmological simulations that do, at least partly, resolve the multi-phase ISM. These calculations, however, generally simulate volumes too small to capture rare, bright galaxies, yet these sources are often the most accessible with current observations. As several representative examples from current state-of-theart work in the literature, we note that Kannan et al. (2022) assumes a fixed HII region gas density across all galaxies and only accounts for stellar particles younger than 10 Myr. Likewise, Hirschmann et al. (2023) adopt a constant gas density and assume all galaxies in the TNG simulation suite have identical stellar spectra. Even the more ambitious approach developed in Katz et al. (2022) only models line emission signals for a sub-set of simulation cells using Cloudy and then relies on a machine learning technique to extrapolate predictions across all cells. Moreover, the SPHINX simulations employed in Katz et al. (2022) span relatively small volumes and hence lack galaxies with star formation rates (SFRs) that are high enough to compare directly with many current ALMA measurements.  \nIn this work we propose a new framework for postprocessing ionized ISM emission lines on top of cosmological simulations which may improve upon some of the above limitations. This framework combines the highz FIRE zoom-in galaxy simulations (Ma et al. 2018, 2019, 2020), the mixture density network (MDN) architecture (Bishop 1994), and an efficient ionized ISM emission line simulation methodHIILines (Yang et al. 2023a). The main scientific goal of this novel hybrid framework is summarized in Figure 1. We aim to use the MDN machine learning architecture to extract information about the small-scale structure of the ISM and the resulting emission line signals, conditioned on the properties of the local ionizing sources, using simulated FIRE highz galaxies that (partly) resolve individual HII regions. We then apply the the conditional probability distribution  \nfunctions (cPDFs) to the stellar particles in cosmological simulations to predict the unresolved line emission around them. In this manner we can robustly postprocess HII region emission line signals among galaxies in cosmological simulations, assuming a small-scale ISM environment identical to that in the FIRE zoom-in galaxies. That is, our model combines the populationlevel statistics of a cosmological simulation of galaxy formation with the small-scale ISM physics captured in FIRE. It is flexible enough for us to explore how varying some of the assumptions here impacts our final results. Ultimately, our framework will help fully exploit revolutionary new data sets from ALMA, JWST, and LIM surveys, across multiple wavebands, extracting their implications for early galaxy formation in its full cosmological context.  \nThe trained model is publicly available at https:// github.com/Sheng-Qi-Yang/GMDN, and is applicable to any numerical or semi-analytical galaxy formation simulation with spatial resolution no finer than that of FIRE. The plan of this paper is as follows. In Section 2 we introduce the FIRE zoom-in galaxies and HIILines , which are used to generate the training and validation datasets for this work. Section 3 introduces the MDN architecture and the model training method. Section 4 explains how the trained MDN can be applied to stellar particles in simulated galaxies. We verify our MDN model by showing that it successfully reproduces the 1D distributions, 2D correlations, and galaxy-wide total line emission signals for the 22 FIRE galaxies. We also test the assumption that the line luminosity given stellar particle property cPDFs from these 22 FIRE galaxies are representative of more general galaxy populations. In Section 5 we apply the trained MDN to well-resolved galaxies in the TNG50, TNG100, and TNG300 simulation suites at z = 6. We then compare our models with current observations, considering multiple different statistical measures, and contrast our results with other modeling approaches in the literature. We summarize the strengths, limitations, and future prospects of our emission line modeling framework in Section 6.", '2.1. FIRE highz galaxies post-processed by HIILines': 'The 22 central galaxies at z = 6 from the FIRE project, made publicly available by Ma et al. (2018, 2019, 2020), cover a wide range of stellar mass M ∗ ,gal ∼ 10 6 M ⊙ -10 10 M ⊙ , gas-phase metallicity Z ∼ 10 -3 Z ⊙ -0 . 6 Z ⊙ , and morphology. Since the FIRE zoom-in simulations can resolve the multi-phase ISM at ∼ 10 parsec spatial resolution, which is comparable to the size of individual HII regions, the ISM environment simulated  \nby FIRE is more reliable than those predicted by largevolume cosmological simulations. The FIRE zoom-in galaxy simulations can hence be combined with spectral synthesis calculations to model ISM emission line signals.  \nAs an example, Figure 2 compares the MZR at z = 6 among the FIRE highz suite (Ma et al. 2018, 2019, 2020), TNG50, TNG100, TNG300 (Pillepich et al. 2018a; Marinacci et al. 2018; Nelson et al. 2018; Springel et al. 2018; Naiman et al. 2018; Pillepich et al. 2018b; Nelson et al. 2019; Pillepich et al. 2019), and the Santa Cruz semi-analytic galaxy formation model (SAM) (Somerville & Primack 1999; Somerville et al. 2008, 2012; Porter et al. 2014; Popping et al. 2014; Somerville et al. 2015), together with direct metallicity measurements from JWST (Nakajima et al. 2023; Curti et al. 2023) and ALMA (Harikane et al. 2020; Yang & Lidz 2020) at z ≥ 6. When the [OIII] auroral line is undetected, metallicities are estimated through the O32+R23 strong line diagnostics using the method proposed by Yang et al. (2023a) (triangular points) (Nakajima et al. 2023). In each case, the metallicity of a simulated galaxy is defined as the SFR-weighted gas-phase metallicity of oxygen averaged over all gas particles, which is close to the metallicity constrained by the ISM emission line measurements (Yang et al. 2023b). Although the scatter among the observational data points is large, the gas-phase metallicities of the TNG galaxies lie above the data for M ∗ ≲ 10 8 . 5 M ⊙ , while the FIRE and SAM MZRs appear in better agreement with current measurements. Since the small-scale ISM properties in FIRE, or at least the gas-phase metallicities, are in reasonable agreement with observations across a broad range of stellar masses, we use FIRE to anchor our modeling framework on small scales.  \nThe next tool in our modeling is the HIILines code from Yang et al. (2023b) which yields efficient [OIII], [OII], H α , and H β line emission predictions, postprocessed on top of zoom-in galaxy formation simulations. HIILines assumes that every stellar particle in a simulated galaxy sources an independent HII region. The stellar radiation field is modeled using the Flexible Stellar Population Synthesis code (FSPS; Conroy et al. (2009); Conroy & Gunn (2010)), which takes the stellar particle age, metallicity, and mass as inputs. The local gas density and metallicity are then estimated through averaging over the 32 nearest gas particle neighbors around the stellar particle of interest. An empiricallymotivated model is used to describe the gas temperature since this quantity is not reliably modeled in the FIRE simulations. This model exploits the fact that the tem-  \nFigure 1. A schematic showing the main goal of this work. First, we use the HIILines code to model HII region emission line signals on top of FIRE zoom-in simulations of high redshift galaxies (top left). Next, an MDN approach is applied to the FIRE models to determine the joint cPDFs for the luminosities across several lines of interest, conditioned on the coarsegrained stellar particle properties (top right). This is then used to populate line emission signals around the stellar particles in lower resolution, yet larger volume, simulations (bottom right). That is, the cPDFs effectively provide a FIRE-informed sub-grid model for the fine-grained HII region and line emission properties, given the surrounding coarse-grained stellar particle properties. The resulting models can then be compared in detail with current and upcoming [OIII], [OII], and hydrogen Balmer line observations from the JWST, ALMA, and SPHEREx ( bottom left ).  \n<!-- image -->  \nrature of the HII region gas is well-correlated with its gas-phase metallicity.  \nSpecifically, we assume that the temperature of all of the OIII and OII gas in a given galaxy is constant although the OIII and OII regions generally differ from each other in temperature - and that these temperatures are determined by the Q HI -weighted gas-phase metallicity of the galaxy, Z Q . Here Q HI is the ionizing photon generation rate of the stellar particles. The OIII temperature versus metallicity relation is fit to highz direct T e measurements, while the OII temperature is fit to local measurements (Yates et al. 2020). HIILines then applies a semi-analytic line emission model, first proposed by Yang & Lidz (2020), to compute the [OIII], [OII], H α , and H β line luminosities of each HII region. The semi-analytic model is ∼ × 1000 times faster than numerical spectral synthesis codes such as Cloudy , and is therefore suitable for post-processing millions of line emitters with different ionizing radiation fields and gas properties. Thanks to its high computational efficiency, HIILines is currently the only model that accounts for  \nspatial variations in the local gas density, metallicity, and incident radiation spectrum, along with their scatter across the HII regions within simulated galaxies. The simulated [OIII], [OII], H β line luminosities and line ratios among the 22 FIRE highz galaxies are in excellent agreement with recent ALMA and JWST measurements at z ≳ 6.  \nIn this work we suppose that the ISM density and metallicity distributions within these 22 FIRE galaxies are representative and apply to any galaxy with similar mass/metallicity, regardless of their morphology. In other words, we assume that the ISM line luminosity cPDFs extracted from these FIRE zoom-in galaxies are universally applicable, and can be applied to arbitrary members of the simulated galaxy populations across large-volume cosmological simulations. This is a strong assumption, yet it lies at the foundation of our method for constructing sub-grid ISM line emission models. We will quantitatively test this assumption in Section 4.1.  \nFigure 2. MZR at z = 6 from FIRE, TNG, and the Santa Cruz SAM, as compared to direct T e measurements. The metallicity of each simulated galaxy is defined as the SFRweighted gas-phase metallicity. The colored bands give the 1 σ dispersion around the MZRs for the TNG and SAM models. The FIRE highz suite contains only 22 central galaxies, and we therefore only show the best fit MZR in this case. Observational measurements are shown by the grey points and error bars. The TNG MZR lies above the observations at M ∗ ≲ 10 8 . 5 M ⊙ .  \n<!-- image -->  \nEach stellar particle in FIRE has a known radiation field, determined through stellar population synthesis calculations, and is assumed to source an independent HII region. We then model the luminosities for 13 of the [OIII], [OII], H α , and H β lines emitted by each HII region. Our goal is to determine the following joint cPDF:  \nϕ = p (log( L [OIII] , [OII] , H α, H β /m ∗ ) | log t, log Z ∗ , log M ∗ ) . (1)  \nHere L [OIII] , [OII] , H α, H β describes the luminosities of the 13 [OIII], [OII], and hydrogen recombination lines emitted by an HII region in a FIRE galaxy, while m ∗ denotes the mass of the stellar particle that sources the particular HII region. We do not fit for the line luminosity distribution directly, but instead fit for the ratio of line luminosity to stellar mass because we want the trained model to be applicable to simulations of arbitrary resolution. The quantities t and Z ∗ are the age and metallicity of the stellar particle, respectively. Given a stellar population synthesis model, these two input parameters fully determine the shape of the local ionizing radiation field. Following Yang et al. (2023b), we model stellar radiation using the FSPS code. Finally, M ∗ is the total stellar mass of a FIRE galaxy. We introduce M ∗ as a third input parameter because it provides information about the gas temperature. As introduced in Section 2, HIILines estimates the galaxy-wide HII region temperature based on Z Q . Since Z Q is somewhat expensive to calculate directly, we instead use the MZR as a short-cut  \nhere: given the correlation between Z Q and M ∗ in FIRE - which agrees well with observations (Figure 2) M ∗ implies information about Z Q , which in turn specifies the model gas temperatures.  \nThe 22 FIRE galaxies contain 3 . 3 × 10 6 HII regions in total, sourced by stellar particles younger than 100 Myr. We randomly sample 80% of the HII regions to form the training dataset. The remaining 20% of the HII regions will be used for model validation. We do not account for HII regions sourced by extremely old stellar populations in this work as they are faint and contribute negligibly to the total galaxy-wide line luminosities (e.g. Yang et al. (2023b)).', '3. GALAXY-MDN (GMDN)': 'In this section we combine the FIRE data postprocessed by HIILines with a machine learning architecture named MDN to extract the connection between stellar particle properties and the corresponding HII region line signals. The combined pipeline here will be referred as Galaxy-MDN (GMDN).  \nMDN is a machine learning architecture that fits the cPDF of the training data with a Gaussian mixture model (GMM, i.e., by a sum of multiple multi-variate Gaussian distributions), conditioned on the input variables (Bishop 1994). Through randomly sampling over the GMM, MDN serves as a generative model and provides outputs with cPDFs approximately matching the distribution of the training data.  \nWe summarize the overall flow of the GMDN in Figure 3. As introduced in Eq 1, each stellar particle is assumed to be surrounded by an HII region and carries a 3-dimensional input array: x = { log t, log Z ∗ , log M ∗ } . This input array, after standardization, will first go through a network layer for linear transformation:  \nX = W in x + b in . (2)  \nHere W in and b in are matrices of size 512 × 3 and 512 × 1 respectively. The size of the network is a hyperparameter. We select 512 due to its balance between computational efficiency and memory usage, which is well-suited for our GPU computing framework. We then add non-linearity to the network result:  \nz h = tanh( X ) . (3)  \nIntroducing non-linearity between each layer of linear transformation is necessary. Without this step, the combination of multiple layers of linear transformation is effectively only a single layer transformation, and the free parameters introduced by each network layer become redundant.  \nFigure 3. GMDN work flow. Given the three input parameters for each stellar particle: age, t ; metallicity, Z ∗ ; and galaxy total stellar mass, M ∗ , GMDN constructs a Gaussian Mixture Model (GMM) with weight, Π, standard deviation, σ , and mean, µ , to describe the output, y . Here y is the standardized line luminosity versus stellar particle mass ratio log( L/m ∗ ) after PCA transformation. GMDN then randomly samples from the GMM to generate line signals for any given stellar particle. The advantage of this framework is that the local ISM environment from FIRE is encoded in the free parameters of the linear networks during the model training. The GMDN can then generate line luminosity models given only our three input parameters (stellar age, metallicity, and the total stellar mass of a galaxy), without requiring further information regarding the surrounding ISM.  \n<!-- image -->  \nSeveral linear networks are then used to compute the weight, mean, and standard deviation for each Gaussian distribution that makes up the GMM. We will use N Gauss separate Gaussian distributions to describe the cPDF of the training data set. Ideally, each Gaussian should be 13-dimensional, corresponding to the 13 [OIII], [OII], and hydrogen recombination lines covered by HIILines . However, the luminosity distributions of these 13 lines are not independent of each other. Instead, they are highly positively correlated. In this case, the GMM requires a covariance matrix to capture the correlations among different lines. In order to simplify the GMM calculations, we perform a principle component analysis (PCA): we determine the eigenvectors of the covariance matrix and rank them according to their eigenvalues. We find that the first 4 PCA components account for more than 99% of the total data variance and that the covariance matrix can therefore be well approximated by including only the first N PCA = 4 eigenmodes. In the eigenmode basis, in which the covariance matrix is diagonal, the GMM may be described as:  \nϕ ( y ) = N Gauss ∑ k =1 Π k (2 π ) N PCA / 2 ∏ N PCA i =1 σ i k ( x ) × exp ( -1 2 N PCA ∑ i =1 ( y i -µ i k ( x )) 2 ( σ i k ( x )) 2 ) . (4)  \nHere y i describes the i th PCA eigenmode: these are suitable linear combinations of the line luminosities, log( L/m ∗ ).  \nThe weight of each Gaussian component Π is computed through:  \nZ Π = W Π z h + b Π , Π = exp( Z Π ) Σexp( Z Π ) . (5)  \nHere W Π and b Π are matrices of size N Gauss × 512 and N Gauss × 1, respectively. The mean and standard deviation of each Gaussian for each PCA component µ and σ are computed through:  \nZ µ = W µ z h + b µ , µ = Z µ , (6)  \nand  \nZ σ = W σ z h + b σ , σ = exp( Z σ ) . (7)  \nHere W µ and W σ are matrices of size N Gauss × 512. The quantities b µ and b σ are matrices of size N Gauss × 1. There are, in total, 8 of these linear networks since each Gaussian is 4-dimensional. The full GMDN contains in total  \ndim( θ ) = 2048 + 4617 N Gauss (8)  \nfree parameters. The number of Gaussian distributions used to build the GMM, N Gauss , is another hyperparameter.', '3.1. Loss function': 'Our goal is to train the GMDN so that, given the age and metallicity of a stellar particle from any simulation, together with the galaxy-wide stellar mass, the GMDN can generate random line luminosities with distributions and correlations identical to the FIRE simulation results. First, the PCA transfers all standardized FIRE HII region line luminosity-to-stellar mass ratios from 13D to 4D. We then define the first part of the loss function as the following:  \nL 1 = ⟨ ω ( t j ) ln ϕ ( y j FIRE ) ⟩ N HII j =1 , (9)  \nfor both the training and validation datasets. Here ⟨⟩ N HII j =1 denotes the mean among all HII regions in the  \ntraining or validation dataset, while ϕ ( y ) is defined in Eq 4. The quantity ω ( t j ) is a normalized weight for an HII region, which is determined by the stellar particle age. Specifically, we define the weight of an HII region to be inversely proportional to the stellar particle age PDF. We make this choice because young stellar populations are rare, yet they source the brightest HII regions, which are major contributors to the galaxy-wide total line luminosities. On the other hand, old stellar populations dominate the galaxy-wide stellar mass, yet they source faint HII regions and make negligible contributions to the total line luminosities. By introducing ω ( t ) in the loss function, we assign higher weights to the young stellar particles to better capture the total, galaxy-wide luminosities.  \nTo ensure that the strong line diagnostics generated by our GMDN are close to the FIRE simulation results, we draw random samples from ϕ and do an inverse PCA transfer to get 13D outputs. We then define the second part of the loss function as:  \nL 2 = 5 ∑ l =1 〈 ω ( t j ) ×| ln R l MDN -ln R l FIRE | 〉 N HII j =1 , (10)  \nwhere R 1 to R 5 denote the following strong line diagnostics: L [OIII] 5007 /L [OIII] 4364 , L [OII] 3729 /L [OII] 3727 , L [OIII] 88 µm /L [OIII] 52 µm , L [OIII] 5007 /L H β , L H α /L H β . The total loss function is:  \nL = L 1 + L 2 . (11)  \nWe introduce line ratio constraints into the loss function because these five line ratios are sensitive diagnostics of ISM properties, including the gas temperatures in the OII and OIII regions, the gas densities, metallicities, and of dust attenuation effects.  \nWe adopt the Adam optimization algorithm with an initial learning rate of 10 -5 to train the GMDN (Kingma & Ba 2014). During the training, we feed the GMDN with FIRE data of batch size 10 5 in each iteration. Each training is stopped after 30,000 epochs so that the loss function change after every 100 epochs of training is less than 0.01 as the training comes to an end. To determine the best value of N Gauss we compute the Bayesian Information Criterion (BIC):  \nBIC = dim( θ ) ln n +2 n L . (12)  \nHere dim( θ ) is the number of free parameters defined in Eq 8, while n is the number of HII regions in the model validation dataset. L is the loss function of the validation dataset for the best fit parameter θ , as defined in Eq 11. Smaller values of BIC generally correspond to more preferred models. We present the BICs under different N Gauss in Figure 4. The BIC test suggests that the optimal choice is N Gauss = 10.  \nFigure 4. The BIC of our GMDN for differing numbers of Gaussian distributions in the mixture model. This suggests that N Gauss = 10 is the optimal choice.  \n<!-- image -->', '4. MODEL VERIFICATION': 'We use the trained GMDN as an [OIII], [OII], and hydrogen recombination line emulator in the following way. Given the age and metallicity of a stellar particle, and the stellar mass of the simulated galaxy, our GMDN gives the PDF for the 4 PCA components ϕ ( y ). Since the GMM is formed by 10 4D Gaussian distributions, we first Gumbel sample to pick out one Gaussian. If the k th Gaussian distribution is sampled, we pick a random number r from a normal distribution N (0 , 1). The 4 PCA components are then computed as:  \nlog( L i PCA /m ∗ ) = σ i k r + µ i k , i = 1 , 2 , 3 , 4 . (13)  \nWe then perform an inverse PCA transform to this 4D random sample to derive the 13D values of log( L/m ∗ ). Finally, we multiple this 13D array by the stellar particle mass to derive the luminosities of the 13 modeled [OIII], [OII], H α , and H β lines emitted by the HII region sourced by this stellar particle. To avoid the rare case where the random sampling process hits the shallow tail of the GMM and gives unrealistically high luminosity predictions, we use the analytical model introduced in Yang et al. (2023a) to estimate maximum values of the [OII] 3729 ˚ A and [OIII] 5007 ˚ A lines:  \nL [OII] 10 , max = ( n O n H ) ⊙ Z Z ⊙ k OII 10 hν OII 10 Q HI α B , HII , L [OIII] 32 , max = ( n O n H ) ⊙ Z Z ⊙ A OIII 32 k OIII 32 hν OIII 32 A OIII 31 + A OIII 32 Q HI α B , HII . (14)  \nHere to compute the upper bounds on line luminosity we assume that the gas density is much lower than the critical densities of all lines. We further assume that the OII or OIII regions span the entire HII region. In Eq 14, ( n O /n H ) ⊙ = 10 -3 . 31 is the oxygen abundance at  \nSolar metallicity (Holweger 2001). The quantities A X ij , k X ij , and µ X ij are the spontaneous decay rate, collisional de-excitation rate, and emitted photon frequency, respectively, when particle X transitions from level i to level j . We assume case-B recombination rates throughout, and α B , HII denotes the case-B HII recombination rate, while α B , H β /α B , HII gives the fraction of HII recombination events that lead to H β line emission. All atomic data are temperature dependent. When estimating the line luminosity upper bounds we first estimate the galaxy Z Q through the FIRE MZR at z = 6:  \nlog( Z Q / [ Z ⊙ ]) = 0 . 39 log( M ∗ / [ M ⊙ ]) -3 . 97 . (15)  \nWe then adopt metallicity-dependent OIII and OII temperature models from Yang et al. (2023a); Yates et al. (2020). If the random sample shows a line luminosity higher than the theoretical maximum value above, we abandon this result and redo the sampling until the HII region line luminosities are all below the upper bounds.  \nIn this section we verify the performance of the best fit GMDN. First, Figure 5 compares the luminosity distributions for the [OII] 3729 ˚ A, [OIII] 5007 ˚ A, and H β lines among all of the HII regions within the largest FIRE highz galaxy, z5m12b. The red contours show the luminosity distributions given by the FIRE simulation, while the blue contours show randomly generated data distributions from the GMDN. Figure 5 demonstrates that the GMDN successfully reproduces the 1D distribution and 2D correlations among [OIII], [OII], and hydrogen recombination lines.  \nFigure 6 further compares the galaxy-wide total (i.e. summed over all of the HII regions in a simulated galaxy) [OII] 3729,29 ˚ A, [OIII] 88 µ m, 4960 ˚ A, 5007 ˚ A, 4364 ˚ A, and H β line luminosities given by the FIRE simulations ( x axis) and the GMDN ( y axis). These lines have been detected by ALMA and JWST at z > 6, with additional observational data forthcoming in the near future. The black dashed line marks equality. The 22 FIRE galaxies used in training the GMDN are presented as black points. Figure 6 shows that our GMDN accurately reproduces the galaxy-wide total line luminosities. The H β line predictions are the most accurate: the maximum log L difference between the GMDN and FIRE simulation results is only 0.04 dex. Furthermore, the GMDN performs well for lines that are weakly sensitive to gas temperature and ionization correction factor effects 1 , including the [OIII] 88 µ m, 4960 ˚ A, and 5007 ˚ A lines. The mean of the absolute differences of log L at  \nFigure 5. Luminosity distributions for the simulated and GMDN-predicted emission in the [OII] 3729 ˚ A, [OIII] 5007 ˚ A, and H β lines for the FIRE galaxy z5m12b. The red contours show distributions from FIRE+ HIILines , while the blue contours show the distributions of random samples generated by our GMDN. Our GMDN can accurately reproduce the distributions and correlations among [OIII], [OII], and hydrogen recombination lines.  \n<!-- image -->  \nlog( L/L ⊙ ) > 4 is 0.1, 0.07, and 0.07 dex for these three lines, respectively. The GMDN predictions for the [OII] 3727,29 ˚ A and the [OIII] 4364 ˚ A auroral line are the least accurate. This is because we have used M ∗ as an input to indirectly trace the gas temperature. However, M ∗ and Z Q are not perfectly correlated. As a result, our GMDN method imperfectly estimates the gas temperature of simulated HII regions and is less successful when applied to temperature-sensitive lines. Overall, however, the trained GMDN accurately emulates the simulated line luminosities: it may therefore help with interpreting highz ALMA and JWST observations, and in modeling LIM signals.', '4.1. Test of Assumptions': "Our aim is to apply the trained GMDN to cosmological simulations. This involves an implicit assumption that the 22 FIRE galaxies we use to model HII region line emission, and their underlying distributions of HII region gas density and metallicity, are representative. We intend to apply the FIRE-trained distributions as a sub-grid model on top of entire galaxy populations in larger volume cosmological simulations. The potential concern here is that the ISM and line emission properties  \nFigure 6. Galaxy-wide total [OII], [OIII], and H β line luminosities for the 22 FIRE highz galaxies as given by HIILines ( x axis) and GMDN ( y axis). The black dashed lines mark equality. Each black point represents one simulated galaxy. Our GMDN successfully reproduces the galaxy-wide [OIII] 88 µ m, 4960 ˚ A, 5007 ˚ A, and H β line luminosities. It is less accurate for the fainter and more temperature sensitive [OII] 3727,29 ˚ A and [OIII] 4364 ˚ A lines.  \n<!-- image -->  \nmay vary more broadly than captured in our relatively small FIRE training sample.  \nAs one quantitative test of this assumption, we can remove galaxies from the training sample and test whether the line emission in these galaxies is well-reproduced by the sources in the new, smaller training set. The galaxies we removed from the training sample are marked with red crosses in Figure 7. These particular sources are excised because their stellar masses, M ∗ , are similar to those of the other FIRE galaxies which are retained in the training set. If the remaining 13 FIRE galaxies share similar ISM properties, then the GMDN should still accurately reproduce the line luminosities of the 9 sources which were excised from the training data. Indeed, Figure 7 illustrates that the reduced training data still reproduces the galaxy-wide line luminosities for the 9 excised galaxies (red crosses). As before, the performance is particularly good for the lines that are less sensitive to gas temperature and ionization correction factor effects, including the [OIII] 88 µ m, 4960 ˚ A, 5007 ˚ A, H α , and H β lines. The average absolute difference in the logarithm of the luminosity between when the smaller set of galaxies and the full set are used in training is: log( L/ [ L ⊙ ]) > 4 are 0.1, 0.1, 0.1, 0.03, and 0.03 dex for these lines, respectively.  \nThis test supports the notion that the conditional line luminosity PDFs involved in our analysis are relatively universal: the FIRE sample includes galaxies with a variety of morphologies, formation histories, and environmental properties, yet their line luminosities may be ro-  \nbustly captured when only a sub-set of the full sample are used in the training data. As an illustration, Figure 8 shows visualizations of the FIRE galaxies z5m11h (included in our model training set) and z5m11g (excluded from the training data). In the visualizations, each galaxy is viewed along an axis that lies in the direction of the galaxy's total angular momentum vector. In each subplot, the blue (red) contours enclose 50% and 64% of the total galaxy stellar mass (H β line luminosity) for each galaxy. In each panel, the grey-scale gives the stellar surface-mass density distribution. For reference, Table 1 gives the stellar mass, H β line luminosity, and half-mass radius for these two example galaxies. Although z5m11h and z5m11g have similar stellar masses, their morphologies and HII region line signal distributions differ markedly. Specifically, z5m11h is relatively compact and contains only one component, while z5m11g contains two main clumps, leading to a much larger half-mass radius. Furthermore, the galaxy-wide H β line luminosities in these two simulated galaxies differ by a factor of 2.6. Nevertheless, our GMDN model - which include z5m11 in its training set - successfully reproduces the HII region line luminosities for z5m11g even when it is excised from the training data. This supports our assumption regarding the universality of the cPDFs employed, even considering galaxies across a broad range of size-scales and morphologies.", '5. APPLYING THE GMDN TO TNG GALAXIES': 'As one example application, we apply the trained GMDNtogalaxies in the publicly available IllustrisTNG  \nFigure 7. Same as Figure 6, but here the GMDN is trained only with the 13 FIRE galaxies marked by the black points. The GMDN trained over this smaller dataset successfully reproduces the galaxy-wide [OIII], [OII], and hydrogen recombination line luminosities for the 9 galaxies that are not used in the model training (red crosses). The model performance is particularly good for the [OIII] 88 µ m, 4960 ˚ A, 5007 ˚ A, and hydrogen recombination lines.  \n<!-- image -->  \nFigure 8. Face-on visualizations of the stellar mass surface density distributions for FIRE galaxies z5m11h and z5m11g. The blue and red contours enclose 50% and 64% of the total stellar mass and H β line flux for each galaxy. The example galaxies, z5m11h and z5m11g, have similar stellar masses but differ strongly in their morphologies. The fact that the GMDN trained on one of these galaxies also reproduces the HII region line luminosities for the other galaxy supports our approach. The conditional line luminosity PDFs given stellar particle age, metallicity, and mass appear relatively universal, and independent of large-scale galaxy morphology and environmental properties. Note the difference in the range of scales on the x and y-axes in the examples of the left/right-hand panels.  \n<!-- image -->  \nTable 1. Stellar mass, H β line luminosity, and half-mass radius for FIRE galaxies z5m11h and z5m11g at z = 6.  \nsimulations. IllustrisTNG is a suite of state-of-theart cosmological galaxy formation simulations (Pillepich et al. 2018a; Marinacci et al. 2018; Nelson et al. 2018; Springel et al. 2018; Naiman et al. 2018; Pillepich et al. 2018b; Nelson et al. 2019; Pillepich et al. 2019). The publicly available IllustrisTNG simulations include cubic volumes with co-moving box lengths of 50, 100, and 300 cMpc, referred to as TNG50, TNG100, and TNG300 (with the smaller volume simulations having higher resolution). The coupled dynamics of dark matter and gas are tracked using the quasi-Lagrangian code AREPO. The simulations include models for gas heating and cooling, star formation, stellar feedback, and the formation and feedback from supermassive black holes, among other ingredients. The range of box size and resolution spanned by TNG50, TNG100, and TNG300 allows one to study both rare, bright galaxies and the numerous, yet smaller and less-luminous members of the galaxy population. These simulations hence provide powerful models for studying the statistical properties of galaxy populations, and this is potentially useful for comparing with JWST and ALMA emission line measurements, and near future line-intensity mapping survey data. As we have emphasized previously, however, TNG does not resolve the multi-phase structure of the ISM, let alone individual HII regions. Fortunately, our trained GMDN calculations provide FIREcalibrated sub-grid models which can be employed to post-process realistic line emission signals on top of the stellar particles in the TNG simulations.  \nWe select galaxies at z = 6 in TNG50, TNG100, and TNG300 that are resolved with at least 100 stellar particles. We also require that the TNG galaxies selected for post-processing contain at least one stellar particle younger than 5 Myr. This is to ensure that recent star formation activity takes place in the selected galaxies. This criterion leaves us with 7057 galaxies in TNG50, 3397 galaxies in TNG100, and 4364 galaxies in TNG300. We note that the 22 FIRE galaxies span a stellar mass range from 10 6 M ⊙ to 10 10 M ⊙ , but 0.04%, 0.7%, and 3% of TNG galaxies show stellar mass higher than this range in TNG50, TNG100, and TNG300, respectively. Since we use galaxy stellar mass as an input parameter to (indirectly) capture the HII region gas temperatures, and observational results show mild temperature variations at high metallicity/mass (Yang et al. 2023a), we set a stellar mass upper limit of 10 10 M ⊙ for the GMDN when it estimates L/M ∗ to avoid the need for any extrapolations.', '5.1. Individual galaxy comparisons': 'Figure 9 compares the post-processed ISM emission line luminosities from each TNG galaxy with current reionization-era observations (Harikane et al. 2020; Witstok et al. 2022; Heintz et al. 2022; Curti et al. 2022; Sanders et al. 2023; Laseter et al. 2023). Dust attenuation effects have been corrected for the JWST galaxies ID4590 and ID10612 from Curti et al. (2022), and all sources reported by Sanders et al. (2023); Laseter et al. (2023), assuming an extinction curve with R V = 2 . 5 (Draine 2011). The other JWST galaxies in the figure are measured to have negligibly small dust extinctions with A V < 0 . 25 mag (and most such sources have less than 0.1 magnitude of extinction). The observational results (grey crosses) generally agree well with the TNG results, as shown in the red, yellow, and blue bands which indicate the line luminosities post-processed on top of the TNG50, TNG100, and TNG300 simulations, respectively. The bands encompass the 16% to 84% range in simulated line luminosities, while the solid lines indicate the median line luminosities. However, the bottom right panel of Figure 9 shows that our method underestimates the [OIII] 4364 ˚ A line luminosities by ∼ 1 σ . This is because, at a fixed SFR, the 22 FIRE galaxies tend to show higher metallicities and lower HII region temperatures than the measured galaxies, and this leads to underestimates of the temperature-sensitive [OIII] 4364 ˚ A line luminosities. In summary, our trained GMDN, which uses FIRE to develop an ISM line emission sub-grid model for TNG, matches current observations fairly well. The main shortcoming is in simulating  \nhighly gas temperature-sensitive lines, where our current method appears less accurate.', '5.2. L (SFR) relations': 'In general, the line luminosities from our simulated galaxies are well-correlated with their SFRs: galaxies with high SFRs produce copious numbers of ionizing photons (via their young and hot O and B stars), which source HII regions and the lines of interest for our study. Of course the trend of line luminosity with SFR and the level of scatter around the mean relation vary from line to line. Here we explore these trends using our simulated TNG line-emitting galaxies, and compare with earlier work in the literature. Ultimately, one related goal here is to use LIM observations to extract the luminosity density as a function of redshift in some of these lines (e.g. Gong et al. (2017); Sun et al. (2021); Bernal & Kovetz (2022)). Provided the L (SFR) relation is understood or empirically-calibrated, these observations can then be used to determine the star formation rate density as a function of redshift, including contributions from low luminosity sources which lie beneath the flux limits of traditional surveys.  \nFigure 10 compares the line luminosity versus SFR relation predicted by TNG+GMDN and Kannan et al. (2022) (hereafter Kannan2022). The blue bands shows the median and 1 σ scatter in the intrinsic line luminosities predicted by our method, while the grey bands show post-processed line signals after accounting for dust attenuation at the wavelengths of interest for each emission line. In order to model dust attenuation we use the observationally-motivated dust Model A from Vogelsberger et al. (2020), which takes the absolute UV magnitude of each simulated galaxy as input. We assume an extinction curve from Cardelli et al. (1989) with an average Milky Way extinction per reddening of R V = 3 . 1. On the other hand, Kannan2022 applies the more sophisticated and expensive Model C of Vogelsberger et al. (2020), where a Monte Carlo radiative transfer code is used to model the dust attenuation.  \nHere we want to re-emphasize the main difference between this work and most other methods that postprocess ISM emission lines on top of cosmological simulations. In most earlier work, including Kannan2022, the authors either directly adopt the often poorly resolved ISM properties in the simulations employed, or adopt simplifications such as assuming constant density HII regions and/or incident radiation fields. These simplifications were necessary in part because of the limited resolution of the simulations and partly because the line emission calculations are generally performed using Cloudy or similar codes. The Cloudy calcu-  \nFigure 9. Comparison of intrinsic galaxy-wide HII region line luminosities between TNG galaxies (color bands) and observations (grey points). The [OIII], [OII] rest-frame optical line signals and H β line luminosities are from recent JWST measurements Heintz et al. (2022); Curti et al. (2022); Sanders et al. (2023); Laseter et al. (2023). The [OIII] 88 µ mline luminosities and galaxy SFRs are from ALMA data (Harikane et al. 2020; Witstok et al. 2022). The TNG galaxies post-processed with our GMDN are generally in agreement with observations, although the GMDN under-predicts the [OIII] 4364 ˚ A auroral line luminosities by ∼ 1 σ .  \n<!-- image -->  \nlations become extremely expensive when each stellar or gas particle in the simulation must be modeled separately, and so this forces simplifications such as assuming that all emitting HII regions have identical gas densities. Our approach refines the earlier line signal modeling efforts through introducing a sub-grid model to capture the small-scale ISM properties which are not directly resolved in the large volume simulations. Specifically, our GDMN allows us to produce random samples of the line emission around stellar particles in a larger volume cosmological simulation, with the statistical properties of this line emission matching those of the 22 FIRE zoomin galaxies. The line emission from the FIRE galaxies is in good agreement with current z ≳ 6 line luminos-  \ny observations. Additionally, the computationally efficient HIILines allows the FIRE galaxy line emission calculations to fully account for variations in the stellar radiation spectral shape, amplitude, ISM gas density, and metallicity. As discussed previously, this is in contrast to previous Cloudy -based analyses which generally require simplifying assumptions (such as assuming a constant gas density and radiation spectrum across all HII regions), owing to their computational expense.  \nWe fit the mean line luminosity (while including a suppression factor to account for dust attenuation at the relevant wavelengths for each line) versus SFR relations given by the GMDN with a double power law (Padmanabhan 2018; Yang et al. 2022):  \nL [ L ⊙ ] = 2 N SFR [ M ⊙ / yr] [ ( SFR / [M ⊙ / yr] SFR 1 ) -α + ( SFR / [ M ⊙ / yr] SFR 1 ) β ] -1 . (16)  \nThe best fit parameters are summarized in Table 2. There are several interesting differences between the  \nbest fit relations from this work (black solid curves) and  \nFigure 10. Line luminosity versus SFR relations given by TNG50, TNG100, and TNG300 in our models, compared to simulation results from Kannan et al. (2022) (red dash dotted curves). The blue bands show the median and 1 σ scatter of our intrinsic (i.e. without dust extinction) GMDN L (SFR) relations, while the grey bands include dust attenuation estimates at the relevant wavelengths. The best fit double power-law relations from our dust-attenuated GMDN L (SFR) models are shown as solid black curves. The lime dashed curves show L (SFR) models based on applying HIILines directly on top of TNG, without including the FIRE-calibrated sub-grid models.  \n<!-- image -->  \nKannan2022 (red dash-dotted curves) in Figure 10. In Kannan2022, all of the L (SFR) relations for the optical lines show shallower slopes at high SFR. This is caused by the stronger dust attenuation effect in metal enriched galaxies. The L (SFR) slopes at high SFR are generally different between these two works. This is likely caused by our different dust model assumptions. In our work, L [OII] (SFR) also bends downward at high SFR partly because the high SFR galaxies tend to contain younger stellar populations which emit harder ionizing radiation spectra. This, in turn, increases the fractional OIII volume while it decreases the OII volume and [OII] luminosity within the HII regions in these galaxies. Furthermore, the high SFR galaxies also generally have higher gas densities, which also leads to larger OIII fractional volumes (Yang et al. 2023b). In the case of the sub-mm [OIII] lines, dust attenuation is negligible, and so Kannan2022 adopt linear fits for L [OIII]88 (SFR) and L [OIII]52 (SFR). However, in our models L [OIII]88 (SFR) and L [OIII]52 (SFR) still show small suppressions (relative to linear relations) at high SFR. In our calculations this occurs again because high SFR galaxies tend to be denser. The higher collisional de-excitation rates at high density lead to transitions without photon emission and a reduction in the line luminosity to SFR ratio (Yang & Lidz 2020). This effect is missed in Kannan2022 be-  \nTable 2. Best fit parameters for the dust attenuated L (SFR) relations parameterized by Eq 16.  \ncause the ISM gas density is fixed in that work to n e =100 cm -3 , and this is lower than the critical densities of these lines. However, in our FIRE-calibrated subgrid model there are numerous HII regions at higher densities, especially within high SFR galaxies. Finally, note that our calculations assume sub-grid metallicities from FIRE, while while Kannan2022 adopt metallicites directly from the TNG simulations. This accounts for some of the differences between our models, as the FIRE and TNG MZRs are disparate (see Figure 2).  \nWe can obtain a more quantitative understanding of the differences between the ISM properties in FIRE and TNG, and how these influence the L (SFR) relations.  \nFirst, we can directly adopt the TNG ISM properties and apply HIILines to predict L (SFR) without including our FIRE-calibrated sub-grid model. These calculations should resemble previous work in the literature, including Kannan2022. In order to explore this, we must first determine some relationships between the rate of ionizing photon production and SFR among TNG galaxies, and their gas phase metallicities. First, the top panel of Figure 11 shows the median TNG galaxy Q HI versus SFR relation, and its 1 σ scatter. For each TNG galaxy, we derive Q HI by summing over all stellar particles younger than 100 Myr in age, after assuming an FSPS stellar population synthesis model and a Chabrier initial mass function (IMF). The resulting best-fit mean Q HI (SFR) relation (black curve) is given by:  \nlog ( Q HI [ s -1 ] ) = 0 . 932 log ( SFR [ M ⊙ / yr] ) +53 . 4 . (17)  \nFor reference, we note that the average TNG Q HI (SFR) relationship is close to the fitting formula from Schaerer (2003) for a constant SFR model with a Salpeter IMF and stellar population ages larger than 6 Myr. (See the magenta band, where the upper and lower bounds adopt constant stellar metallicities between 0.1 and 1 Z ⊙ , respectively.) Further, the average SFR-weighted gas-phase metallicity, which serves to mimic the Q HI -weighted metallicity, for all TNG galaxies is shown in the bottom panel of Figure 11. The best fit result is:  \nlog ( Z SFR [ Z ⊙ ] ) = 0 . 0411 log ( SFR [ M ⊙ / yr] ) 2 +0 . 202 log ( SFR [ M ⊙ / yr] ) -0 . 527 . (18)  \nWe can then estimate the HII and OIII region gas temperatures with the empirically-motivated relation (Yang et al. 2023a):  \nT 4 = 0 . 824(log Z SFR ) 2 +0 . 101 log Z SFR +1 . 08 . (19)  \nHere Z SFR is in units of Z ⊙ , and T 4 is the ionized ISM gas temperature in units of 10 4 K.  \nUsing these relationships as input, and further assuming a constant ISM gas density of n e = 100 cm -3 and that the volume of each OIII region matches the HII region volume, V OIII /V HII = 1, we compute the intrinsic L (SFR) relations for all of the [OIII] and hydrogen recombination lines from HIILines . The TNG+HIILines results, which ignore our FIRE sub-grid modeling steps, are shown by the lime dashed curves in Figure 10. First, note that the GMDN predictions for the H α and H β luminosities overlap strongly with the simplified TNG+HIILines models. This occurs because the hydrogen recombination lines are fairly insensitive to local  \nISM properties, such as the HII region gas density and metallicity. Instead, the luminosities in the Balmer lines are mostly determined by the Q HI of each galaxy, and this is well-captured by the stellar particle information from TNG alone (and the stellar population synthesis model). Nevertheless, the Kannan2022 L (SFR) relations for these two lines are slightly higher than our predictions. This might be caused by the fact that we have assumed a different IMF for the simulated stellar populations. On the other hand, the TNG+HIILines L (SFR) relations for the [OIII] 88 and 52 micron lines are significantly higher than the GMDN results. This mainly results because the former method adopts the TNG gasphase metallicities directly, while the GMDN results are based on the FIRE simulation metallicities. Figure 2 shows that the TNG metallicities are higher than those in FIRE by as much as ∼ 1 order of magnitude, leading to the discrepant line luminosity predictions. Finally, the TNG+HIILines L (SFR) relations for the [OIII] restframe optical lines are also higher than the GMDN postprocessing results, although the difference is less drastic compared with the [OIII] sub-millimeter lines. This is because, with the higher TNG gas-phase metallicities, Eq 19 predicts lower ISM temperatures for each galaxy, and this suppresses the [OIII] optical line luminosities. We have not yet, however, found an explanation as to why the Kannan2022 [OIII] L (SFR) results differ from our TNG+HIILines predictions even at low SFR.', '5.3. Luminosity Functions': 'The line luminosity functions, ϕ ( L ), encode interesting information regarding ISM properties, stellar radiation fields, and also the population-level statistics of the emitting galaxies. More specifically, in the case of the Balmer lines, the line luminosity functions directly trace the overall ionizing emissivity of the emitting galaixes. Assuming photo-ionization equilibrium and case-B recombinations (which depend weakly on temperature), the rate of ionizing photon production (asides for those that escape the galaxy to ionize atoms in the IGM) may be translated directly into a Balmer line luminosity. Hence, after accounting for dust attenuation and the likely small average escape fraction of ionizing photons, the Balmer line luminosity function can be used to infer the rate of ionizing photon production per unit volume. This ionizing emissivity is a key quantity for understanding reionization, and so the Balmer line luminosity functions should provide important empirical guidance here. The [OIII] and [OII] line luminosity functions yield information regarding the chemical enrichment history of the universe, and the properties of the ISM, especially when compared to the Balmer line luminosity functions.  \n<!-- image -->  \nFigure 11. Correlation between each of the galaxy-wide hydrogen ionizing photon generation rate, Q HI (top panel), and the SFR-weighted metallicity, Z SFR (bottom panel), and SFR. The red, yellow, and blue lines show the median relations for the TNG50, TNG100, and TNG300 simulations, respectively, while the bands enclose the relations for 16% and 84% of the simulated galaxies. The best fits for the average Q HI (SFR) and Z SFR (SFR) relations are shown by the black curves. For comparison, Q HI (SFR) from the fitting formula of Schaerer (2003), within the stellar metallicity range of 0 . 1 ≤ Z/ [ Z ⊙ ] ≤ 1, is shown by the magenta band.  \n<!-- image -->  \nFurthermore, the relative luminosity functions of [OIII] and [OII] depend on the typical shape of the ionizing radiation spectrum and should thereby provide insights into the statistical properties of the stellar populations in the emitting galaxies. Furthermore, current and upcoming measurements, including new JWST observations, will deliver on the scientific promise here. In addition, line-intensity mapping surveys effectively probe integrals over the luminosity functions, including the impact of low-luminosity sources (e.g. Bernal & Kovetz 2022).  \nMotivated by these exciting prospects, the black curves in Figure 12 show z = 6 line luminosity function results from our new GMDN models on top of TNG50 (solid), TNG100 (dashed), and TNG300 (dotted). The combined luminosity functions are fit to Schechter func-  \ntion models:  \nϕ ( L ) = dn dL = ϕ ∗ ( L L ∗ ) α exp( -( L/L ∗ )) 1 L ∗ , (20)  \nwhere dn denotes the co-moving number density of galaxies in a given line luminosity bin, i.e. with line luminosity between L and L + dL . The best-fit results are given by solid curves in Figure 12 and the best-fit parameters are summarized in Table 3.  \nWe also compare the simulated luminosity functions with observations from De Barros et al. (2019); Matthee et al. (2023); Sun et al. (2023a). Our simulated [OIII] 5007 ˚ A luminosity function model agrees with observational results to within their 1 σ uncertainties in most luminosity bins, although the model lies consistently beneath the observational data. Likewise, the simulated Hα and [OIII]4960 , 5007 + H β luminosity functions fall below the current observational results by more than 1 σ . In fact, previous studies (without our sub-grid modeling framework) generally found more pronounced discrepancies between these measurements and models (e.g. Hirschmann et al. 2023).  \nOne possibility is that part of the discrepancies between luminosity function model and measurements owes to sample/cosmic variance. That is, the measurements are performed over small regions of the sky and so the cosmic variance uncertainties can be large, yet these are not generally incorporated into the observational error budgets. If the observed fields happen to be centered near over-dense regions, they could yield larger luminosity function estimates than expected based on our ensemble-averaged models. For example, the H β +[OIII] luminosity function measurement comes from the EIGER survey (Matthee et al. 2023). This survey spans 25.9 arcmin 2 over the redshift range from 5 . 33 < z < 6 . 96, corresponding to a volume of 1 . 06 × 10 5 cMpc 3 . This is similar to the volume spanned by the TNG50 simulation box. In order to assess the level of cosmic variance in the EIGER luminosity function measurements, we hence divide the TNG100 and TNG300 simulation boxes into 8 and 216 cubical sub-volumes to roughly match the EIGER survey volume. The spread in the luminosity function estimates across these subvolumes then provides an estimate of the cosmic variance contributions to the luminosity function uncertain-  \nes. 2 The minimum to maximum luminosity function ranges are shown as cyan bands in Figure 12. As can be discerned from the figure, the cosmic variance contributions to the error budget are important and can explain away much of the discrepancies with our models. Larger volume surveys will be required to enable stronger tests of these models in the future.  \nA virtue of our modeling framework is that it is efficient and flexible enough to allow us to explore variations around the input assumptions and thereby gain insights into which features of the models are most important. In order to obtain a better understanding of how the FIRE and TNG simulation properties influence our line luminosity function model predictions, we recalculate the luminosity functions after varying the assumed MZR model. As illustrated in Figure 2, the TNG MZR shows a higher normalization (i.e. larger metallicities at a given stellar mass) than in FIRE by about ∼ 0 . 4 dex. We hence explore the impact of manually boosting the metallicities among both the stellar and gas particles in FIRE by 0.4 dex. We then model the line emission from the modified FIRE HII regions, retrain the GMDN, and use this to predict modified postprocessed emission line signals from the TNG simulation. The resulting modified luminosity function models are shown as green curves in Figure 12. As expected, the [OIII] 88 and 52 µ m luminosity functions are shifted towards higher luminosities by a factor of ∼ 2 . 5 owing to the enhanced metallicities in the case of the modified MZR model. The luminosity functions in the rest-frame optical [OIII] lines also show luminosity enhancements, but the shifts are smaller in these lines. This occurs, because while increasing the metallicity of each galaxy boosts its luminosity, there is a counteracting decrease since the more metal rich gas leads to lower gas temperatures in our models. Finally, the Balmer lines are insensitive to variations in the MZR model as hydrogen lines are metallicity independent, and the only small effect arises through the metallicity-dependent shift in the gas temperature, which weakly impacts the Balmer line strengths.  \nOther factors that can influence the simulated luminosity functions include switching from the FSPS single star model to binary star radiation spectra, or adopting different dust attenuation treatments. We will explore  \nTable 3. Best fit parameters for the parameterized luminosity function Eq 20. The unit of ϕ ∗ is [cMpc ergs -1 s]. The unit of L ∗ is ergs/s. [OIII] + H β stands for the total luminosity of [OIII]4960, 5007, and the H β lines.  \nthese possibilities quantitatively in future work. Overall, it will be instructive to compare the upcoming HII region line luminosity function measurements with the simulation predictions: discrepancies should help refine the simulation models.', '6. SUMMARY AND DISCUSSION': 'Current and forthcoming emission line measurements from the JWST, ALMA, and SPHEREx allow new tests of state-of-the-art simulations of galaxy formation. This requires, however, producing detailed models of ISM emission lines across cosmological volumes. A difficulty here is that the line emission is shaped by the multi-phase properties of the ISM, yet this small-scale structure remains unresolved in large-volume cosmological simulations. To overcome this challenge, we introduced a novel, multi-scale modeling framework which combines state-of-the-art zoom-in galaxy simulations, a semi-analytic HII region emission line model, and the MDN machine learning architecture.  \nIn a previous work, Yang et al. (2023b), we applied a spectral synthesis model, HIILines , to 22 FIRE central galaxies at z = 6 to model the [OIII], [OII], and hydrogen recombination line luminosities among all of the individual HII regions across each simulated galaxy. The FIRE simulations provide self-consistent galaxyscale models at ∼ 10 pc resolution, partly capturing the properties of the multi-phase ISM. In addition, these 22 FIRE galaxies span a wide range in mass, metallicity, and morphology, and are therefore representative of a wide variety of the galaxies found in larger volume cosmological simulations. The strength of HIILines is its high computational efficiency: this allows it to account for variations in the stellar particle and local ISM properties, avoiding simplified treatments where the gas density and/or radiation field are assumed constant across  \nFigure 12. Comparisons between our GMDN model luminosity functions (black transparent curves) and current observations (De Barros et al. 2019; Matthee et al. 2023; Sun et al. 2023a). Schechter function fits to the simulation results are shown as solid black curves, and the corresponding best-fit parameters are summarized in Table 3. The red transparent curves show luminosity function models in which the FIRE galaxy stellar and gas particle metallicities are boosted by 0.4 dex. This model is intended to mimic the TNG MZR, in contrast to our fiducial case where the MZR matches that from the FIRE simulations.  \n<!-- image -->  \nall HII regions and galaxies. The galaxy-wide [OIII], [OII], hydrogen recombination line luminosities, line ratios, luminosity versus SFR relations, MZRs, and UV luminosity functions given by the FIRE highz simulations agree well with current measurements (Yang et al. 2023b; Sun et al. 2023b). The FIRE simulations hence serve as the basis for our sub-grid line emission models.  \nHere we used the 22 FIRE galaxies as a training set and the GMDN technique to extract joint cPDFs for the [OIII] [OII], H α , and H β line luminosities per stellar particle mass, conditioned on the stellar particle age, metallicity, and total galaxy stellar mass. In this trained GMDN model, the local ISM properties among the FIRE galaxies are effectively folded into the free parameters in the network. Assuming that independent sub-grid HII regions exist around each stellar particle in a larger volume cosmological simulation or semi-analytic model, the GMDN acts as a generative model to populate such sub-grid HII regions with line luminosities. It can be applied to any simulation or semi-analytic model which tracks stellar particles and has a coarser resolution than the FIRE simulations.  \nAs an example application, we use our GMDN approach to post-process HII region line signals on top of well-resolved galaxies in the TNG50, TNG100, and TNG300 simulations at z = 6. We present detailed comparisons between the models and observations as well  \nas with other modeling efforts, considering the galaxywide line emission signals, correlations between line luminosity and SFR, and galaxy line luminosity functions. We find that the multi-line GMDN L (SFR) relationships differ from those in other work in the current literature. The differences partly trace to variations in the smallscale ISM properties across different galaxies, which are captured in our scheme, but neglected in most other current studies. These differences highlight the necessity of a multi-scale simulation modeling framework for these lines, as important small-scale ISM property variations are otherwise neglected when predicting line emission statistics from large volume cosmological simulations. Although our GMDN model generally agrees with current observational measurements, it tends to underpredict the [OIII] 4364 ˚ A line luminosity. This line is highly sensitive to the HII region gas temperatures, and future versions of our model may benefit from improvements in our temperature models.  \nAlthough we find excellent agreement between simulated and observed line luminosity measurements, and line luminosity versus SFR relations, the simulated [OIII], H α , and H β line luminosity functions fall consistently lower than current measurements. We find that these discrepancies may owe to cosmic variance in the small fields spanned by current line luminosity function measurements.  \nIn summary, our machine learning based modeling framework is useful for exploring how galaxy formation model variations influence the observable ISM emission line statistics. By combining the galaxy population and stellar particle statistics sampled across large-volume cosmological simulations, with zoom-in galaxy simulations capable of resolving HII region properties, we avoid some of the simplifying assumptions made in previous work. We nevertheless aim to improve and extend our analysis in the near future. First, we hope to refine our treatment of the gas temperature in HII/OIII/OII regions. Second, our current study focused entirely on simulated galaxies at z = 6, and we hope to extend the FIRE+ HIILines calculations and GMDN modeling to lower redshifts to help interpret emission line obser-  \nvations at Cosmic Noon and in the more nearby universe. Third, it will be useful to extend HIILines to cover other HII region lines including OII recombination lines, [NII], [SII], and [SIII] lines, some of which will be observed by ALMA, the Roman Space Telescope, and multiple LIM surveys.', '7. ACKNOWLEDGEMENTS': "S.Y. thanks Jialian Li for beneficial discussions about generative machine learning methods. S.Y. acknowledges support from the Director's postdoctoral fellowship funded by LANL LDRD project No. 20240863PRD2. AL acknowledges support through NASA ATP grant 80NSSC20K0497. H.L. acknowledges the support by LANL/LDRD program."}