Predictions of the L_{[C ii]}–SFR and [Cii] Luminosity Function at the Epoch of Reionization

Deep field surveys carried out with the Hubble Space Telescope have enabled the detection of galaxies out to the cosmic dawn at z ≳  11 and provide stringent constraints on the bright end of rest-frame UV luminosity functions (LF) of galaxies during the epoch of reionization (EOR; e.g., Finkelstein 2016; Oesch et al. 2016; Song et al. 2016). While measuring distributions of galaxy properties, such as the LF, provides important constraints on how galaxies evolve over cosmic time, it is also useful to target individual sources at high redshift in order to conduct detailed case studies. Such studies allow us to address open questions such as properties of the interstellar medium (ISM), how metal enrichment in galaxies proceeded in the early universe, how much star formation is dust obscured or how much dust is present in these early systems, and the physical conditions under which star formation and stellar mass assembly took place. Follow-up observations with the Atacama Large Millimeter/submillimeter Array (ALMA) probing the dust continuum and line emission from ions and molecules in the dense ISM of high-redshift galaxies appear extremely promising for characterizing the star-forming gas and dust properties of early galaxies and helping to constrain the physical processes at play (see a review by Hodge & da Cunha 2020).

Owing to the relative brightness of the fine-structure line from singly ionized carbon at rest frame 157.7 μm and accessibility of the line with ALMA at high redshift, there is currently great interest in using [Cii] line emission as a tracer to study galaxies during the EOR. At the moment, most of these follow-up observations target some of the brightest objects selected at other wavelengths (e.g., Capak et al. 2015; Marrone et al. 2018; Smit et al. 2018). [Cii] has been detected in a handful of normal star-forming galaxies at the EOR (e.g., Capak et al. 2015; Carniani et al. 2018) and even resolved on kiloparsec scales to study the gas kinematics for a handful of sources (e.g., Jones et al. 2017; Matthee et al. 2017; Carniani et al. 2018; Hashimoto et al. 2019). However, due to the small field of view of ALMA, it is extremely challenging and expensive to carry out flux-limited (untargeted) surveys over significant areas. The largest-area untargeted survey to date that probes [Cii] at z ∼ 6−8 is the ALMA Spectroscopic Survey in the HUDF (Large Program; ASPECS), which covers 4.6 arcmin² (Aravena et al. 2016; Walter et al. 2016).

While these detections are exciting, interpretation of the [Cii] line luminosity and its connection with galaxy properties remains complex and poorly understood. [Cii] is the dominant coolant in the ISM in nearby star-forming galaxies, and its luminosity is expected to be correlated with the star formation rate (SFR). Indeed, such a correlation is seen in local galaxies (De Looze et al. 2014; Herrera-Camus et al. 2015), but with hints that the L_{[C ii]}/SFR ratio depends on other galaxy properties such as metallicity and dust temperature (e.g., Malhotra et al. 2001; Luhman et al. 2003; Díaz-Santos et al. 2014). This is expected, as the dust content affects the degree of shielding against hard ionizing photons, as well as the amount of photoelectric heating, which affect the ionization state of the line-emitting gas and the collision rate. In addition, theoretical modeling has shown that L_{[C ii]}/SFR is also expected to depend on the pressure of the ISM environment, if molecular cloud sizes depend on the ambient pressure, and on the density distribution on small scales within the ISM (Popping et al. 2019).

Another motivation for studying the physics of line emission and its connection to galaxy properties and the physics of galaxy formation is the development of multiple experiments that will carry out line intensity mapping (LIM) studies. Different LIM experiments will enable the detection of various tracers of dense and more diffuse gas in and around galaxies, including Lyα, Hα, 21 cm, CO, and [Cii] (see Kovetz et al. 2017 for a review). LIM promises to map the statistical signal from emission in galaxies over very large volumes (on the order of megaparsec to gigaparsec scales), with the trade-off of not resolving individual galaxies. As such, LIM experiments also probe emission from faint galaxies. While LIM holds great promise for studying galaxy evolution and cosmology (e.g., the cosmic star formation history, evolution of the ISM and intergalactic medium [IGM], and physical processes during the EOR; Kovetz et al. 2017), the power spectra depend on the line luminosities of different galaxy populations within the volume sampled (or the intrinsic line LF). Therefore, physically motivated models that can self-consistently predict line luminosities are vital for strategizing LIM experiments and interpreting LIM data. As shown by Yue & Ferrara (2019), a shallower L_{[C ii]}–SFR relation would imply that the [Cii] LF drops quickly at the bright end and most of the IM signal comes from faint galaxies. This in turn determines the detection depth needed for [Cii] LIM experiments. At the moment, most LIM forecasts for tracers such as CO and [Cii] have been made using a series of empirical scaling relations (e.g., Gong et al. 2012; Uzgil et al. 2014; Keating et al. 2015; Chung et al. 2019).

Carrying out detailed and realistic predictions of the [Cii] emission for a large cosmologically representative sample of galaxies is extremely challenging. The [Cii] line can arise from all phases of the ISM, by being collisionally excited by electrons, atoms, or molecules. Hence, the line strength depends strongly on the density and kinematic temperature of these species. Modeling has repeatedly shown that different ISM phases are all important to consider when deriving [Cii] emission for a galaxy (e.g., Olsen et al. 2017; Pallottini et al. 2019; Lupi & Bovino 2020). However, state-of-the-art cosmological simulations of volumes larger than 100 Mpc do not resolve particle masses below ∼10⁶−10⁷ M_⊙ (at z = 0; e.g., Davé et al. 2019; Nelson et al. 2018), which corresponds to hydrogen densities below n_H < 100 cm⁻³ and temperatures above 10⁴ K. In particular, the molecular ISM phase—with typical densities above a few hundred cm³ and temperatures below 100 K—is not tracked in cosmological simulations, but knowledge about it is critical in order to reliably simulate [Cii] line emission. Currently, all large-volume cosmological simulations adopt phenomenological "subgrid" recipes to treat unresolved processes such as star formation, stellar feedback, and chemical enrichment, and these carry significant uncertainties (see review by Somerville and Davé 2015). Additionally, a second type of "subgrid" modeling is required to describe the detailed structures of molecular gas on cloud scales, which must be input into the radiative transfer (RT) and ionization state modeling tools, described in more detail below.

Previous works that have attempted to model [Cii] for galaxy populations mainly fall into three categories: (i) extremely simple, empirical mappings between L_{[C ii]} and halo mass (Visbal & Loeb 2010; Gong et al. 2012), (ii) semianalytic (SAM) models to predict galaxy properties coupled with machinery to predict the [Cii] emission in post-processing (Popping et al. 2016; Lagache et al. 2018; Popping et al. 2019), or (iii) small samples of cosmological zoom-in simulations again post-processed to compute the RT and line emission (Olsen et al. 2015; Narayanan et al. 2015; Vallini et al. 2015; Katz et al. 2017; Olsen et al. 2017; Pallottini et al. 2017b, 2017a). A few recent studies have coupled on-the-fly RT, nonequilibrium chemistry modeling, and line spectral synthesis with ultra−high-resolution hydrodynamic zoom-in simulations (Katz et al. 2017, 2019; Pallottini et al. 2019).

