Simulating emission line galaxies for the next generation of large-scale structure surveys Open Access

ABSTRACT

We investigate emission line galaxies across cosmic time by combining the modified L-Galaxies semi-analytical galaxy formation model with the JiuTian cosmological simulation. We improve the tidal disruption model of satellite galaxies in L-Galaxies to address the time dependence problem. We utilize the public code cloudy to compute emission line ratios for a grid of H ii region models. The emission line models assume the same initial mass function as that used to generate the spectral energy distribution of semi-analytical galaxies, ensuring a coherent treatment for modelling the full galaxy spectrum. By incorporating these emission line ratios with galaxy properties, we reproduce observed luminosity functions for H α, H β, [O ii], and [O iii] in the local Universe and at high redshifts. We also find good agreement between model predictions and observations for autocorrelation and cross-correlation functions of [O ii]-selected galaxies, as well as their luminosity dependence. The bias of emission line galaxies depends on both luminosity and redshift. At lower redshifts, it remains constant with increasing luminosity up to around |$\sim 10^{42.5}\, {\rm erg\, s^{-1}}$| and then rises steeply for higher luminosities. The transition luminosity increases with redshift and becomes insignificant above z = 1.5. Generally, galaxy bias shows an increasing trend with redshift. However, for luminous galaxies, the bias is higher at low redshifts, as the strong luminosity dependence observed at low redshifts diminishes at higher redshifts. We provide a fitting formula for the bias of emission line galaxies as a function of luminosity and redshift, which can be utilized for large-scale structure studies with future galaxy surveys.

1 INTRODUCTION

In the past few decades, large-scale surveys have propelled revolutionary developments in astronomy. Various surveys, such as the Sloan Digital Sky Survey (SDSS; York et al. 2000; Gunn et al. 2006), the Dark Energy Survey (DES; Dark Energy Survey Collaboration 2016), and the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP; Aihara et al. 2018; Kawanomoto et al. 2018) have significantly enhanced our understanding of the universe and the underlying theories of galaxy formation. These surveys have contributed to the solidification of the ΛCDM cosmological framework, precise measurements of the expansion rate of the universe, and the provision of extensive data on the large-scale structure of the universe, including the mapping of millions of galaxies. In addressing evolving scientific challenges, currently ongoing and upcoming large-scale next-generation surveys, such as the Chinese Space Station Telescope (CSST; Zhan 2011), the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration 2022), Euclid (Laureijs et al. 2011), the Large Synoptic Survey Telescope (LSST; Ivezić et al. 2019), the Nancy Grace Roman Space Telescope (Roman; Dressler et al. 2012; Spergel et al. 2013), and the Subaru Prime Focus Spectrograph (PFS; Takada et al. 2014; Tamura et al. 2016) aim to achieve substantial breakthroughs in fundamental issues such as the origin and evolution of dark matter and dark energy, as well as the origin and evolution of galaxies and black holes. These state-of-art instruments and surveys have lower detection limits and better spatial resolution, thereby extending observational data towards fainter objects and higher redshifts.

Emission line galaxies (ELGs) serve as the main tracer at z ∼ 1 for large-scale surveys, playing a crucial role in estimating photometric and spectroscopic redshifts and properties related to galaxy formation. In recent years, numerous theoretical works have also focused on the formation of ELGs. Some studies have utilized the Halo Occupation Distribution (HOD) method to assign emission line luminosities to dark matter-only simulations (Geach et al. 2012; Avila et al. 2020; Gao et al. 2022; Reyes-Peraza et al. 2023; Rocher et al. 2023a, b). Approaches involving more physical models include semi-analytical models (Orsi et al. 2014; Merson et al. 2018; Stothert et al. 2018; Izquierdo-Villalba et al. 2019; Zhai et al. 2019; Favole et al. 2020; Baugh et al. 2022; Knebe et al. 2022) and hydrodynamical simulations (Park et al. 2015; Shimizu et al. 2016; Katz 2022; Smith et al. 2022; Tacchella et al. 2022; Hirschmann et al. 2023; Osato & Okumura 2023; Yang et al. 2023).

Numerical simulations (Kauffmann et al. 1999; Springel et al. 2001, 2005; Guo et al. 2011; Henriques et al. 2015; Schaye et al. 2015; Pillepich et al. 2018) have been successful in investigating contemporary galaxy formation and cosmology. These simulations have shown the ability to reproduce various observed galaxy characteristics across different periods in cosmic history. Their capability in predicting, interpreting, and optimizing observational outcomes has rendered them valuable tools for comprehending the processes underlying galaxy formation and the progression of large-scale structures. Hydrodynamic cosmological simulations such as EAGLE (Schaye et al. 2015) and IllustrisTNG (Pillepich et al. 2018) have shown enhanced capacity in studying baryonic physics processes, especially gas dynamics. Recently, Hirschmann et al. (2023) computed optical and NUV emission lines originating from various sources such as star clusters, narrow-line regions of AGN, post-asymptotic giant branch stars, and fast radiative shocks for galaxies in the IllustrisTNG simulation up to redshift 7 via post-processing, providing valuable predictions for JWST. Similarly, Osato & Okumura (2023) constructed mock H α and [O ii] catalogues based on IllustrisTNG to investigate the clustering of emission line galaxies. However, such approaches are hindered by substantial CPU time requirements and computational constraints, precluding the simultaneous achievement of both high precision and large simulation volumes.

By combining physically motivated recipes of galaxy formation with N-body cosmological dark matter simulations, semi-analytical models (e.g. Kauffmann et al. 1999; Cole et al. 2000; Guo et al. 2011; Lacey et al. 2016; Lagos et al. 2018) provide a cost-effective alternative. This approach allows simultaneous consideration of high precision and large simulation volumes, as compared to hydro simulations. Izquierdo-Villalba et al. (2019) combined the semi-analytical model l-galaxies (Henriques et al. 2015) with the emission line model from Orsi et al. (2014) to construct a light cone for the J-PLUS survey. Baugh et al. (2022) applied the pre-computed grid of emission line luminosity released by Gutkin, Charlot & Bruzual (2016) within the semi-analytical galaxy formation code GALFORM (Lacey et al. 2016) to reproduce the observed locus of star-forming galaxies on standard line ratio diagnostic diagrams. Favole et al. (2020) applied the emission line model described by Orsi et al. (2014) to three different semi-analytical models: SAG (Cora et al. 2018), SAGE (Croton et al. 2016), and GALACTICUS (Benson 2012) and concluded that utilizing average star formation rates is a feasible method to generate [O ii] luminosity functions. However, these works used different stellar population synthesis (SPS) models for computing stellar components and H ii regions, introducing additional inconsistencies in the final results.

In this study, we implement the state-of-art semi-analytic model l-galaxies (Henriques et al. 2015) onto the merger trees extracted from a large-box-size, high-resolution N-body dark matter simulation to produce a galaxy catalogue for upcoming large-scale surveys. We improve the satellite disruption model in l-galaxies to address a theoretical issue with varying time resolutions. We record the complete star formation history (SFH) for each individual galaxy, enabling the computation of photometric magnitudes for any given filters as post-processes. We combine a grid of H ii region models with the public radiation transfer code cloudy to derive emission line ratios using the same SPS model employed for the stellar components, ensuring the self-consistency of our predictions.This guarantees consistent treatment between the stellar SED and the emission line luminosities. The grid of H ii regions cover a wider parameter space compared to many previous work. By combining this with the semi-analytical galaxy output, we calculate the luminosities of the 13 most frequently utilized NUV and optical emission lines.

This paper is organized as follows. Section 2 provides an overview of the N-body dark matter simulations Jiutian-1G and Mini-Hyper and the details of our semi-analytic model and emission line models. Section 3 presents our model predictions for various galaxy properties, while Section 4 shows the properties of emission line galaxies. We conclude with a summary and discussion in Section 5.

2 DATA AND METHOD

In this section, we give a brief description about the dark matter merger trees, semi-analytic models, and emission line models. Two sets of N-body cosmological simulations are adopted, a large simulation, JiuTian-1G, and four sets of merger trees exacted from a small run with different numbers of snapshots. Details are listed in Table 1. Then we modify the l-galaxies model (Henriques et al. 2015) and apply it on the Jiutian-1G dark matter simulation. We then utilize a publicly available radiation transfer code to determine the luminosity of emission lines by implementing a photoionization model surrounding the star formation regions.