A variety of tools are brought to bear to compute the emergent line emission in the literature. A stellar population synthesis code such as Starburst99 is used to derive the amount of intrinsic radiation from stellar emission (Leitherer et al. 2014). Codes such as Skirt or Powderday are used to perform dust RT (Camps & Baes 2015; Narayanan et al. 2015), and tools such as RADMC-3D, Despotic, Lime, Cloudy, or MAPPINGS are used to compute line RT in the ISM (Allen et al. 2008; Krumholz 2014; Ferland et al. 2017). For a detailed summary and comparison of these codes, we refer interested readers to Olsen et al. (2018).

Some of the biggest differences among the aforementioned codes are the density range considered/permitted, the geometry, and the species included in the chemical network. These differences also imply different demands on computational time and memory, each with different benefits and trade-offs. The accuracy needed for a given galaxy simulation and the emission line of interest typically determine the method used. For instance, a photodissociation region (PDR) code such as Despotic can be used to calculate line emission from the neutral ISM; however, Despotic does not simulate line emission originating in the ionized phase of the ISM.

In addition to different methods being adopted for RT and line spectral synthesis, the type of simulation used also sets limits on the galaxy sample size obtainable and the level of realistic physics that can be adopted. SAMs are computationally inexpensive and can easily generate galaxy catalogs of statistically significant sample sizes (e.g., Popping et al. 2016; Lagache et al. 2018; Popping et al. 2019) and thus are excellent for testing physical recipes and exploring wide ranges of parameter space, but they do not provide any detailed information on subgalactic structure. In contrast, the highest-resolution zoom-in hydrodynamical simulations can numerically resolve the ISM down to scales of ∼10 pc at high redshifts (e.g., Katz et al. 2017; Pallottini et al. 2017a, 2019) but are computationally expensive. The samples of galaxies with simulated line emission based on numerical hydrodynamical simulations are thus limited in number to 1–30 per study, and therefore they also probe a limited parameter space of galaxy properties. Previous studies commonly targeted the most massive halos and/or a handful of sources with properties resembling some known properties of a given observed z ∼ 6 galaxy. Over the years, the resolution of large-volume cosmological hydrodynamics simulations has increased significantly, reaching ≳100 pc resolution even before carrying out additional refinement using zoom-in techniques. This produces a statistically significant and unbiased sample of galaxies while reaching down to a relevant spatial scale to simulate emission emerging from the ISM (though subgrid models are still required). Inoue et al. (2020) gave a successful demonstration of calculating CO line emission in direct post-processing of the cosmological IllustrisTNG simulation using a simulation box size of 75 Mpc.

In this work, we leverage the new simba suite of cosmological hydrodynamic simulations (Davé et al. 2019) to select a galaxy sample at z ∼ 6 for line emission post-processing that is unprecedented in its size (11,137 galaxies) and dynamic range (halo mass  ∼10⁹−10¹² M_⊙;  stellar mass  ∼10⁷−10¹¹ M_⊙). The sample is drawn from a set of volumes that vary in resolution, such that the sample from each volume is representative of the population that is well resolved in the simulation. We apply an updated version of the sígame package (Olsen et al. 2017), which includes subresolution modeling of the ISM, RT, and line spectral synthesis. We use this calculation to provide predictions of how the [Cii] luminosity at this redshift is related to other galaxy and DM halo properties and of the [Cii] LF, and we compare our results with available observations and with other models in the literature.

The paper is structured as follows. In Section 2 we describe the simba simulations, and in Section 3 we describe the method used to simulate [Cii] line emission. We present the results in Section 4 and discuss the limitations in Section 5. Finally, conclusions and implications of our findings are presented in Section 6. Throughout this paper, we adopt a concordance cosmology, with total matter, vacuum, and baryonic densities in units of the critical density Ω_Λ = 0.693, Ω_m  = 0.307, Ω_b  = 0.048, a dimensionless Hubble parameter h = 0.678, scalar spectral index of n = 0.96, and power spectrum normalization of σ₈ = 0.823 (Planck Collaboration et al. 2016).

2.1. Cosmological Hydrodynamics Simulation: simba

We use galaxies from the simba cosmological galaxy formation simulations (Davé et al. 2019) for this study. simba is a suite of gizmo-based simulations using meshless finite mass hydrodynamics, incorporating state-of-the-art feedback modules that provide very good agreement with a wide range of lower-redshift observables. The suite consists of random cubical volumes of 100, 50, and 25 cMpc h⁻¹ on a side, each with 1024³ dark matter particles and 1024³ gas elements. By combining the results of these simulations, we achieve unprecedented dynamic range, with the highest resolution equaling, e.g., that in a study of far-infrared lines using zoom simulations presented in Olsen et al. (2017). The smoothing lengths and the initial gas mass resolutions for the 100, 50, and 25 cMpc h⁻¹ volumes are _min = 0.5, 0.25, and 0.125 h⁻¹ kpc and m_gas = 1.8 × 10⁷ M_⊙, 2.3 × 10⁶ M_⊙, and 2.9 × 10⁵ M_⊙, respectively (see Table 1 of Davé et al. 2019). These runs all use identical input physics, begin at z = 249, and assume a Planck-concordant cosmology.

simba is the successor of the mufasa simulations (Davé et al. 2016), and details of the improvements in simba are provided in Davé et al. (2019). Among the various updates, key ones relevant to this work are as follows: (i) simba uses the grackle-3.1 package to model radiative cooling and photoionization heating, updated from mufasa to apply radiative processes via isochoric substep cycling, and also computing the neutral hydrogen content accounting for self-shielding on the fly via the prescription in Rahmati et al. (2013); (ii) simba explicitly models the growth and feedback of supermassive black holes (SMBHs) residing in galaxies, with the growth of the SMBH set by torque-limited accretion of cold gas (Anglés-Alcázar et al. 2017a) and Bondi accretion of hot gas, while black hole feedback is modeled via bipolar kinetic outflows and injection of X-ray energy; (iii) simba includes a subgrid model to form and destroy dust within the ISM of galaxies during the simulation run, with the dust mainly produced by Type II supernovae (SNe), asymptotic giant branch (AGB) stars, and condensation from metals and destroyed predominantly via sputtering (including SN shocks) and consumption by star formation (Li et al. 2019); (iv) simba employs ejective star formation feedback like mufasa, but with scalings updated to reflect particle tracking results from the Feedback in Realistic Environment (FIRE) simulations (Anglés-Alcázar et al. 2017b), with minor modifications to better reproduce EOR galaxy properties (see Wu et al. 2020, for details).

simba has been compared to a wide range of observations across cosmic time and is found to provide reasonable agreement, including for the galaxy stellar mass function (GSMF) and mass–metallicity relation (Davé et al. 2019), black hole properties (Thomas et al. 2019), dust properties (Li et al. 2019), galaxy sizes and profiles (Appleby et al. 2020), and cold gas contents including CO (J = 1  →  0) LFs from z = 0 to 2 (Davé et al. 2020). Minor disagreements with observations include an overproduction of the very highest mass galaxies at z ≲ 2, too large size for low-mass quenched galaxies at z ≲ 1, and an underproduction of the dust mass function at z ∼ 2. Relevant to this work where we focus on z = 6, Wu et al. (2020) examined the EOR properties of simba galaxies and found good agreement at z = 6 with the UV LF, UV slope measurements when a Calzetti (2001) dust model is assumed, and galaxy sizes down to the faintest limits. Hence, simba provides a realistic platform to examine the far-infrared line properties of EOR galaxies, which is the goal of this study.

In Figure 1, we show the GSMF and SFR function at z ∼ 6. We point out the robust numerical convergence in the GSMF and SFR function (and galaxy properties; see Davé et al. 2019 and Section 2.2) without refining or fine-tuning the parameters in the subgrid models of simba between the various simulation boxes at the redshift studied in this work. Note that we do not impose any stellar mass limits in the GSMF plotted. This illustrates the exceptional convergence reached across the boxes. This is crucial, as we make predictions of the [Cii] LF in the luminosity range of 5.5 < log (L_{[C ii]}/L_⊙) < 8.5 by combining galaxies in the 25 cMpc h⁻¹, 50 cMpc h⁻¹, and 100 cMpc h⁻¹ boxes (hereafter Simba-25, Simba-50, and Simba-100, respectively).

Zoom In Zoom Out Reset image size

2.2. Main Sample: 11,137 Galaxies at z ≃ 6

Galaxies from the simulation suite are identified using a galaxy finder that adopts a six-dimensional friends-of-friends algorithm (caesar). For the purpose of this work, we include only galaxies that have at least 64 stellar and gas particles, respectively, to ensure that they are resolved in the simulation. As illustrated in Figure 1, these mass requirements correspond to $\mathrm{log}({M}_{\star ,\min }/{M}_{\odot })=7.24$ for Simba-25, $\mathrm{log}({M}_{\star ,\min }/{M}_{\odot })=8.15$ for Simba-50, and $\mathrm{log}({M}_{\star ,\min }/{M}_{\odot })=9.05$ for Simba-100. Similarly, we impose thresholds on the SFR averaged over 10 Myr based on the turnover seen in the SFR function indicating incompleteness in SFR. This corresponds to $\mathrm{log}\left({\rm{SFR}}/{{\rm{M}}}_{\odot }\,{{\rm{yr}}}^{-1}\right)\gt -1.9$, −0.8, and 0.4 for Simba-25, Simba-50, and Simba-100, respectively. In addition, we only include galaxies with a molecular gas mass of at least M_mol = f_H2 M_gas > 10⁵ M_⊙, as the subgrid model has to form giant molecular clouds (GMCs) of at least 10⁴ M_⊙ each by sampling the cloud mass from a GMC mass function (see Section 3). After applying these criteria, we have a sample of N_tot = 11,137 galaxies for a single snapshot at z = 6. The ranges of their physical properties are listed in Table 1.

Table 1. Parameter Space Probed by Our simba Galaxy Sample at z ∼ 6

Properties	Units	Ranges
$\mathrm{log}{M}_{\mathrm{halo}}$	M⊙	[9.02 12.36]
SFR10	M⊙ yr−1	[0.01, 708]
SFR100	M⊙ yr−1	[0.007, 329]
$\mathrm{log}{{\rm{\Sigma }}}_{\mathrm{SFR}}$	M⊙ yr−1 kpc−2	[−2.60, 1.31]
log M*	M⊙	[7.18, 10.72]
log Mgas	M⊙	[7.27, 10.46]
〈Zgas〉SFR	Z⊙	[0.06, 1.86]

Download table as:  ASCIITypeset image

To illustrate the range of physical properties sampled by the simba galaxies, Figure 2 shows the relations between the specific SFR (sSFR), SFR-weighted metallicity (〈Z_gas〉_SFR), molecular gas-to-stellar mass ratio (M_mol/M_*), and the stellar mass-weighted age in our z = 6 sample. The three clumps of points represent galaxies from our three simulation volumes and generally show reasonable convergence. The weighted quantities are indicated with the 〈...〉 notation (e.g., 〈Z_gas〉_SFR), as defined as follows:

$$\begin{eqnarray}&&\langle x\rangle \equiv \displaystyle \frac{{\displaystyle \sum }_{i}{\rho }_{i}{x}_{i}}{{\displaystyle \sum }_{i}{\rho }_{i}},\end{eqnarray} \tag{ 1 }$$

where x is the variable and ρ_i is the volume of each fluid element i. The Σ_SFR in this paper is defined as

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{SFR}}=\displaystyle \frac{{\rm{SFR}}}{\pi {R}_{1/2}^{2}},\end{eqnarray} \tag{ 2 }$$

where R_1/2 is the half-mass radius of gas of a given galaxy. The scaling relations shown are similar to the mass–metallicity relation and fundamental metallicity relation (a.k.a. MZR and FMR; e.g., Maiolino et al. 2008; Mannucci et al. 2010), which are commonly used to gain insights into the interplay between star formation, gas accretion, and feedback during the evolution of a galaxy and are shown to illustrate the range of physical properties sampled by the simba galaxies. As can be seen in the second panel, most of the simba galaxies are rich in molecular gas. The trend of decreasing M_mol/M_* with increasing M_* arises owing to stellar feedback, which preferentially suppresses the stellar mass in lower-mass systems, thereby increasing M_mol/M_*. The middle panel also shows how, for a given stellar mass bin, the SFR increases with the molecular gas mass fraction, as expected. That said, there are certainly jumps in M_mol/M_* fractions between the different simulation volumes indicating less-than-ideal convergence, although this becomes less apparent when plotting M_mol against M_* as shown in the left panel of Figure 12 in Appendix A. The bottom panel displays the simba galaxies on top of the "star-forming main sequence" (SFMS; e.g., Speagle et al. 2014; Iyer et al. 2018), which most of our galaxies follow at $\mathrm{log}{M}_{* }\,\gtrsim 9$ M_⊙—which is also the stellar mass limit of the observational data. At lower stellar mass, the sSFR of simba galaxies falls below the SFMS extrapolation. ¹¹ The color-coding in the three panels shows that, at a given stellar mass, galaxies that are higher metallicity have lower gas contents, galaxies that have higher gas contents have higher sSFR, and galaxies that have higher sSFR have lower mean stellar age.

Zoom In Zoom Out Reset image size