2.1 Dark matter simulation

JiuTian Simulations comprise a series of cosmological N-body simulations, ranging in box size and resolution. The JiuTian-1G (hereafter JT1G) simulation is a large dark matter only N-body simulation within the framework of Λ cold dark matter (ΛCDM) cosmology designed for next-generation surveys. Utilizing the L-Gadget3 code (Springel 2005), the JT1G simulation tracks 6144³ dark matter particles within a cubic simulation box with a side length L_box = 1 Gpc h⁻¹. This box length is twice as large as the previous Millennium Simulation (MS; Springel et al. 2005), resulting in an eight-fold increase in volume compared to MS. The particle mass is |$3.72 \times 10^{8} \, {\rm M_\odot }\, h^{-1}$|⁠, almost three times smaller than the original MS. The JT1G simulation stores 128 snapshots ranging from redshift 127 to 0, with an average time gap of approximately 100 Myr. This time resolution is chosen for weak lensing studies and is twice as large as MS with 64 snapshots. We adopt the cosmological parameters from Planck Collaboration I (2020) as follows: σ₈ = 0.8102, |$H_0=67.66{\rm \, km \, s^{-1}\, Mpc^{-1}}$|⁠, |$\Omega _\Lambda =0.6889$|⁠, Ω_m = 0.3111, Ω_b = 0.0490(f_b = 0.1575). Following Springel et al. (2005), dark matter haloes and subhaloes are identified using the friends-of-friends (FOF) and subfind (Springel et al. 2001) algorithms. Additionally, to establish the merger trees, the subhaloes are linked with their unique descendants employing the B-Tree code. Details about JT1G are referred to Han et al. (in preparation) and Li et al. (in preparation).

Although previous works have extensively examined the particle mass resolution (see Crain & van de Voort 2023, for a review), limited attention has been given to time resolution. Benson et al. (2012) concluded that 128 snapshots are necessary for the GALACTICUS semi-analytical model (Benson 2012) to achieve convergence within a 5 per cent level in stellar mass. We conduct four sets of merger trees with different temporal resolutions from a smaller simulation, Mini-Hyper, to assess the time convergence. The parent simulation has a box size of 125 |${\rm Mpc}\, h^{-1}$| and a particle mass similar to JT1G, |$3.674 \times 10^{8} \, {\rm M_\odot }\, h^{-1}$|⁠. The total particle number is 768³ which is stored in 513 distinct snapshots. We employ slightly different cosmological parameters compared to JT1G: σ₈ = 0.826, |$H_0=67.3{\rm \, km \, s^{-1}\, Mpc^{-1}}$|⁠, |$\Omega _\Lambda =0.693$|⁠, Ω_m = 0.307, Ω_b = 0.04825(f_b = 0.1572).

Four merger trees are then constructed accordingly with varying time intervals using the B-Tree code with different ‘SnapSkipFac’. This approach yields four sets of merger trees with comparable merger tree structures but differing numbers of snapshots. Table 1 shows the parameters for all our simulations. In practice, we fix the first and last snapshots across all simulations. Then, we select every second snapshot to construct a simulation with 257 snapshots. Further skipping every other snapshot results in a simulation with 129 snapshots. Following the same methodology, simulations with 65/33 snapshots were constructed, with the number of snapshots decreasing by a factor of 2 each time. Therefore, we obtain four sets of merger trees with similar tree structures, wherein the time intervals between two snapshots increase by a factor of 2 each time. It is worth noting that these four sets of merger trees have the same dark matter halo properties at common redshifts. Fig. 1 depicts how we create merger trees with different snapshots.

2.2 Semi-analytical model: l-galaxies

We use the version of the l-galaxies code as described in Henriques et al. (2015) (hereafter H15) and make modifications to solve the time convergence problem. H15 includes physical prescriptions for various baryonic processes, such as shock heating, gas cooling, star formation, supernova feedback, formation and growth of supermassive black holes (SMBH), AGN feedback, metal enrichment, etc. Details about the physical recipes and parameters can be found in their supplementary material. We have modified the satellite disruption procedure to address the issue of time convergence and have readjusted the parameters to replicate the abundance of SMBH in the local Universe.

2.2.1 Time convergence problem

We utilize the original H15 code to examine whether similar galaxy properties can be acquired on dark matter merger trees with varying time intervals. The black lines in the first row of Fig. 2 show the statistical properties of galaxies at z ∼ 0 predicted by the original H15 model, including the galaxy stellar mass function (SMF), galaxy abundance as a function of star formation rate (SFRF) and the supermassive black hole mass function (BHMF). Distinct line styles represent results from merger trees with different time intervals: solid lines for Mini-Hyper-33, dashed lines for Mini-Hyper-65, dotted lines for Mini-Hyper-129, and dotted–dashed lines for Mini-Hyper-257 merger trees. Surprisingly, we observe substantial differences in galaxy properties across simulations. The left panel reveals that simulations with smaller time intervals tend to have more massive galaxies. The difference is remarkably large, reaching several orders of magnitude at |$M_*\gt 10^{11}\, {\rm M_\odot }$| between Mini-Hyper-33 and Mini-Hyper-257. Larger offsets in SFRF and BHMF are evident, as the Mini-Hyper-257 showcases significantly higher numbers of highly star-forming galaxies and considerably less SMBH compared to the other simulations. These substantial differences with different time intervals strongly suggest a challenge in time convergence within the H15 code.

2.2.2 Relevant physical processes

Further investigation into the H15 code reveals that the primary cause of the time convergence problem is linked to the fate of orphan satellite galaxies, which lost their subhaloes due to physical processes or numerical effects. In H15, following the disruption of its subhalo, a merging clock is set simultaneously based on an estimated time (t_friction) for the orphan to spiral into the central galaxy due to dynamical friction. An orphan galaxy will ultimately either undergo disruption or merge. It could suffer tidal disruption on its way spiralling into the centre, relying on the competition between self-gravity and the tidal force from the main dark matter halo. In practice, a comparison is made between the baryonic (cold gas and stellar mass) density within the half-mass radius and the dark matter density of the main halo at the assumed pericentre (R_peri) of the orphan’s orbit.

If the tidal force exceeds the bounding gravity, prior to Δt = t_friction where Δt is the time since merger clock is set, the orphan galaxy will be completely disrupted and no longer undergo merging. Their SMBH, gas, and stars are all distributed in the halo of the central galaxy. SMBH in the central galaxy will not grow in this scenario. Conversely, if tidal force never exceeds the bounding gravity (until Δt = t_friction), the orphan galaxy will merge into the central galaxy at Δt = t_friction. During galaxy mergers, the SMBH in the central galaxy could grow by swallowing the SMBH from the satellite galaxies that are being merged and by experiencing strong gas accretion, which is triggered by the merger process.

A residual time-dependent treatment was adopted in the original H15 model, where the merger recipe and the satellite disruption recipe are called at different times. H15 divides the time between two adjacent snapshots into 20 sub-steps and calls the ‘merger’ recipe in each sub-step, while the ‘disruption’ recipe is only called at the end of each snap gap. This prioritizes the occurrence of merger, while the disruption is delayed until the final sub-step, regardless of meeting the disruption criterion earlier. Larger intervals of time between two consecutive snapshots increase the chances of mergers occurring. As a result, there are fewer disrupted orphan galaxies, and the SMBH become more massive. The increased mass of SMBH in turn leads to more effective active galactic nucleus (AGN) feedback, resulting in smaller central galaxies. We conducted experiments to preserve time resolution in the substep by adjusting the number of substeps in different Mini-Hyper simulations, but this did not solve the time resolution problem.

2.2.3 Modifications to the disruption model

To address this time convergence problem, we make modifications to the disruption model in the code. We replace the assumed peri-centre distance with the actual distance of orphan satellite galaxies from the central galaxy R_orphan, and call the disruption function at each substep. We track the most bound particle, which was defined at the point just before the orphan galaxy lost its substructure. By multiplying its distance from the central galaxy R_mostbound by a factor that accounts for the impact of dynamical friction, we can estimate the distance of the orphan galaxy R_orphan as follows.