Galaxies of similar stellar masses and SFRs may have different sizes, surface densities, gas contents, metallicities, interstellar radiation field strengths, structural properties, and gas dynamics, all of which would produce varying [Cii] luminosities (see, e.g., Kaufman et al. 1999; Vallini et al. 2015; Olsen et al. 2017). As such, comparing the physical properties of observed (i.e., [Cii]-detected ones in the context of this work) and simulated galaxies is pertinent to establishing the reliability of model predictions. In other words, comparing these global properties of simba galaxies with those of observed galaxies enables one to place the observed ones, given their [Cii] luminosities, in a theoretical framework. It is beyond the scope of this paper to perform a detailed comparative study since the information available for observed galaxies at these redshifts (see Section 4.1) remains limited and inhomogeneous. For comparisons of simba galaxies with observations at lower redshifts and discussion on redshift evolution in these relations, we refer interested readers to Davé et al. (2019). At the EOR, the size–luminosity relation of simba galaxies agrees with observations (Kawamata et al. 2018; Wu et al. 2020).

We use an updated version of sígame (Olsen et al. 2015, 2017) to post-process the simba simulation outputs. For details of the code, we refer interested readers to Olsen et al. (2017). Here we briefly summarize the salient points of sígame and updates made to the code as part of this work. For each gas fluid element, sígame divides the molecular gas mass (i.e., f_H2, i m_gas, i) into GMCs by sampling the Galactic GMC mass function over the mass range of 10⁴–10⁶ M_⊙ (${dn}/{dM}\propto {M}_{\mathrm{GMC}}^{-1.8};$ see, e.g., the review by McKee & Ostriker 2007; Blitz et al. 2007; Kennicutt & Evans 2012). The remaining mass of the parent fluid element is assumed to be in the diffuse gas phase and is subsequently distributed into diffuse ionized and neutral gas phases. This division is determined by the boundary at which the inner neutral region transitions to the outer ionized region of the diffuse clouds as computed from RT calculations within Cloudy. That is, the neutral gas phase corresponds to the region beyond a radius where the neutral fraction x_HI = n_HI/(n_HI + n_HII) >  0.5, such that it is dominated by neutral hydrogen—here n is the number density. As such, sígame accounts for line emission from three distinct ISM phases. The smoothing length of the parent fluid element is adopted as the size of the diffuse gas clouds, whereas the size of each GMC is derived from a pressure-normalized mass–size relation, following

$$\begin{eqnarray}&&\displaystyle \frac{{R}_{\mathrm{GMC}}}{\mathrm{pc}}={\left(\displaystyle \frac{{P}_{\mathrm{ext}}/{k}_{{\rm{B}}}}{{10}^{4}{\mathrm{cm}}^{-3}{\rm{K}}}\right)}^{-1/4}{\left(\displaystyle \frac{{M}_{\mathrm{GMC}}}{290{M}_{\odot }}\right)}^{1/2},\end{eqnarray} \tag{ 3 }$$

where k_B is the Boltzmann constant and the external cloud pressure P_ext is defined assuming midplane hydrostatic equilibrium within the galaxy.

Each GMC radial density profile is assumed to follow a truncated logotropic profile. Both GMCs and diffuse gas phases inherit the metallicity of their parent fluid element. The far-UV (FUV) luminosity of each star is calculated based on its age and metallicity and is determined by interpolating over a grid of starburst99 stellar population synthesis models (Leitherer et al. 2014; with the default Kroupa 2002 initial mass function). Each GMC is irradiated by a local FUV radiation (6–13.6 eV), where the strength of the radiation field (G₀ in Habing units; 1.6  × 10⁻³ erg cm⁻² s ⁻¹) is determined by summing up the FUV flux from nearby stellar particles and by assuming that the flux falls off as 1/r². In the diffuse gas phase, the FUV radiation field is determined based on the SFR surface density of the galaxy. For our sample, the SFR surface density ranges between ≃1 and 6200 times the Milky Way.

We use the photoionization code cloudy version 17.01 (Ferland et al. 2017) to simulate the thermochemistry in the three distinct ISM phases tracked by sígame by performing detailed balance calculations of the various species, taking into account physical processes such as H₂ photo processes, dust physics (grain-atom/ion charge transfer), and cosmic-ray (CR) ionization. The line luminosities are then derived from the cooling rates for different line transitions. For computational purposes, lookup tables are generated for the GMC and the diffuse (neutral and ionized) gas phase, respectively. The FUV radiation field impinging on the gas phases is assumed to have the same spectral shape as in the solar neighborhood.

CRs are added, with an ionization rate equal to that of the Milky Way scaled linearly by a factor of (G_0,gas/G_0,MW). For the GMC models, the clouds are in theory completely embedded within diffuse gas, and thus H-ionizing radiation is turned off in the Cloudy models (see Olsen et al. 2017).

The main parameters in the GMC phase of the cloudy models considered are the G_0, GMC of the radiation source, radius of the cloud (R_GMC), and cloud density profile as a function of cloud radius (n_H(R_GMC)). Turbulent velocity is added to the GMC models according to the velocity dispersion calculated from the cloud radius and pressure, assuming that clouds are virialized. For the diffuse gas phase, the main model parameters are gas density (n_H), gas kinetic temperature (T_k ), diffuse cloud size, (R_dif), and metallicity (Z).

Compared to Olsen et al. (2017), the main updates made to sígame used in this work are as follows:

• 1.  
  Instead of fixing the number and width of shells used by cloudy to model each GMC, we now allow cloudy to determine the optimal quantities to ensure convergence. This modification leads to a more accurate calculation of the grain photoelectric heating of the gas and increases the importance of gas heating due to this mechanism in the GMC models, which is the main excitation mechanism for [Cii] emission. Namely, [Cii] is collisionally excited such that higher kinetic temperature leads to more molecular motions and collisions inside GMCs and PDRs, the main sites for [Cii] emission in galaxies.
• 2.  
  To ensure good sampling of the parameter space for both GMCs and diffuse gas clouds, the number of cloudy models used to create lookup tables is significantly increased from 1296 to 4096 models by using eight grid points in each parameter space dimension rather than six as in Olsen et al. (2017). The lookup tables are further described in Section 3.1.
• 3.  
  The dust content of the ISM is a crucial factor in setting the [Cii] luminosity. As often done, we will assume here that dust scales with metallicity via the dust-to-metal (DTM) ratio, but instead of using a solar DTM value of ∼0.46 (as done in Olsen et al. 2017), we take a DTM of 0.25 based on the mean value of our simba galaxies (see Section 3.1).

3.1. GMCs and Diffuse Gas-phase Properties and cloudy Model Grids

3.2. Dust and Elemental Abundances

While simba tracks dust in the simulation, we do not create a different cloudy lookup table for each DTM ratio found in the simulated galaxies since this becomes computationally intractable (i.e., it corresponds to a hyperparameter where each DTM ratio would have a separate set of 4096 cloudy models). Instead, we adopt a DTM ratio based on the median of simba galaxies at z ∼ 6, corresponding to ξ_DTM = 0.25, which is defined as

$$\begin{eqnarray}&&{\xi }_{\mathrm{DTM}}=\displaystyle \frac{{M}_{\mathrm{dust}}}{{f}_{Z}\,{M}_{\mathrm{gas}}+{M}_{\mathrm{dust}}},\end{eqnarray} \tag{ 4 }$$

where M_dust and M_gas are the dust mass and total gas mass in solar mass units, respectively, and f_Z is the mass fraction of metals (i.e., f_Z M_gas yields the mass of metals in gas phase). The dust content of each cloud is then set to scale linearly with its metallicity through this DTM ratio. The default set of lookup tables in Cloudy assumes a DTM ratio of ξ_DTM = 0.46 at solar metallicity. A more commonly adopted expression for the DTM ratio is

$$\begin{eqnarray}&&{\rm{DTM}}\equiv \displaystyle \frac{{\rm{DGR}}}{{\rm{Z}}}.\end{eqnarray} \tag{ 5 }$$

In the Milky Way, Z = Z_⊙ and DGR is ∼0.01, yielding log DTM = −2.

In cloudy, one can supply the total metallicity, and the abundance of each of the elements is scaled correspondingly assuming solar composition (i.e., abundance ratios of the Sun). In order to account for abundance patterns of galaxies that differ from solar, we use the abundance of each metal tracked in simba (i.e., He, C, N, O, Ne, Mg, Si, S, Ca, and Fe). For elements that are not tracked in the simulation, we use solar abundance ratios. Since the mass fraction of each element tracked varies as a function of metallicity, we fit a spline curve to the running means across all gas fluid elements. This provides a function that maps a given metallicity to an abundance pattern, such that the relative elemental abundances in the cloudy input are scaled according to the metallicity of the cloud. In Appendix B we show the particle distribution and scalings found for our sample of simba galaxies. We note that some elements, like carbon and nitrogen, can be quite far from their solar abundance value, which in simba is a result of including enrichment from Type II SNe, Type Ia SNe, and AGB stars, with separate yield tables for each class of star as described in Oppenheimer and Davé (2006).

4.1. [Cii]—SFR Relation at z ≃ 6

In Figure 3, we plot the simulated L_{[C ii]} and SFR ¹² of the simba galaxies, together with measurements from existing observations at z ≃ 6, local measurements, and other model predictions at z ≃ 6. The L_{[C ii]}–SFR relation converges across the different simulation volumes as does the L_{[C ii]}–M_mol relation as seen in the right panel of Figure 12 in Appendix A. For clarity, information shown in this figure is also plotted across three panels in Figure 4. We fit a linear model to L_{[C ii]} and SFR in log–log space to facilitate comparison with literature work and obtain

$$\begin{eqnarray}&&\mathrm{log}{L}_{[{\rm{C}}{\rm\small{II}}]}=(6.82\pm 0.08)+(0.66\pm 0.01)\times \mathrm{logSFR},\end{eqnarray} \tag{ 6 }$$

where L_{[C ii]} is in units of L_⊙ and SFR is in units of M_⊙ yr⁻¹.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The [Cii] luminosities of the simba galaxies are consistent with existing upper limits and a handful of detections from targeted observations (e.g., Ouchi et al. 2013; Kanekar et al. 2013; Ota et al. 2014; González-López et al. 2014; Maiolino et al. 2015; Schaerer et al. 2015; Capak et al. 2015; Willott et al. 2015; Inoue et al. 2016; Pentericci et al. 2016; Knudsen et al. 2016; Bradač et al. 2017; Decarli et al. 2017; Knudsen et al. 2017; Smit et al. 2017; Carniani et al. 2018). In addition, our results are in agreement with the latest upper limits placed at z ≃ 6 by the ALMA large program ASPECS, which is an untargeted survey, placing an upper limit of L_{[C ii]} < 2  × 10⁸ L_⊙ for galaxies with UV-derived SFR of ∼ 0.25–50 M_⊙ yr⁻¹ (Walter et al. 2016, Uzgil et al. 2020, in preparation). ¹³ In particular, the ASPECS sources with upper limits on L_{[C ii]} shown in Figure 3 with blue squares are a combination of Lyα emitters (LAEs) with spectroscopic redshifts from the MUSE survey (Inami et al. 2017) and Lyman break galaxies (LBGs; Bouwens et al. 2015). The Lyα luminosities of the LAEs are L_Lyα  = (0.7−1.5) × 10⁴² erg s⁻¹, with UV-based SFR < 4 M_⊙ yr⁻¹ and stellar mass of $\mathrm{log}({M}_{* }/{M}_{\odot })$ = 8.04–8.75 (note that only two of the six MUSE LAEs have stellar mass constraints). The LBGs have H-band magnitudes of H₁₆₀ = 27.5–30.9 mag, corresponding to a UV-based SFR of 0.25–48 M_⊙ yr⁻¹, and have stellar masses of $\mathrm{log}({M}_{* }/{M}_{\odot })$ = 7.99–9.37. In general, the running mean in [Cii] luminosity of the simba sample is lower than the average of existing detections of targeted observations at z ∼ 6.

Our results are consistent with those based on the Serra simulation suite by Pallottini et al. (2017b) and Pallottini et al. (2017a), which is a suite of cosmological zoom-in adaptive mesh refinement (AMR) simulations that resolve the gas down to 10 pc scales at z ≃ 6. It is also in reasonably good agreement with the sample of 30 mufasa galaxies analyzed in Olsen et al. (2017), thus broadly confirming these results using a larger sample from its successor simulation (simba), while reaching comparable resolution over cosmological volumes. That said, our results yield a flatter L_{[C ii]}–SFR relation than other models at z ≃ 6, such as those based on SAMs and semiempirical models (Vallini et al. 2015; Lagache et al. 2018; Popping et al. 2019). We discuss the potential causes of the differences seen between our results and other models in the literature in Section 5.

4.2. [Cii] Luminosity Function at z ≃ 6

In Figure 5, we show predictions for the [Cii] LF at z ≃ 6 based on the simulated L_{[C ii]} of galaxies in Simba-25, Simba-50, and Simba-100. We note that previously LF predictions were only possible in models that made more simplified assumptions to connect L_{[C ii]} to dark matter halos. Nonetheless, our results are in agreement with those based on SAMs by Lagache et al. (2018) and Popping et al. (2019), and with constraints from the latest limits placed using data from ASPECS (Uzgil et al. 2020, in preparation; see footnote 13).

Zoom In Zoom Out Reset image size

Using the Capak et al. (2015) targeted sample, Hemmati et al. (2017) report a volume density that is almost an order of magnitude higher than our results at the bright end at $\mathrm{log}$ L_{[C ii]} ≳ 8.5 L_⊙, although they are consistent within the error bars. Note that the actual uncertainties on the [Cii] LF constrained by the Capak et al. (2015) sample are likely to be larger than those reported by Hemmati et al. (2017), as incompleteness and selection bias are not corrected for. As shown in Figure 3, most of the simba galaxies have L_{[C ii]} = 10⁶–10⁸ L_⊙. Thus, it is unsurprising to see a discrepancy in the [Cii] LF between the Capak et al. (2015) sample and our sample owing to the lack of overlap in terms of L_{[C ii]}.