We introduce a minimum distance at which disruption could occur, set as the scale radius of the gas disc in the central galaxy, R_gas. The choice of R_gas as the minimum distance is justified by the notion that if an orphan galaxy has reached the region of R_gas, it should be considered as having entered the region of central galaxies and is undergoing a strong interaction. In such situations, it is more appropriate to treat the event as a merger rather than a disruption.

We implement the modified model to the Mini-Hyper simulations with different snapshots without making any changes to the parameters from the initial H15 version. The stellar mass function, star formation rate function, and black hole mass function of the modified model are shown as the red lines in the first row of Fig. 2, indicating good agreement among different time resolutions. The quantification of the differences with respect to the Mini-Hyper-257 simulation is presented in the second row and shows that the differences in SMF and SFRF are within 5 per cent in most cases. At higher values, it could suffer from the limited sample size. The variation in BHMF seems larger, especially at the low-mass and at the high-mass end. This suggests that the growth of SMBH could be more sensitive to the time interval between consecutive snapshots. Furthermore, varying numbers of snapshots could also result in difference within the generated trees. For example, Wang et al. (2016) showed that using more snapshots typically results in shorter branches for most tree builders. This is because more linking errors may occur when processing more snapshots, as tree builders are prone to resolution or flip-flop problems. Han et al. (2018) also showed that the subfind and dtree combination could produce a substantial amount of fragmented branches, which in turn impacts the properties of resulting galaxies.

2.2.4 Parameter adjustment

It has been noticed that the mass of SMBH is lower by an order of magnitude compared to current observations (e.g. Yang, Hu & Li 2019) at z = 0 (see also the upper panel in Fig. 9, the black line is the result of H15, the purple symbols are observational data). To determine the parameters of our new model, we incorporate the black hole mass function at z = 0 as constraint, in addition to the observed galaxy stellar mass function at z = 0, 1, 2, and 3, and the passive fraction at z = 0.4. The observational data we use in this study are listed in Table A1. In accordance with Henriques et al. (2013), we establish a representative subset of subhalo merger trees and employ the MCMC method as described in Henriques et al. (2009), to thoroughly explore the multidimensional parameter space with the updated disruption model in l-galaxies. In short, we divide the haloes into I halo mass bins with a width of 0.5 dex, and the galaxies into J stellar mass bins with a same width. We randomly select n_i haloes from a total N_i haloes in the ith halo mass bin. The number n_i is determined by a set of linear constraint equations, |$\sum _{i=1}^I \frac{N_i^2}{n_i^2} \phi _{i j}=F^2 \Phi _j^2$|⁠, where ϕ_ij is the average number of galaxies in the jth stellar mass bin for haloes in the ith halo mass bin, Φ_j is the total number of galaxies in the jth stellar mass bin. F is the uncertainty of the stellar mass function, which we set to be F < 0.05. Our final representative subsample is ∼ 1/512 of the whole box. Details about the sample construction is available in Henriques et al. (2013), APPENDIX B.

For comparison with the default H15 parameters, the best-fitting parameters are enumerated in Table 2. Our modified code maintains compatibility with the default H15 parameters. The best-fitting parameters are then applied to the full volumes of the JT1G simulation to generate the SAM galaxy catalogue.

2.2.5 Mass-resolution convergence test

We have evaluated the mass resolution effect by applying the same SAM models (our modified model) and parameters (the best-fitting parameters for JT1G) to the merger trees extracted from the rescaled Millennium-II simulation (Boylan-Kolchin et al. 2009; Henriques et al. 2015). The MSII simulation has a mass resolution of |$7.69\times 10^{6}\, {\rm M_\odot }\, h^{-1}$|⁠, which is ∼ 50 times higher than JT1G, although the volume is smaller by 1000 (see Henriques et al. 2015, and reference therein for more details about the rescaled MS ii). MSII has been shown to resolve the smallest halo capable of hosting a detectable galaxy (e.g. Guo et al. 2011). Fig. B1 illustrates that the stellar mass of the galaxy converges at |$10^{9}\, {\rm M_\odot }$| between JT1G and MSII within a 10 per cent from z = 0–3. This difference increases to 30 per cent at |$10^{8}\, {\rm M_\odot }$|⁠, partly attributed to slight differences in the cosmological parameters and possibly stemming from distinctions in their initial conditions (Li et al. in preparation). At high masses, the disparity between JT1G and MSII is mainly due to varying mass resolutions. The mass resolution in MRII is about 50 times higher than in JT1G, allowing it to resolve more small haloes that can later merge into larger systems. These mergers lead to larger supermassive black holes (SMBHs) in MRII, resulting in more effective AGN feedback that suppresses star formation and ultimately leads to a smaller stellar mass. The difference in stellar mass in massive halos is not significant, typically around 0.15 dex. However, due to the steep slope at the high mass end of the SMF, this slight variation in stellar mass can result in a notable difference in the stellar mass function at high masses. Along with mass resolution, cosmic variance stemming from the smaller volume of MSII could also be a factor. Similar comparisons are performed on the abundance of the galaxy as a function of SFR in Fig. B2. It shows a very good convergence between the JT1G simulation and MSII at z = 0. At higher redshifts, the convergence is somehow larger, but all within 10 per cent.

2.3 Emission line model

In this section, we first briefly describe the process of generating the galaxy spectral energy distribution (SED). Then, we explain in detail how to use the expected SED as the input ionizing spectrum to calculate the luminosities of emission lines. In practice, we employ the radiative transfer code cloudy to calculate the relative strength of emission lines on a grid of H ii region models. Meanwhile, we adopt empirical relations to establish connections between the general properties of galaxies and the parameters that describe the ionization regions. According to the general properties of a given galaxy, the luminosity of each emission line is derived by performing interpolation within a pre-computed grid of line luminosities. During this process, we utilize the same Stellar Population Synthesis model and initial mass function as those employed in calculating the galaxy SED, ensuring the self-consistency of galaxy SED and emission line calculation.

2.3.1 Galaxy spectral energy distributions

Stellar population synthesis (SPS) models serve as essential tools in astrophysics, enabling the generation of synthetic SEDs with detailed information about the star formation history (SFH), metallicity, and initial mass function (IMF). These models provide valuable insights into the formation and evolution of galaxies. In this work, our default SPS model is based on the Bruzual & Charlot (2003) framework, with a Chabrier IMF (Chabrier 2003).

The entire SFH of both the disc and bulge components is stored in 22 distinct bins, as detailed by Shamshiri et al. (2015). Consequently, the SED can be produced using any desired SPS model and IMF as post-processes.

2.3.2 Grid of emission line models

We use the c17.04 release of the photoionization code cloudy (Ferland et al. 2017), a widely used tool, to address radiation transfer in photoionized regions. To solve the radiation transfer equation, we specify the most important input physical properties of the gas cloud and spectrum of ionizing sources, including the intensity and spectrum of the ionizing photons, gas geometry, gas metallicity Z, and hydrogen density n_H.

We adopt the BC03 (Bruzual & Charlot 2003) SPS model with a Chabrier IMF (Chabrier 2003) to calculate the intensity and spectrum of ionizing photons. These SPS model and IMF are also used in calculating the galaxy SED. It is noteworthy that previous studies did not consistently employ the same stellar population models for computing the galaxy SED and for the photoionization modelling of nebular emission. For example, Baugh et al. (2022) used the M05 model (Maraston 2005) to calculate galaxy stellar SED, yet adopted the emission line models based on BC03 in the H ii regions. Izquierdo-Villalba et al. (2019) adopted different stellar synthesis models to calculate the emission line grid (Levesque, Kewley & Larson 2010) and the stellar SED (Bruzual & Charlot 2003) when generating the emission line galaxy mock catalogue for J-PLUS. This difference results in a lack of self-consistency between the two distinct components of the combined galaxy spectrum (i.e. stellar and emission). In contrast, our work ensures internal consistency by employing the same SPS model (BC03) and IMF (Chabrier) for both the stellar component and the nebular model within galaxies.

For the gas geometry, we follow Byler et al. (2017) to assume a spherical shell geometry with the ionizing source located in the centre, and to adopt an inner radius R_inner at 10¹⁹ cm.