Miller et al. (2020) derive a [Cii] LF using the Bolshoi–Planck dark-matter-only simulation catalog from Behroozi et al. (2013), the abundance-matched SFR from Hayward et al. (2013), and the empirical L_{[C ii]}–SFR relations established at z ∼ 0 by De Looze et al. (2014). Consistent with our results, Miller et al. (2020) report a [Cii] LF that underpredicts the observational constraints placed by Hemmati et al. (2017) using the Capak et al. (2015) sample and that placed based on a blind search of five deep fields centered on IR-bright galaxies and quasar host galaxies at z ∼ 6. By simulating only regions with L_{[C ii]} matched to the central galaxies observed in the deep fields (8.7 $\lt \,\mathrm{log}$ L_{[C ii]}/L_⊙ < 9), Miller et al. (2020) find a good agreement between the [Cii] LF and the observational constraints. On this basis, they argue that the [Cii]-detected sources in the deep fields are indeed in biased overdense regions. This may partially explain the discrepancy seen between the observed and the simba-based [Cii] LF—since the largest simulation box of simba is 100 cMpc, and it may not contain these rare highly biased regions.

4.3.  L_{[C ii]}–M_halo Relation

Forecasts for upcoming [Cii] LIM surveys have been obtained using scaling relations between [Cii] luminosity and halo mass, where the latter quantity is obtained from large-volume N-body simulations (e.g., Silva et al. 2015; Kovetz et al. 2017). The large cosmological volumes provided by N-body simulations are needed to make mock light cones for LIM (e.g., Yang et al. 2020); however, most previous work in this area has made simple empirical assumptions regarding the relationship between dark matter halo properties and [Cii] line luminosity. Here we study the L_{[C ii]}–M_halo relation based on the central galaxies in our cosmological hydrodynamic simulation.

In Figure 6, simba galaxies are shown in L_{[C ii]} versus M_halo with a fit made following the formalism of Silva et al. (2015), where SFR is expressed in terms of M_halo:

$$\begin{eqnarray}&&\mathrm{SFR}={M}_{0}\times {\left(\displaystyle \frac{{M}_{\mathrm{halo}}}{{M}_{a}}\right)}^{a}{\left(1+\displaystyle \frac{{M}_{\mathrm{halo}}}{{M}_{b}}\right)}^{b}.\end{eqnarray} \tag{ 7 }$$

Together with the expected linear relation between L_{[C ii]} and SFR, the following equation relates M_halo to L_{[C ii]}:

$$\begin{eqnarray}&&\mathrm{log}{L}_{[{\rm{C}}{\rm\small{II}}]}={M}_{0}^{{\prime} }+a^{\prime} \mathrm{log}\left(\displaystyle \frac{{M}_{\mathrm{halo}}}{{M}_{a}^{{\prime} }}\right)+b^{\prime} \mathrm{log}\left(1+\displaystyle \frac{{M}_{\mathrm{halo}}}{{M}_{b}^{{\prime} }}\right),\end{eqnarray} \tag{ 8 }$$

where L_{[C ii]} is in units of L_⊙ and SFR is in units of M_⊙ yr⁻¹. With the current set of parameters adopted in the subgrid model (see Section 3), the fitted parameters are ${a}^{{\rm{{\prime} }}}=0.97$, ${b}^{{\rm{{\prime} }}}=-9.85$, ${M}_{0}^{{\rm{{\prime} }}}=3.62$, ${M}_{a}^{{\prime} }$ = 1.50  × 10⁷, and ${M}_{b}^{{\prime} }$ = 2.90  × 10¹³.

Zoom In Zoom Out Reset image size

The L_{[C ii]}–M_halo relation of central galaxies in simba has a scatter of ≃0.5 dex around the fit. We also show the running mean in the same figure that follows the parametric form. In the same figure, we show a comparison with the L_{[C ii]}–M_halo relation from a SAM (Popping et al. 2019; Yang et al. 2020), which is steeper than both Silva et al. (2015) and this work but is, however, consistent within the scatter of simba galaxies.

As discussed in Section 4, our predicted [Cii] LF is consistent with that predicted from other models based on SAMs, and the simba galaxies lie in the same region of L_{[C ii]}–SFR as the galaxies studied by one of the most detailed cosmological zoom-in AMR simulations at the same redshift (Pallottini et al. 2017b, 2017a). However, compared to other models, our model underpredicts the [Cii] luminosity at the bright end (and high SFR; Figure 4) and yields a flatter L_{[C ii]}–SFR relation. This discrepancy may arise from, for instance, different ranges of properties for galaxies in the samples, different predicted scaling relations between galaxy properties in different galaxy evolution models, a limited number of massive halos in our simulation, and the different subgrid treatments of physical processes in simba and in sígame compared to other approaches adopted in the literature used for comparison here (Table 2). For instance, Vallini et al. (2015) simulate the [Cii] line emission by post-processing the UV radiation field of an smooth particles hydrodynamics (SPH)-based simulated galaxy using Licorice and calculate the [Cii] emission using a combination of an analytical model and the PDR code ucl_pdr (Bayet et al. 2009 and references therein; to account for contributions from PDR). While the subgrid modeling of Popping et al. (2019) follows a similar approach to that of sígame, the former is based on a SAM while the latter is applied to hydrodynamical simulations. As a result, there are relatively subtle differences between the two, such as the assumption of exponential gas disks for all galaxies in the former. In addition, the former approach assumes that all GMCs in each galaxy share the same metallicity based on the global metallicity and adopts despotic instead of cloudy in performing the thermochemistry calculation. As mentioned in Section 1, despotic does not account for line emission in the ionized phase, while the latter does. In contrast to Popping et al. (2019), Lagache et al. (2018) use cloudy to post-process their SAM. While using cloudy is the same approach as in sígame (including this work) and works by, e.g., Katz et al. (2019) and Pallottini et al. (2019), Lagache et al. (2018) adopt different subgrid approaches and assumptions compared to those made for hydrodynamical simulations. In particular, their model does not account for [Cii] emission coming from regions outside of PDRs, the complexity of the multiphase ISM, and the detailed structures of molecular clouds. This illustrates the various differences between existing models that can contribute to the discrepant L_{[C ii]}–SFR slopes. Using SAMs, Popping et al. (2019) experiment with different assumptions made in the subgrid approaches and indeed find differences in the resulting [Cii] luminosity.

Lagache et al. (2018) report that after selecting galaxies with the same range of stellar mass, SFR, and gas-phase metallicities as the Mufasa sample studied by Olsen et al. (2017)—i.e., with M_* ∈  (0.7–8) × 10⁹ M_⊙, SFR ∈ [3–23] M_⊙ yr⁻¹, and Z_gas ∈ [0.15–0.45] Z_⊙—the L_{[C ii]}–SFR relation of their SAM galaxies is flatter and more consistent with that of Olsen et al. (2017). Yet this relation remains flatter than that obtained when applying the same set of criteria to the SAM galaxies of Popping et al. (2019) (9,653 galaxies after selection). In Figure 7, we show the L_{[C ii]}–SFR for a subset of simba galaxies selected based on these criteria (1136 galaxies) compared to Olsen et al. (2017) and Popping et al. (2019). A relation with a flatter slope than the Popping et al. (2019) SAM-based models persists for the simba galaxies.

Zoom In Zoom Out Reset image size

As shown and discussed in the literature (e.g., Vallini et al. 2015; Olsen et al. 2017; Lagache et al. 2018; Pallottini et al. 2019; Popping et al. 2019), a lower metallicity can result in a lower L_{[C ii]} at given SFR. In fact, as shown in Figure 8, while the distributions in SFR between the subset of galaxies from the Popping et al. (2019) sample and this work are comparable after applying the Olsen et al. (2017) selection cut, the former are more metal-rich, with higher molecular gas masses compared to the latter. The trend of decreasing molecular gas mass and metallicity with L_{[C ii]} is also seen in Figures 9 and 10, where we show the different L_{[C ii]} distributions when the full set of 11,137 simba galaxies are selected using the first and third quartiles in molecular gas mass and metallicity, respectively. To account for the effect of SFR in L_{[C ii]}, we bin the galaxies into three SFR bins, ∈[0.01, 10.3], ∈[0.3, 10.3], and ∈[10.3, 329] M_⊙ yr⁻¹. ¹⁴ At least in the lowest SFR bin, the variation seen in the [Cii] luminosity is mostly driven by the lower metallicity and molecular gas mass (see also Narayanan & Krumholz 2017). Across all SFRs, galaxies with lower L_{[C ii]} correspond to those with the least molecular gas mass, metallicity, and SFR surface density (Figures 9, 10, and 11).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As mentioned in Section 1, models in the literature adopt different approaches and assumptions in post-processing the simulations. This could also yield different simulated line luminosities. A steeper slope, more compatible with other models in the literature, namely, Vallini et al. (2015) and Popping et al. (2019), would result from adopting the local ISM abundance ratios ¹⁵ instead of the abundance pattern tracked in simba (see Section 3.2 and Olsen et al. 2017). Both Vallini et al. (2015) and Popping et al. (2019) adopt a solar abundance pattern (of relevance to this work is the C/H ratio) and scale the C/H ratio according to the gas-phase metallicity of each galaxy to determine the [Cii] emissivity. That is, their models do not consider abundance patterns that differ from solar. The resulting L_{[C ii]} could differ significantly owing to the amount of cooling via C⁺. On the other hand, Lagache et al. (2018) accounted for abundance ratios that differ from solar. Specifically, they scale the element abundances based on the median of their sample. Unsurprisingly, the slopes of the L_{[C ii]}–SFR relation of mufasa and simba galaxies are the most compatible to the Lagache et al. (2018) model. As mentioned by Lagache et al. (2018) and Katz et al. (2019), details in modeling the interstellar radiation field intensity, self-shielding, and the choice and implementation of stellar feedback are all effects that can cause differences between existing models (see also Pallottini et al. 2017b).

In contrast to mufasa, where dust is not tracked in the simulation, causing Olsen et al. (2017) to adopt a constant DTM ratio, simba tracks dust in the simulation; however, we do not create a different cloudy lookup table for each DTM ratio since it would become computationally intractable. The mean DTM of simba galaxies is ξ_DTM = 0.25, which yields a [Cii] luminosity ≲0.5 dex higher than for models with a DTM ratio set to the solar value of ξ_DTM = 0.46 (see Appendix C; see also Olsen et al. 2017, who find that L_{[C ii]} only increased by ∼0.15 dex at a given SFR when decreasing the DTM by a factor of 2 from solar DTM for their Mufasa SFMS sample ¹⁶ ). Note that, in reality, the DTM ratio varies from galaxy to galaxy (and within galaxies themselves) and is correlated with the galaxy metallicity (e.g., Rémy-Ruyer et al. 2014; Popping et al. 2017; De Vis et al. 2019). Since the standard deviation of the DTM ratio of the simba sample studied here is σ ∼ 0.15, we do not expect the simulated [Cii] luminosity to deviate more than 0.5 dex as a result of variations in the DTM ratio.

A final important caveat is that the effect of AGNs was not included in the present modeling with cloudy, although cloudy does have the capability to do so, and this has been shown to be relevant at least for CO line emission at high redshift (Vallini et al. 2019), and hence likely also relevant for [Cii]. The effect of AGNs is an additional feature that we wish to include in the future.

In this work, we presented the first prediction of the [Cii] LF during the EOR (z ≃ 6) based on large-volume cosmological hydrodynamical simulations (simba) coupled with RT and line spectral synthesis calculations. We simulate the [Cii] line luminosity for a sample of 11,137 galaxies identified in the combined 25, 50, and 100 cMpc h⁻¹ boxes (with 2 × 1024³ particles each) at z ≃ 6. The runs for the three simulation boxes have identical input physics and produce converged GSMF and SFR functions without fine-tuning the parameters in the subgrid models of simba. In addition, both GSMF and SFR functions are in good agreement with observations at this redshift. This is crucial, as we make predictions for the [Cii] LF in the luminosity range of 5.5 < log(L_{[C ii]}/L_⊙) < 8.5 by combining the boxes.

We use an updated version of sígame to post-process the simba output. Three major improvements are implemented relative to the previous version of sígame presented in Olsen et al. (2017) to produce the results presented: (i) We do not fix the number and width of shells for each GMC model, but instead allow cloudy to determine the optimal quantities to ensure convergence. This modification leads to more accurate calculation of the grain photoelectric heating of the gas and increases the importance of gas heating due to this mechanism in the GMC models—the main excitation mode for [Cii] emission. (ii) The number of cloudy models used to create lookup tables is significantly increased from 1296 to 4096 models to ensure better sampling of the physical properties of the ISM. (iii) A DTM ratio of 0.25 is adopted based on the mean value of the simba galaxy sample rather than using a solar DTM value. Finally, Olsen et al. (2017) simulated line emission for a subset of 30 zoom-in Mufasa galaxies along the SFMS, whereas with simba we are able to expand the parameter space in M_halo, SFR, M_*, M_gas, SFR surface density, and metallicity at comparable resolution without the need for "zooming in" on specific galaxies.

We summarize the main results of this paper in the following:

• 1.  
  The simulated L_{[C ii]} is consistent with the range observed in z ∼ 6 galaxies, with a spread of ≃0.3 dex at the high SFR end of >100 M_⊙ yr⁻¹, which increases to ≃0.6 dex at the lower end of the SFR. The predicted L_{[C ii]}–SFR is consistent with targeted observed samples within the uncertainties due to selection and incompleteness effects. On the other hand, our model does not produce galaxies with values of L_{[C ii]} as high as those for some galaxies observed in targeted heterogenous samples reported in the literature, at a given SFR.