In contrast to Byler et al. (2017), who assume a constant gas density, we take into account the gas density within various ranges. In practice, we calculate the emission line ratios for a grid of H ii region models by sampling |$\mathcal {U}$|⁠, Z, and n_H as follows:

where |$\mathcal {U}$| is the ionization parameter, a dimensionless parameter defined as the ratio of hydrogen ionizing photons to the total hydrogen density

where n_γ is the volume density of the ionizing photons. Z is the metallicity of cold gas, and n_H is the hydrogen density. The output is a grid of line ratios relative to the H α luminosity with different parameters. Using linear interpolation in |$\log _{10} \mathcal {U}$|⁠, log₁₀Z/Z_⊙ and log₁₀n_H, we extract line ratios for each star-forming galaxy. If the values for |$\mathcal {U}$|⁠, Z, and n_H fall outside the range of the grid, we use the closest limiting values from the grid to align with the prediction of the catalogue.

Fig. 3 shows the classic BPT diagram (Baldwin, Phillips & Terlevich 1981) of the pre-computed grid of emission line luminosities. The x-axis depicts |${\rm [N\, {\small II}]\lambda 6584/H\,\alpha }$|⁠, and the y-axis depicts |${\rm [O\, {\small III}]\lambda 5007/H\,\beta }$|⁠. The grey grid with blue symbols shows the grid with a hydrogen density of n_H = 10 cm⁻³, varying the metallicity Z and the ionization parameter |$\mathcal {U}$|⁠. Additionally, we present the grid of emission line ratios from Byler et al. (2017) as black lines with red symbols. We employ a hydrogen density, metallicity, and ionization parameter grid akin to that utilized in Byler et al. (2017), albeit with a different spectrum of ionizing photons. Furthermore, we use the Chabrier IMF, while they employ the Kroupa IMF (Kroupa 2001). Our BPT diagram bears a resemblance to theirs, with a slight offset observed at higher metallicity. The |${\rm [O\, {\small III}]\lambda 5007/H\,\beta }$| ratio exhibits a monotonically increasing trend with increasing |$\mathcal {U}$|⁠, while the |${\rm [N\, {\small II}]\lambda 6584/H\,\alpha }$| ratio demonstrates a monotonically increasing pattern with increasing Z. We notice that as Z becomes sufficiently large, an overlapping within the grid itself becomes evident, reflecting the degeneracy of Z and |$\mathcal {U}$|⁠.

2.3.3 H α luminosity

cloudy provides line ratios relative to the H α luminosity, which is closely related to the star formation rate. For a nebula that is absolutely optically thick for ionizing photons and optically thin for redward photons (case B, Osterbrock & Ferland 2006), we can theoretically derive the relation between the intensity of a specific hydrogen recombination line and the ionizing photon rate by quantum mechanics. The relation between the luminosity of H α and the ionizing photon rate is expressed as

where |$\alpha _{\mathrm{H} \,\alpha }^{\mathrm{eff}}$| is the effective recombination coefficient at H α, and α_B is the case B recombination coefficient. Q_H is the total number of ionized photons emitted per second, which could be related to the star formation rate assuming the Chabrier IMF (Chabrier 2003):

The H α luminosity can then be expressed as a function of SFR (Falcón-Barroso & Knapen 2013):

It is noteworthy that we use the same IMF in this L(H α) − SFR relation as employed in both the stellar component and the photoionization model, which ensures the overall consistency of our model.

When paired with the line-ratio produced by cloudy as outlined in Section 2.3.2, we are able to calculate the luminosity of any specific emission line at λ_j.

where R(λ_j, q, Z, n_H) is the cloudy predicted ratio of the desired emission line at wavelength λ_j to H α with a given set of |$\mathcal {U}$|⁠, Z and n_H.

We include the 13 most widely emission lines in NUV and optical ranges in the final catalogue,¹ as listed in Table 3. Additional emission lines can be provided upon request.

2.3.4 Emission line luminosity for semi-analytical galaxies

With the grid of emission lines, the final step in predicting the emission lines for semi-analytical galaxies is to establish a connection between the general properties of the galaxies and the parameters that determine the emission lines. These parameters include gas metallicity, hydrogen density, and ionization parameters.

The gas metallicity is obtained directly from the semi-analytical catalogue. Here, we use gas metallicity Z_cold for the whole galaxy. Calculations including individual heavy elements will be performed in further work.

We adopt the empirical relations based on local observations (Kashino & Inoue 2019) to link the hydrogen density n_H to the stellar mass and specific star formation rate as follows,

where sSFR is the specific star formation and M_* is the stellar mass. Both can be obtained directly from the semi-analytical catalogue.

The joint effect of the ionizing spectrum and its intensity can be described by the ionization parameter, |$\mathcal {U}$|⁠. We adopt the empirical relations based on local observations (Kashino & Inoue 2019) to link the ionization parameter to the specific star formation rate, gas metallicity and hydrogen density as follows,

So far, we have obtained the input parameters for the corresponding galaxy. By interpolating within a pre-computed grid of line luminosities, as explained in the previous section, one can determine the luminosity of every emission line for each model galaxy.

The procedure for generating emission lines is generally similar to previous works (Orsi et al. 2014; Izquierdo-Villalba et al. 2019; Favole et al. 2020) but varies in details. Our pre-calculated emission line models rely on gas density, metallicity, and ionization parameter. These previous works assumed a fixed |$n_{\rm H} = 100 \rm cm^{-3}$|⁠, with the photoionization parameter depends on gas metallicity. Consequently, their emission line models are solely influenced by metallicity. Furthermore, we employ the same BC03 SPS model and Chabrier IMF for determining both stellar SED and emission line luminosities, thus ensuring internal consistency within the model. Orsi et al. (2014), Izquierdo-Villalba et al. (2019), and Favole et al. (2020) utilize the Salpeter IMF for computing emission line ratios but switch to the Kroupa IMF for calculating the H α luminosity. Baugh et al. (2022) introduced a distinct SPS model for stellar SED (M05) and nebular regions (BC03). The main difference between our approaches and Merson et al. (2018) lies in the use of different links between galaxy properties and the input parameters for cloudy. For example, they use an average gas density, the total gas mass divided by the cubic radius, whereas we implemented a empirical scaling relation associating the SFR and stellar mass with typical gas density in the SFR vicinity.

2.4 Dust extinction

The extinction induced by dust significantly influences the observed spectra of galaxies, absorbing UV/optical photons and re-emitting them at longer wavelengths. Consequently, galaxies rich in dust tend to have red colours even if they have a high SFR. In this work, we follow the dust model in H15, considering both the extinction from diffuse ISM (Devriendt, Guiderdoni & Sadat 1999) and from star-forming molecular clouds (Charlot & Fall 2000). The optical depth of dust, as a function of wavelength, is independently computed for each component; and a slab geometry is assumed to compute the total extinction of the relevant populations.

First, we calculate the extinction caused by the diffuse ISM. The wavelength dependent optical depth of galaxy disks is assumed as follows:

where |$\left(\frac{A_\lambda }{A_v}\right)_{Z_{\odot }}$| is the extinction curve for solar metallicity taken from Mathis, Mezger & Panagia (1983), (1 + z)⁻¹ represents the redshift dependence, Z_gas is the metallicity of cold gas, s = 1.35 for λ < 2000 Å and s = 1.6 for λ > 2000 Å. 〈N_H〉 is the mean column density of hydrogen

where R_{gas, d} is the scale-length of the cold gas disk, and a = 3.36. It should be noted that in previous work (Guo et al. 2011; Henriques et al. 2015), although they claimed to adopt a = 1.68, the actual factor used in their code was approximately a ∼ 3.36. Therefore, we use 3.36 instead of 1.68 in our calculations.

Another source contributing to the extinction is the molecular cloud around young stars. Following Charlot & Fall (2000), we assume that only young stars born within the last 10 Myr will suffer from such effects. The optical depth is calculated as follows

where μ is given by a Gaussian distribution with a mean of 0.3 and a standard deviation of 0.2, truncated at the boundaries of 0.1 and 1.

The final extinction in magnitude for each component is given by:

where θ is the inclination angle between the angular momentum of the disc and the z-direction of the simulation box, and τ_λ is the optical depth of the corresponding component. Young stars (age less than 10 Myr) suffer from both extinction components, while older stars are affected by diffuse ISM only. In the case of emission lines, we only consider the extinction from molecular clouds, as we only calculate the emission of star-forming regions.

In the following sections, all results are calculated using the new model and parameter settings applied to the merger trees from JT1G, unless otherwise specified.

3 GENERAL GALAXY PROPERTIES

This section presents a comprehensive comparison of various galaxy properties between our model predictions and observational data. We categorize the comparison into two classes: one for properties utilized as constraints in our MCMC parameter adjustment, and properties directly related; and the other for properties that are not utilized as constraints.

3.1 Comparison with observational galaxy properties used as MCMC constraints

3.1.1 Galaxy stellar mass functions

Fig. 4 shows the galaxy stellar mass functions span redshifts from 0 to 3. The red lines depict our model results, while the black lines represent the results of MS using default H15 model. The purple symbols with error bars are the observational data points used in our MCMC procedures that were produced by combining various observational studies in an attempt to estimate systematic uncertainties in the constraints. Further details can be found in Appendix A of H15. Our predicted SMF aligns successfully with observational data across the entire stellar mass range, spanning from local to high redshift. Notably, at very low masses, our results exhibit a slightly steeper slope compared to observations. Both the new set of parameters and the higher resolution of the JT1G dark matter simulation compared to MS, where H15 was conducted, could potentially contribute to the observed differences. Our results outperform H15 for stellar masses beyond the knee of the SMF. This is because we have a much larger volume, eight times larger to be exact, which enables us to include more galaxies with high masses.

3.1.2 Star formation rate

The star formation rate (SFR) represents a fundamental statistical measure within galaxy formation theory. In this work, a passive galaxy is defined as having an sSFR less than |$10^{-11}\rm yr^{-1}$|⁠, where sSFR = SFR/M_* denotes the specific star formation rate. We utilize the passive fraction at z = 0.4 as a constraint when adjusting the parameters.

Fig. 5 displays the probability distribution function (PDF) of sSFR for different stellar mass bins at redshift 0. In accordance with H15, we assign a random Gaussian sSFR centred at log(sSFR) = −0.3log (M_*) − 8.6 with a dispersion of 0.5 for model galaxies with |${\rm log}({\rm sSFR}) \le -12 \rm yr^{-1}$|⁠. The shaded regions represent results from SDSS DR7 (Brinchmann et al. 2004; Salim et al. 2007), while the black lines are results from H15. Our new model predicts an sSFR distribution similar to H15 across the entire stellar mass range and shows a higher peak in sSFR for the subset of star-forming galaxies at low masses, which is more consistent with the SDSS data compared to H15.

In Fig. 6, we show the SFRF at redshift 0, with various observational data points (Mauch & Sadler 2007; Robotham & Driver 2011; Patel et al. 2013; Gruppioni et al. 2015; Marchetti et al. 2016; Zhao et al. 2020). It is noteworthy that there exists a significant discrepancy between different observational data sets, with variations of up to three orders of magnitude at the high SFR end. Our model results fall within the range covered by these observations, emphasizing the need for more accurate measurements to better constrain theoretical models.

3.1.3 Red and blue galaxies

sSFR is closely correlated with galaxy colours, which are widely employed to distinguish the star formation states of galaxies. Therefore, we further study the stellar mass functions of red galaxies and blue galaxies at z = 0 and 0.4. We use the same colour cut as H15 for segregating galaxies into red and blue by: u − r = 1.85 − 0.075 × tanh((M_r + 18.07)/1.09). The upper and bottom panels in Fig. 7 show the SMF of red and blue galaxies, and the left and right columns present the findings at redshifts 0 and 0.4, respectively. Both our model results and those of H15 successfully reproduce the overall shape of the observations (Bell et al. 2003; Baldry et al. 2004; Ilbert et al. 2013; Muzzin et al. 2013). However, our model has slightly more red galaxies at the low mass end at z = 0. Furthermore, the SMF of red galaxies experiences a noticeable decline at |$M_*\sim 10^{10}\, {\rm M_\odot }$| for H15 at z = 0, which is absent in our model.

The discrepancy between model predictions and observations becomes more pronounced when expressed in terms of the red fraction as a function of stellar mass (see Fig. 8). The upper and bottom panels illustrate the red fraction as a function of stellar mass at redshift 0 and 0.4, respectively, with observational data points sourced from Bell et al. (2003), Baldry et al. (2004), Muzzin et al. (2013), and Ilbert et al. (2013). At redshift 0, our model’s results align marginally with the observations. We overestimate the proportion of red galaxies at the lower mass end and underestimate the red fraction at intermediate masses. H15 performs better at low masses but deviates more from the observations at |$M_*\sim 10^{10}\, {\rm M_\odot }$|⁠. At redshift 0.4, both our model and H15 are in line with the observations, although at intermediate masses, our model exhibits slightly fewer passive galaxies compared to H15. These findings underscore the need for further investigation into the quenching mechanisms, particularly around the knee of the stellar mass function.

3.1.4 Supermassive black holes

It has been noted that H15 underestimates the mass of SMBH by approximately an order of magnitude (Wang et al. 2018; Yang et al. 2019), suggesting that in the default H15 model the black hole growth rate is too slow. As illustrated in the upper panel of Fig. 9, our model agrees well with the observed BHMF taken from Shankar et al. (2009), which is used in the MCMC procedures, while H15 significantly underestimates the BHMF by an order. The new parameter set significantly enhances the growth rate during mergers, while also reducing the efficiency of AGN feedback. This adjustment allows the SMBH to attain larger masses, suppressing star formation activities in massive galaxies. Modifications in these aspects were essential to achieve a more precise representation of the observed BHMF.

The bottom panel of Fig. 9 illustrates the SMBH mass versus bulge mass relation. Observational data are obtained from Häring & Rix (2004), and the pink shaded region represents the relation from Kormendy & Ho (2013). Our new model predicts larger SMBH masses at a given bulge mass than H15 and aligns more closely with the observations. This alignment is a direct consequence of our enhanced growth rate during mergers.

Some earlier studies also focus on improving the modeling of SMBH growth in l-galaxies. Izquierdo-Villalba et al. (2020) introduces a time delay for Eddington accretion, promotes galaxy mergers, and incorporates additional pathways for SMBH growth, like disc instabilities. Spinoso et al. (2023) improve the modelling of SMBH seeds through various formation channels. Izquierdo-Villalba et al. (2024) construct a new framework by including the super-Eddington accretion events. All these advancements have resulted in better alignment of SMBH in l-galaxies with observational data.

3.2 Comparison with observational galaxy properties not used as MCMC constraints

3.2.1 Galaxy versus dark halo relations

The relationship between galaxy and dark halo properties represents one of the most fundamental connections in the field of galaxy formation. Fig. 10 illustrates these scaling relations by highlighting the correlation between galaxy stellar mass and both the maximum halo velocity and the maximum virial mass. The maximum velocity and virial masses were determined at z = 0 for central galaxies, while for satellite galaxies, they were determined at the last infall time. In general, both of these scaling relations demonstrate good agreement with the H15 models, particularly in terms of their slopes. With regard to the correlation between stellar mass (M_*) and virial mass (M_vir), the newly proposed model showcases somewhat enhanced stellar masses for JT1G galaxies in comparison to the H15 models. This suggests a higher galaxy formation efficiency in the new model. It is noteworthy to mention that despite the minor variations, both model forecasts agree with direct measurements obtained using local data (Ristea et al. 2024) and abundance matching methodologies (Moster et al. 2018; Behroozi et al. 2019).

3.2.2 Gas metallicity

Fig. 11 depicts the relationship between gas metallicity and stellar mass. The metallicity of the gas (Z_gas) is determined by the ratio of the metal mass in the cold gas to the mass of the cold gas. Then it is converted to oxygen abundance using the formula 12 + log₁₀(O/H)_gas = 8.69 + log₁₀(Z_gas/Z_⊙). It is evident that our newly developed model exhibits somehow higher gas metallicity compared to H15 across a wide range of masses. This characteristic enables it to be more in line with the observations (Tremonti et al. 2004) made at high masses.

3.2.3 Evolution of SFR

In Fig. 12, we compare the model-predicted star formation rate density (SFRD) as a function of redshift to various observational results (Sanders et al. 2003; Takeuchi et al. 2003; Schiminovich et al. 2005; Wyder et al. 2005; Dahlen et al. 2007; Reddy & Steidel 2009; Karim et al. 2011; Magnelli et al. 2011; Robotham & Driver 2011; Bouwens et al. 2012; Cucciati et al. 2012; Gruppioni et al. 2013; Magnelli et al. 2013; Schenker et al. 2013) from redshift 0 to 6. The grey line represents the best-fitting result from Madau & Dickinson (2014). The observed SFRD peaks at a redshift of 2 gradually, which diminishes in magnitude when moving towards either higher or lower redshift values. The results of our simulation generally reproduce this trend, with specific values falling within the observational constraints. This suggests that our model effectively traces the evolutionary history of the SFR in the universe. More stars are formed in the new model on JT1G compared to H15 on MS below redshift 4.

3.2.4 Galaxy correlation functions

The study of galaxy correlation functions encompasses the analysis of the spatial distribution of galaxies, which in turn provides valuable information on the distribution of matter. We use the Landy–Szalay estimator (Landy & Szalay 1993; Szapudi & Szalay 1998) to calculate the redshift-space cross-correlation functions between different subsamples

where r_p is the distance along the projected direction and r_π is distance along the line-of-sight. D_xD_y, D_xR_y, D_yR_x, and R_xR_y are the normalized pair counts for data x-data y, data x-random y, data y-random x and random x-random y. Here, x and y indicate different subsamples. When x = y, it becomes an autocorrelation function. The real-space projected correlation function |$w_{\mathrm{p}}^{x y}\left(r_{\mathrm{p}}\right)$| is calculated by integrating ξ_xy(r_p, r_π) along the line-of-sight (Davis & Peebles 1983)

where r_{π, max} = 40 |${\rm Mpc}\, h^{-1}$|⁠. We adopt the same r_p and r_π values as in Li et al. (2006) to facilitate a direct comparison of our results with theirs. Specifically, we employ 28 r_p bins spanning from 0.1–50 |$\rm Mpc\,\mathit{ h}^{-1}$| with equal logarithmic intervals and 40 r_π bins covering from 0–40 |$\rm Mpc\,\mathit{ h}^{-1}$| with equal linear intervals. The Corrfunc code (Sinha & Garrison 2020) is utilized for the calculation of the correlation function.

Given that galaxy clustering can be strongly influenced by their mass, we have categorized our model galaxies into various bins based on their stellar masses. Fig. 13 illustrates the comparison between the projected autocorrelation functions of our simulated galaxies and those observed in SDSS DR7. Each panel within the figure represents a distinct stellar mass range, shown in the upper right corner, with the black lines representing the results for the entire galaxy sample. Notably, the results for red and blue galaxies are depicted by red and blue lines, respectively. The symbols accompanied by error bars in the figure signify the measurements obtained from Li et al. (2006).

Our catalogue effectively reproduces the projected autocorrelation function across the majority of stellar mass ranges while also capturing its dependence on colour. However, certain discrepancies between the simulation and the observations are evident. Specifically, at larger radii, we observe deviations for galaxies smaller than |$10^{10}\, {\rm M_\odot }$|⁠. This could potentially be attributed to the limited volume occupied by these faint galaxies, thereby experiencing a notable cosmic variance effect.

Furthermore, we note differences between the predictions of our model and the observations at smaller radii for galaxies exceeding |$10^{10.5}\, {\rm M_\odot }$|⁠, indicating a significant proportion of satellite galaxies. By segregating red galaxies from blue galaxies, we find that the more strongly clustering of red galaxies contributes more prominently to the overall excess observed at smaller radii.

In summary, our refined model, incorporating adjusted parameters, excels in reproducing a wide range of observational properties: the stellar mass functions from local to redshift 3, the bimodal colour distributions of galaxies, the black hole mass functions, as well as other galaxy–halo relations and the clustering of galaxies, etc.

4 PROPERTIES OF VARIOUS EMISSION LINES

In this section, we first demonstrate that our galaxy catalogue meets the requirement for the next generation of large-scale surveys. We then present the model predictions of emission lines in comparison with observations. In Section 4.2, we show the luminosity function of H α, [O ii], [O iii], and [O iii] + H β, comparing our model results with a set of observations from local to high redshift. In Section 4.3, we further compare the projected correlation function of |$L_{\rm [O\, {\small II}]}$|-selected samples with observations. Finally, in Section 4.4, we show the model predicted bias of galaxies as a function of luminosity for different emission lines.

4.1 Luminosity functions for the next generation of large-scale surveys

Predictions for the luminosity function (LF) of various surveys can provide more direct forecasts for future observations. We convolve the filter functions of various surveys with the simulated galaxy SED to obtain the observed-frame LF in different bands in Fig. 14. It includes the luminosity functions in the Euclid Y, J, and H bands, in the LSST u, g, r, i, z, and y bands, and in the CSST NUV, u, g, r, i, z, and y bands across redshifts from 0 to 3. All magnitudes account for both the two components of dust extinction mentioned in Section 2.4 and include the contributions from various emission lines. The bottom x-axis displays apparent magnitude, while the top x-axis represents absolute magnitude. The two vertical dashed grey lines on the figures signify the detection limits of the CSST main survey (apparent magnitude of 26.5) and the deep field survey (apparent magnitude of 28). The transition point in the luminosity function, where the number density starts to decrease towards the dimmer end, serves as an approximate indicator of completeness. It happens at about 26.5 mag at z = 0.24, and around 28 mag at z = 0.5. Our simulation is deemed complete for the main survey above a redshift of 0.3 and the deep field survey above a redshift of 0.5. Similar completeness is expected for other surveys. We note that most tracers are above z ∼ 0.5. Meanwhile, the large box size of JT1G (1 Gpc h⁻¹) makes it suitable for studying large-scale structures. Therefore, our results offer a reasonably comprehensive prediction for the outcomes of the next generation of large-scale surveys.

4.2 Luminosity Function of Various Emission Lines

Fig. 15 shows the luminosity function of H α, [O ii], [O iii], and [O iii] + H β from redshift 0 to 3. The model predictions incorporate dust attenuation from young birth cloud as mentioned in Section 2.4 to align with direct observational symbols obtained from various sources. For clarity, we selected four representative redshifts to display for each line, ensuring the inclusion of corresponding observational data. Each observational point is presented in panels corresponding to the nearest redshift values.

The first row in Fig. 15 presents the LF of H α at z ∼ 0, 0.5, 1, and 2. Observational data points are collected from Ly et al. (2007), Gilbank et al. (2010), Drake et al. (2013), Sobral et al. (2013, 2015), Hayashi et al. (2018, 2020), and Khostovan et al. (2020). Our model effectively reproduces the observed H α LF across a wide luminosity spectrum, ranging from |$10^{39} $| to |$10^{43} \, {\rm erg\, s^{-1}}$|⁠, covering redshifts from 0 to 2. Since the luminosity of H α is directly related to the SFR by equation (5), the successful reproduction of the H α LF suggests that our model can accurately emulate the overall SFRF from redshift 0 to 2.

The second row in Fig. 15 displays the LF of [O ii] at z ∼ 0.5, 1, 1.5, and 2. The [O ii] luminosity is computed as the sum of the [O ii] doublet. Observed data points are collected from Ly et al. (2007), Zhu et al. (2009), Gilbank et al. (2010), Ciardullo et al. (2013), Drake et al. (2013), Comparat et al. (2015), Khostovan et al. (2015, 2020), Sobral et al. (2015), Hayashi et al. (2018, 2020), and Cedrés et al. (2021). The overall agreement between observations and simulation is satisfactory, although there is a slight overestimation of the LF at the knee at low redshift. It is worth noting that limitations arise from the survey volume, which only ranges from tens of thousands to several hundred thousand |$\rm Mpc^3$|⁠. Due to this constraint, the survey is susceptible to cosmic variance at these luminosities.