• 2.  
  The [Cii] LF is consistent with the upper limits placed by the only existing untargeted flux-limited [Cii] survey at the EOR (ASPECS) and those predicted by semianalytic models.
• 3.  
  Our model yields an L_{[C ii]}–SFR relation similar to mufasa (Olsen et al. 2017) but is flatter compared to some other models in the literature. The flatter slope results from different galaxy properties sampled by different simulations, implementation, and assumptions of subgrid recipes between different simulations and the post-processing steps (see Table 2).
• 4.  
  At a fixed SFR, galaxies with higher molecular gas mass, metallicity, and SFR surface density have higher [Cii] luminosity.
• 5.  
  We present the L_{[C ii]}–M_halo relation for the central galaxies in simba at z ∼ 6. Our relation is steeper than those based on N-body simulations and SAMs but is consistent within the scatter of ≃0.5–0.6 dex.

The differing results presented in the literature on simulating [Cii] line emission at the EOR highlight the challenges modelers face in this field. As discussed in this paper, SAMs are computationally efficient and can simulate the line emission for a statistically significant sample of galaxies, but SAMs have their limitations. They do not contain information regarding the structure of the ISM of galaxies, or the 3D distribution and morphology of galaxies, to name a few. Cosmological hydrodynamic simulations, on the other hand, can provide more detailed information on the temperature and density of the IGM and ISM, the 3D structures of galaxies, and the local properties of galaxies (e.g., each gas element has a different metallicity) but are computationally demanding. In addition, large-volume cosmological simulations still lack the resolution needed to resolve the multiphase ISM in detail. Our convergence tests among different resolution simulations indicate good convergence in stellar masses and SFRs but less-than-ideal convergence in the molecular gas fractions, indicating that some subgrid models may still need refinement to improve convergence properties. Zoom-in simulations have been used to achieve higher resolution, but the computational cost limits the number of galaxies (and thus the galaxy parameter space studied) that can be "resimulated" with the zoom-in approach in reasonable time. Calibrating parameters based on higher-resolution simulations with resolved ISM properties will be necessary to test the subgrid models and assumptions made, for example, adopting the distributions based on parsec-scale hydrodynamic simulations (e.g., Tress et al. 2020), and such work is underway by sígame group members.

Deep galaxy surveys over the past decade have provided constraints on the cosmic star formation history and SMBH growth history (e.g., Madau & Dickinson 2014; Hickox & Alexander 2018; Aird et al. 2019 ), while surveys and intensity mapping experiments planned in the next decade, with facilities such as the James Webb Space Telescope, WFIRST, Euclid, LSST, CONCERTO, HERA, SPHEREx, EXCLAIM, and TIM (see reviews by Kovetz et al. 2017; Cooray et al. 2019), will expand and sharpen our view of galaxy evolution by discovering galaxies in new ways, obtaining photometric and spectroscopic redshifts, and probing a wider parameter space in galaxy properties and the large-scale environment. Due to the brightness of the fine-structure line [Cii] and its accessibility at high redshift, it is one of the main spectral lines that may be observed at the EOR to study the ISM properties of galaxies and to secure their spectroscopic redshifts. LIM experiments such as CONCERTO, TIME, and CCAT-p will measure the [Cii] line power spectrum from galaxies at the EOR; however, the detection limits and interpretation of the observed power spectrum depend on the [Cii] line luminosities of different EOR galaxy populations within the volume sampled. The fainter populations are below the current detection limit of existing facilities, but their signal can be predicted, highlighting the importance of building theoretical frameworks to simulate the [Cii] line at the EOR using cosmological simulations.

We thank Bade Uzgil for providing the relevant numbers from the ASPECS observations and the referee for providing constructive comments that improved the clarity of this manuscript. T.K.D.L acknowledges support from the Simons Foundation and the hospitality of the Cosmic DAWN Center and the Danmarks Tekniske Universitet (DTU-Space). The Flatiron Institute is supported by the Simons Foundation. T.K.D.L, K.P.O, and T.R.G thank NORDITA for its hospitality during the NORDITA 2019 program "Zoom-in and Out: from the Interstellar Medium to the Large Scale Structure of the Universe." This research has made use of NASA's Astrophysics Data System Bibliographic Services. We acknowledge use of the Python programming language (Van Rossum & de Boer 1991), Astropy (Astropy Collaboration et al. 2013), Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), SciPy (Jones et al. 2001), and yt (Smith et al. 2009; Turk et al. 2011). We thank Robert Thompson for developing the Python package caesar.

The left panel of Figure 12 is an alternative way to show the information displayed in the middle panel of Figure 1, by replacing the plotting method with hexbin contours and showing molecular gas mass on the y-axis instead of molecular to stellar mass fraction. The median molecular gas mass is shown on top of the hexbin contours for each simulation box separately. In Figure 1, the different simulation box sizes seem inconsistent with each other, but when looking at Figure 12, the discontinuities largely go away, and we instead see a slightly increasing molecular gas mass fraction with stellar mass. Unfortunately, we do not properly cover the stellar mass range >10¹⁰ M_⊙ to compare with observations showing a slight downward trend of M_mol/M_* fractions with stellar mass in this high mass range (Saintonge et al. 2011). The right panel of Figure 12 shows how the [Cii] luminosity per molecular gas mass corresponds to a common [Cii] molecular gas mass conversion factor across the different simba volumes. The use of L_{[C ii]} as a gas mass tracer with these data is currently being investigated in a separate paper.

Zoom In Zoom Out Reset image size

In cloudy, one can supply the total metallicity, where the abundance of each element is then scaled assuming the solar composition. We account for abundance patterns of galaxies that differ from solar using the abundance ratio for each element tracked in simba (i.e., He, C, N, O, Ne, Mg, Si, S, Ca, and Fe). Since the mass fraction of each element varies as a function of metallicity, we determine its running means for each element using all the gas fluid elements of the simba galaxies studied here. This yields a function that maps a given metallicity to an abundance pattern (see Figure 13), such that for a given metallicity of the clouds, the relative elemental abundances in the cloudy input are scaled to reflect the (average) abundance ratio from simba. For more details on this procedure, we refer interested readers to Olsen et al. (2017).

Zoom In Zoom Out Reset image size

In Figure 14, we show the effects of adopting a different DTM ratio on the L_{[C ii]}–SFR relation, in particular, using a solar DTM following Olsen et al. (2017) (ξ_DTM = 0.46) instead of the mean DTM ratio found in the simba sample (ξ_DTM = 0.25). Using the former ratio, the [Cii] luminosity is ≲0.5 dex lower (see Olsen et al. 2017, who find that the L_{[C ii]} is only ∼ 0.15 dex lower at the same SFR for their Mufasa SFMS sample).

Zoom In Zoom Out Reset image size