The third row in Fig. 15 illustrates the LF of the [O iii] doublet at z ∼ 0.5, 0.76, 1, and 1.5. Observational data points are collected from Ly et al. (2007), Drake et al. (2013), Hayashi et al. (2018, 2020), and Khostovan et al. (2020). The last row shows the luminosity function of [O iii]5007 + H β at z ∼ 0.75, 1.5, 2, and 3, with observational points collected from Khostovan et al. (2015, 2020) and Sobral et al. (2015). Our model predictions demonstrate consistency with the observed [O iii] doublet and [O iii]5007 + H β LFs within a wide range of luminosity and redshift.

In summary, our study successfully reproduces the observed luminosity function for a set of emission lines from the local universe to redshift 3. This achievement suggests that both our semi-analytic galaxy model and the emission line model align with the actual universe.

4.3 Clustering of ELG

We investigate the clustering of emission line galaxies by calculating their projected auto- and cross-correlation functions for various subsamples. Following the methodology outlined in Gao et al. (2022), we partition our simulated galaxies at redshift 0.76 into distinct subsamples based on their [O ii] luminosity and stellar mass. In practice, we create four samples based on [O ii] luminosity, denoted as L0, L1, L2, and L3. For the cross-correlation functions, we generate four samples based on their stellar masses, denoted as M0, M1, M2, and M3. Detailed information is provided in Table 4.

We present the projected auto and cross-correlation functions in Fig. 16. Each panel is for one particular [O ii] luminosity bin. Observed data points with error bars are obtained from Gao et al. (2022), derived from the VIMOS Public Extragalactic Redshift Survey (VIPERS; Scodeggio et al. 2018). In each panel, the black line represents the autocorrelation function of the specified [O ii] luminosity sample. Our model successfully replicates the autocorrelation functions of galaxies within different [O ii] luminosity bins across a broad range of radii, including one-halo and two-halo terms. The cross-correlation functions between the subsamples selected based on [O ii] luminosity and subsamples selected based on stellar mass are depicted by the cyan/yellow/green/purple lines. It should be noted that these cross-correlation lines have been appropriately shifted by a factor of 2ⁿ, where n corresponds to the specific designation of the sample, M0, M1, M2, M3. Results from different stellar mass selected samples are presented using different colours as denoted in the bottom left corner of each panel. Overall, our model predictions exhibit an excellent agreement with observations for all subsamples. The only exception appears for the projected cross-correlation functions between M3 and L3, where the model prediction is higher at small scales compared to observations, implying a stronger star formation close to massive central galaxies.

4.4 Bias of different emission line tracers

Emission line galaxies are one of the main targets for the next generation of large-scale surveys. It is essential to understand their bias relative to the matter distribution. Galaxy bias is often used in various cosmological probes, including baryonic acoustic oscillation, redshift distortion, etc. The bias is defined using:

where ξ_gg is the 3D correlation function of the galaxy sample, and ξ_mm is the correlation function of the total matter. A value of b = 1 signifies that the particular galaxy sample traces the distribution of the total matter, while when b is greater or less than 1, it means that the particular galaxy sample is more or less strongly clustered than the total mass.

Here, we calculate the bias of emission line galaxies as a function of luminosity and redshift. We split the simulated ELGs into bins of 0.2 dex within the luminosity range from 10⁴¹ to |$10^{44} \, {\rm erg\, s^{-1}}$|⁠. We focus on four of the most luminous lines, H α, H β, [O ii], and [O iii]. Bins containing fewer than 1000 galaxies are excluded. We compute the average bias over the range 20 to 60 |${\rm Mpc}\, h^{-1}$| where the bias is relatively constant and denote it as |$\overline{b}_{\rm gm}$|⁠. Fig. 17 presents |$\overline{b}_{\rm gm}$| as a function of luminosity at different redshift 0 to 3 for H α, H β, [O ii], and [O iii] selected galaxies. The errors are Poisson errors.

Fig. 17 shows that for Hα selected galaxies at z < 1 the bias exhibits a flat slope with luminosity for less luminous objects (⁠|$\lt 10^{42.5}\, {\rm erg\, s^{-1}}$|⁠) and undergoes a sharp increase at high luminosity. The luminosity transition increases as redshifts increase and disappears at z > 1, where the bias shows a more gradual increase with luminosity across the entire range of luminosity. At low luminosities, the bias monotonically increases with redshifts. However, at high luminosities, the bias decreases with redshifts up to z ∼ 1 and then increases with redshift above z ∼ 1. This is due to the sharp increase in luminosity at low redshifts which disappears at high redshifts.

Similar features are observed in galaxies selected based on other emission lines, although the exact redshift and luminosity dependence differ among different line selections. For instance, the transition luminosity is considerably smaller in galaxies selected based on [O iii] emission line compared to H α selected galaxies.

The variation in line selections can be understood as different line emissions corresponding to different physical conditions. Lines with the same luminosity may be hosted by galaxies with different properties. For example, for H |$\alpha \sim 10^{43} \, {\rm erg\, s^{-1}}$|⁠, the typical galaxy stellar mass is |$10^{11.43}\, {\rm M_\odot }$|⁠, SFR is |$51.26\, {\rm M_\odot }{\rm yr^{-1}}$|⁠, and halo mass is |$10^{13.49}\, {\rm M_\odot }$|⁠. In contrast, for [O iii] with the same luminosity, the typical galaxy stellar mass is |$10^{11.02}\, {\rm M_\odot }$|⁠, SFR is |$65.81\, {\rm M_\odot }{\rm yr^{-1}}$|⁠, and halo mass is |$10^{13.06}\, {\rm M_\odot }$|⁠.

Detailed comparison among different emission lines is presented in Fig. 18. At redshift 0, it is evident that for H α and [O ii] selected galaxies, there is almost no dependence on emission line luminosity up to |$10^{42.5} \, {\rm erg\, s^{-1}}$|⁠. The typical halo mass is |$10^{12.6}\, {\rm M_\odot }$| for those with luminosity around the transition. However, at higher luminosities, the bias increases rapidly, reaching ∼4 at |$10^{43} \, {\rm erg\, s^{-1}}$|⁠, corresponding to a typical halo mass of |$10^{13.5}\, {\rm M_\odot }$|⁠.

At z = 0, [O ii] selected galaxies always have a similar bias compared to H α selected galaxies. H β selected galaxies exhibit a comparable bias to H α selected and [O ii] selected galaxies with luminosities below |$10^{42} \, {\rm erg\, s^{-1}}$|⁠. The corresponding halo mass at the transition point is |$10^{12.55}\, {\rm M_\odot }$|⁠, which is also similar to the halo mass of the H α selected and [O ii] selected galaxies. Beyond this luminosity threshold, the H β selected galaxies experience a rapid increase in bias and surpass the bias of the H α selected and [O ii] selected galaxies, reaching a bias of approximately 2 at |$10^{42.5} \, {\rm erg\, s^{-1}}$|⁠. The transition luminosity happens at even smaller luminosity for [O iii] selected galaxies. Surprisingly, we find that for galaxies with the least luminous [O iii] emissions, the bias is well below 1, also much lower than H α, H β and [O ii] selected galaxies of the same luminosity. The corresponding typical stellar mass is |$10^{10.24}\, {\rm M_\odot }$|⁠, SFR is |$3.24\, {\rm M_\odot }{\rm yr^{-1}}$|⁠, and halo mass is |$10^{11.88}\, {\rm M_\odot }$|⁠. Such relations are also observed at other redshifts.

The relative bias of [O iii] selected galaxies compared to other line-selected galaxies varies with luminosities and redshifts. At high luminosities, H α and [O ii] selected galaxies have a lower bias than [O iii] selected galaxies at all redshift of interest. Conversely, at low luminosities, the [O iii] selected galaxies have a lower bias below z = 1. The transition luminosity increases with increasing redshifts. At higher redshifts, [O ii] selected galaxies have a similar bias compared to [O iii] selected galaxies.

For less luminous ELGs, especially at low redshifts, their bias dependence on luminosity is weak, posing challenges in modelling their spatial distribution using halo occupation models and abundance matching techniques. To quantify the relationship between |$\overline{b}_{\rm gm}$| and luminosity at different redshift, we employ a combined power law and exponential function to fit |$\overline{b}_{\rm gm}(L)$|⁠:

where Φ_L, L_*, α, and β are free parameters. We set a minimum value of log(L_*) = 40 since we do not use data from lower luminosity bins.

The fitting results are presented in Figs 17 and 18 by solid lines. Our fitting formula consistently reproduces the observed bias correlation at all luminosities across all redshifts. The only exception appears at the highest luminosity end for some redshifts, where the sample size is relatively too small. The fitting parameters as a function of redshift are listed in Table 5.

In summary, in combination with the SAM and emission line models, we successfully reproduce most of the observed properties of emission line galaxies properties, including their luminosity functions, correlation functions, and their evolution with redshifts. At high redshifts, their bias increases with both redshift and luminosity, while at low redshift, the bias is a decreasing function of redshift, attributed to the stronger dependence on luminosity towards lower redshifts. Across all redshifts, the luminosity dependence of galaxy bias is weak below |$10^{42} \, {\rm erg\, s^{-1}}$|

5 CONCLUSIONS

Our investigation focuses on generating a simulated galaxy catalogue for next-generation surveys, especially to include emission line galaxies in a self-consistent way. We solve the time convergence issue in the widely used semi-analytic model l-galaxies and employ it in the JiuTian-1G simulation. Furthermore, we compute the luminosity of various emission lines, enabling predictions about emission line galaxies. We further study the clustering and bias of different emission lines and provide a fitting formula for bias as a function of luminosity and redshift. Our catalogue successfully reproduces various observational properties. The main conclusions of our study can be summarized as follows:

• We observe a significant convergence problem in the l-galaxies model presented by Henriques et al. (2015) when applied to dark matter merger trees with varying time intervals. This issue stems from the disruption model being exclusively implemented at the end of each snapshot gap. Therefore, merger trees with fewer snapshots tend to experience reduced disruption and an increased number of mergers. As mergers predominantly contribute to the growth of SMBH, a higher frequency of mergers results in more massive SMBHs, leading to more efficient AGN feedback and, consequently, less massive galaxies. By modifying the disruption model in l-galaxies, We successfully achieve excellent convergence in simulations with merger trees of varying time intervals.

• Our adapted model is applied to the JiuTian-1G simulation, and the corresponding parameters were readjusted. Our catalogue successfully reproduces numerous statistical observational properties and accurately captures the clustering patterns of diverse galaxy samples. In particular, it has been able to replicate the SMBH mass function, which was underestimated by H15 by an order of magnitude.

• We demonstrate that in combination with the high resolution large boxsize JT1G simulation, l-galaxies has successfully generated a galaxy catalogue that fulfils most of the requirements of next-generation large-scale surveys.

• We compute the luminosity of 13 commonly used NUV and optical emission lines. Our model effectively reproduces the observed luminosity function of H α, H β, [O ii], and [O iii]. Additionally, the projected correlation of [O ii] ELGs shows good agreement with observations.

• We further explore the bias of emission line galaxies as a function of luminosity and redshift. The dependence varies with luminosity ranges and redshifts. We observe that at low redshift, the bias of galaxies with low luminosity shows insensitivity to luminosity, while it increases rapidly at the high luminosity end. At high redshift, bias gradually increases with luminosity. Above z = 1, galaxy bias increases monotonically with redshift, while below z = 1, such a monotonic increase only holds for low luminosity galaxies. At high luminosities, galaxy bias decreases with redshift up to z = 1 and then increases with redshift.

• We offer fitting formulas that capture the dependence of bias on both luminosity and redshifts.

In conclusion, our adapted model successfully replicates various galaxy observational properties, and the predictions from our emission line model align well with observations. The bias has a complex dependence on luminosity and redshift, which varies with luminosity range and redshift range.

Due to the limitation of storage, we only present photometric magnitude from several surveys and the luminosity of 13 emission lines. Additional photometric magnitude, emission lines, and full SED are available upon request.

ACKNOWLEDGEMENTS

We thank Professor Simon White for his suggestions. This work is supported by the National SKA Program of China (Nos. 2022SKA0110201 and 2022SKA0110200), and CAS Project for Young Scientists in Basic Research grant no. YSBR-062, the National Natural Science Foundation of China (NSFC) (grant nos. 12033008, 12273053, 11988101, 11622325), the K.C. Wong Education Foundation, the science research grants from China Manned Space Project with No.CMS-CSST-2021-A03 and No.CMS-CSST-2021-B03, and the Strategic Priority Research Program of Chinese Academy of Sciences, Grant no. XDB0500203. QG thanks European Union’s HORIZON-MSCA-2021-SE-01 Research and Innovation programme under the Marie Sklodowska-Curie grant agreement number 101086388. JH acknowledges support from the Yangyang development fund.

DATA AVAILABILITY

The data produced in this paper are available upon reasonable request to the corresponding author.

Footnotes

References

APPENDIX A: OBSERVATIONAL CONSTRAINTS USED IN MCMC

In the MCMC procedure, we use the stellar mass function at z = 0, 1, 2, and 3, the passive fraction at z = 0.4, and the black hole mass function at z = 0 as constraints to find the best-fitting parameters. The observational data we use in this study are listed in Table A1. The constraints of the SMF at z = 0, 1, 2, and 3 are combined from the listed data set (see Appendix A in Henriques et al. 2015, for details), while the passive fraction and BHMF only rely on data from a single observation.

APPENDIX B: CONVERGENCE OF OUR NEW MODEL

In order to assess the convergence of our model across simulations with varying dark matter particle masses, we conducted our analyses on the high-precision rescaled Millennium-II (MSII) simulation (Boylan-Kolchin et al. 2009). The dark matter particle mass in the MSII simulation is set to |$7.69\times 10^{6}\, {\rm M_\odot }\, h^{-1}$|⁠, with precision nearly 50 times greater than that of the JT1G simulation. Furthermore, the MSII simulation has been scaled to the cosmological parameters of Planck14, which differ slightly from those used in our Planck18-based model. Therefore, comparing the performance of our model across the MSII and JT1G simulations allows for a more nuanced understanding of its convergence concerning varying resolutions and cosmological parameters. Fig. B1 illustrates the performance of our model in predicting the SMF from redshift 0 to 3 on both the JT1G and MSII simulations. The red curves represent the JT1G simulation, while the black curves represent the MSII simulation. The upper row depicts the SMF, and the lower row illustrates the ratio of JT1G to MSII. Overall, our model exhibits a robust agreement between JT1G and MSII within the range of |$10^{9}$| to |$10^{10.5}\, {\rm M_\odot }$|⁠, with differences consistently below 10 per cent. At the lower mass end, there is a noticeable overestimation due to the precision limitations of the JT1G simulation. Conversely, at the higher mass end, the MSII simulation, constrained by its box volume, lacks massive galaxies.

Fig. B2 is the same as Fig. B1, but presenting the Star Formation Rate Function (SFRF) across redshifts 0 to 3. At lower redshifts, the agreement between JT1G and MSII results is excellent, with a deviation of less than 5 per cent. However, at higher redshifts, JT1G deviates significantly from MSII at lower SFR values, indicating a limitation in resolution. Notably, at SFR values greater than 1 |${\rm M_\odot }\,{\rm yr^{-1}}$|⁠, JT1G and MSII results remain in close agreement.

In summary, our model demonstrates outstanding temporal resolution convergence and satisfactory convergence across simulations with different particle masses. The JT1G simulation provides comprehensive and resolution-independent coverage for galaxies with stellar masses above |$10^{9}\, {\rm M_\odot }$| and SFR exceeding 1 |${\rm M_\odot }\,{\rm yr^{-1}}$|⁠, up to redshift 3.

APPENDIX C: SEDS OF TWO TYPICAL GALAXIES

In general, we can calculate the galaxy SED with emission lines for each galaxy in the catalogue. Here, the SEDs of a typical star-forming galaxy and a passive galaxy are presented in Fig. C1 using the method described in Secion 2.3. We include the 13 emission lines in the SED for star-forming galaxies. It illustrates the declining feature with increasing wavelength for the star-forming galaxy, the D4000 break and UV-upturn for the passive galaxy.