A Self-consistent Framework for Multiline Modeling in Line Intensity Mapping Experiments

TL;DR

This paper introduces a physically-motivated, self-consistent modeling framework for line intensity mapping that connects galaxy infrared luminosity to dark matter halos, enabling better interpretation of LIM data across multiple lines.

Contribution

It develops a new modeling framework for LIM that incorporates analytic galaxy-halo connections and applies it to multiple emission lines, improving physical inference from LIM observations.

Findings

Models consistent with current observational constraints.

Demonstrates extraction of ISM properties from auto and cross-correlation analyses.

Shows potential for studying ISM evolution over cosmic time.

Abstract

Line intensity mapping (LIM) is a promising approach to study star formation and the interstellar medium (ISM) in galaxies by measuring the aggregate line emission from the entire galaxy population. In this work, we develop a simple yet physically-motivated framework for modeling the line emission as would be observed in LIM experiments, which is done by building on analytic models of the Cosmic Infrared Background that connect total infrared luminosity of galaxies to their host dark matter halos. We present models of the HI 21cm, CO(1-0), [C II] 158$\mu$m, and [N II] 122 and 205$\mu$m lines consistent with current observational constraints. With four case studies of various combinations of these lines that probe different ISM phases, we demonstrate the potential for reliably extracting physical properties of the ISM, and the evolution of these properties with cosmic time, from auto and cross-correlation analysis of these lines as measured by future LIM experiments.

Tables

Table 1. Table 1: Fiducial Parameters of CIB Model

Parameter	Description	Value	Reference
$L_0$	$L_{\mathrm{IR}}$ normalization	$0.0135L_{\odot }/M_{\odot }$	Eq. 1
$s$	$z$ evolution of $L_{\mathrm{IR}}$	3.6	Eq. 2
$T_0$	$T_{\mathrm{dust}}$ at $z=0$	$24.4$ K	Eq. 4
$\alpha$	$z$ evolution of $T_{\mathrm{dust}}$	0.36	Eq. 4
$\beta$	RJ-side index	1.75	Eq. 5
$\gamma$	Wien-side index	1.7	Eq. 5, 6
$M_{\mathrm{eff}}$	effective halo mass	${10}^{12.6}{M}_{\odot }$	Eq. 8
${\sigma }_{L/M}$	log scatter	0.5	Eq. 8

Table 2. Table 2: Physical Parameters of the Reference ISM Model.

Signal	Parameter	Symbol	Value
–	Molecular gas fraction$\star$	$f_{H_2}$	0.2
–	Ionized gas fraction$\star$	$f_{H\mathrm{ii}}$	0.1
CO	$L$–$M$ conversion$\star$	${\alpha }_{\mathrm{CO}}$	$\frac{4.4M_{\odot }}{K\mathrm{km}{s}^{-1}{\mathrm{pc}}^2}$
	Excitation temperature	$T_{\mathrm{exc}}$	10 K
	$H_2$ number density	$n_{H_2}$	$2\times {10}^3{\mathrm{cm}}^{-3}$
[C ii]	PE efficiency$\star$	${ϵ}_{\mathrm{PE}}$	$5\times {10}^{-3}$
[N ii]	Gas temperature$\star$	$T_{\mathrm{gas},H\mathrm{ii}}$	${10}^4$ K
	Electron number density$\star$	$n_{e,H\mathrm{ii}}$	${10}^2{\mathrm{cm}}^{-3}$

Table 3. Table 3: Reference Instrumental Parameters for Case Studies

Parameter	$t_{\mathrm{obs}}$	$D_{\mathrm{ap}}$	${\Omega }_{\mathrm{survey}}$	$N_{\mathrm{feeds}}$	{${\nu }_{\min }$, ${\nu }_{\max }$}	$\Delta z$	${\delta }_{\nu }$	$V_{\mathrm{vox}}$	NEFD ($T_{\mathrm{sys}}$)	${\mathrm{SNR}}_{\mathrm{tot}}$
Units	[hr]	[m]	[deg2]		[GHz]		[GHz]	[Mpc3]	[$\mathrm{mJy}{s}^{1/2}$] ([K])
	Case I: $P_{Hi}$ and $P_{\mathrm{CO}}$ at $z=2$
H i	2000	1000	1	100	{0.44, 0.52}	$\pm 0.25$	0.003	22.8	(50)	7.6
CO	500	10	1	100	{35, 42}	$\pm 0.25$	0.05	7.11	(40)	7.7
	Case II: $P_{C\mathrm{ii}}$, $P_{Hi}$, $P_{\mathrm{CO}}$, $P_{C\mathrm{ii},Hi}$ $P_{Hi,\mathrm{CO}}$, and $P_{C\mathrm{ii},\mathrm{CO}}$ at $z=2$
H i	2000	1000	1	100	{0.44, 0.52}	$\pm 0.25$	0.003	22.8	(50)	7.6, H i$\times$[C ii]: 5.0
[C ii]	1000	12	1	100	{585, 691}	$\pm 0.25$	2	0.04	50	8.0, H i$\times$[C ii]: 6.6
CO	500	10	1	100	{35, 42}	$\pm 0.25$	0.05	7.11	(40)	7.7, H i$\times$CO: 5.1
	Case III: $P_{122}^{N\mathrm{ii}}$, $P_{205}^{N\mathrm{ii}}$ and $P_{122\times 205}^{N\mathrm{ii}}$ at $z=2$
[N ii] 122	2000	10	1	400	{757, 894}	$\pm 0.25$	1	0.02	10	5.2
[N ii] 205	2000	10	1	400	{450, 531}	$\pm 0.25$	1	0.07	10	4.7, cross: 6.7
	Case IV: $P_{C\mathrm{ii}}$ and $P_{205}^{N\mathrm{ii}}$ at $z=2$
[C ii]	1000	12	1	100	{585, 691}	$\pm 0.25$	2	0.04	50	8.0
[N ii] 205	2000	10	1	400	{450, 531}	$\pm 0.25$	1	0.07	10	4.7

Equations

$$L_{\mathrm{IR},\left(1+z\right)\nu }\left(M,z\right)=L_{\mathrm{IR},0}\Phi \left(z\right)\Sigma \left(M\right)\Theta \left[\left(1+z\right)\nu \right],$$

$$\Phi \left(z\right)={\left(1+z\right)}^s,$$

$$I_{\nu }\propto {\nu }^{\beta }{B}_{\nu }\left[T_{\mathrm{d}}(z)\right],$$

$$T_{\mathrm{d}}(z)=T_0(1+z{)}^{\alpha },$$

$$\Theta \left[\left(1+z\right)\nu \right]=A\left(z\right)\times \{\begin{array}{ll}{\nu }^{\beta }{B}_{\nu }\left[T_{\mathrm{d}}(z)\right]& \text{ν<ν0}\\ {\nu }^{-\gamma }& \text{ν≥ν0}\end{array},$$

$$\frac{\mathrm{d}\ln \left\{{\nu }^{\beta }{B}_{\nu }\left[T_{\mathrm{d}}(z)\right]\right\}}{\mathrm{d}\ln \nu }=-\gamma$$

$$\int \Theta \left(\nu ,z\right)d\nu =1$$

$$\Sigma(M)=M\frac{1}{\sqrt{2\pi\sigma^{2}_{L/M}}}\mathrm{exp}\Big{[}-\frac{(\mathrm{log}_{10}M-\mathrm{log}_{10}M_{\mathrm{eff}})^{2}}{2\sigma_{L/M}^{2}}\Big{]}~{},$$

$$L_{\mathrm{IR}}(M,z)={\int }_{300\mathrm{GHz}}^{37.5\mathrm{THz}}\mathrm{d}\nu L_{(1+z)\nu }(M,z).$$

$${\dot{M}}_{\star }(M,z)={\mathcal{K}}_{\mathrm{IR}}{L}_{\mathrm{IR}},$$

$${\dot{\rho }}_{\star }(M,z)={\int }_{M_{\min }}^{M_{\max }}\frac{\mathrm{d}N}{\mathrm{d}M}{\dot{M}}_{\star }(M,z),$$

$$L_{\mathrm{IR}}\left(M,z\right)=P_0{M}_{\mathrm{d}}\left(M,z\right){\left[\frac{T_{\mathrm{d}}(z)}{T_0}\right]}^{4+\beta },$$

$$P_0$$

$$\simeq 110L_{\odot }/M_{\odot },$$

$$M_{\mathrm{d}}\propto \Sigma (M)(1+z{)}^{s-\alpha (4+\beta )}.$$

$${\Omega }_{\mathrm{d}}\left(z\right)=\frac{1}{{\rho }_{\mathrm{crit},0}}\int \mathrm{d}M\frac{\mathrm{d}N}{\mathrm{d}M}{M}_{\mathrm{d}}(M,z),$$

$$M_{\mathrm{H}}(M,z)=K(z)\Sigma (M)\Phi (z),$$

$$M_{\mathrm{H}}=M_{\mathrm{H}\text{mathsc}\mathrm{i}}+M_{{\mathrm{H}}_2}+M_{\mathrm{H}\text{mathsc}\mathrm{ii}}.$$

$$M_{{\mathrm{H}}_2}(M,z)$$

$$M_{\mathrm{H}\text{mathsc}\mathrm{ii}}(M,z)$$

$$M_{\mathrm{H}\text{mathsc}\mathrm{i}}(M,z)$$

$$\frac{Z}{Z_{\odot }}(z)\sim 100\frac{M_{\mathrm{d}}(M,z)}{M_{\mathrm{H}}(M,z)}.$$

$${\bar{T}}_{\mathrm{H}\text{mathsc}\mathrm{i}}={\mathcal{C}}_{\mathrm{H}\text{mathsc}\mathrm{i}}{\bar{\rho }}_{\mathrm{H}\text{mathsc}\mathrm{i}}(z)=\frac{3h{c}^3{A}_{21}}{32\pi k_{\mathrm{B}}{m}_{\mathrm{p}}{\nu }_{21}^2}\frac{(1+z{)}^2}{H(z)}{\bar{\rho }}_{\mathrm{H}\text{mathsc}\mathrm{i}}(z),$$

$${\bar{\rho }}_{\mathrm{H}\text{mathsc}\mathrm{i}}(z)=\int \mathrm{d}M\frac{\mathrm{d}N}{\mathrm{d}M}{M}_{\mathrm{H}\text{mathsc}\mathrm{i}}(M,z).$$

$$L_{[\mathrm{C}\text{mathsc}\mathrm{ii}]}=\left(1-f_{{\mathrm{H}}_2}\right){ϵ}_{\mathrm{PE}}{L}_{\mathrm{IR}},$$

$$L_{[\mathrm{N}\text{mathsc}\mathrm{ii}]}\simeq n_{\mathrm{e},\mathrm{H}\text{mathsc}\mathrm{ii}}{n}_{{\mathrm{N}}^{+}}{q}_{\nu }h{\nu }_{[\mathrm{N}\text{mathsc}\mathrm{ii}]}V,$$

$$Q_0=n_{\mathrm{e},\mathrm{H}\text{mathsc}\mathrm{ii}}{n}_{{\mathrm{H}}^{+}}{\alpha }_{\mathrm{B}}\left(T_{\mathrm{gas},\mathrm{H}\text{mathsc}\mathrm{ii}}\right)V,$$

$$Q_0(M,z)=\frac{N_{\mathrm{ion}}{\dot{M}}_{\star }}{m_{\mathrm{p}}/(1-Y)}\simeq 1.14\times 10^{53}\left[\frac{{\dot{M}}_{\star }(M,z)}{M_{\odot }/\mathrm{yr}}\right]{\mathrm{s}}^{-1},$$

$$L_{[\mathrm{N}\text{mathsc}\mathrm{ii}]}\simeq \frac{q_{\nu }h{\nu }_{[\mathrm{N}\text{mathsc}\mathrm{ii}]}}{{\alpha }_{\mathrm{B}}\left(T_{\mathrm{gas},\mathrm{H}\text{mathsc}\mathrm{ii}}\right)}\frac{n_{{\mathrm{N}}^{+}}}{n_{{\mathrm{H}}^{+}}}\frac{N_{\mathrm{ion}}{\dot{M}}_{\star }}{m_{\mathrm{p}}/(1-Y)}$$

$$L_{[\mathrm{N}\text{mathsc}\mathrm{ii}]}^{\mathrm{tot}}=9\times 10^6{L}_{\odot }{\left(\frac{T_{\mathrm{gas},\mathrm{H}\text{mathsc}\mathrm{ii}}}{10^4\mathrm{K}}\right)}^{0.76}\times \frac{{\dot{M}}_{\star }}{M_{\odot }/\mathrm{yr}}\times \frac{Z}{Z_{\odot }},$$

$$R_{[\mathrm{N}\text{mathsc}\mathrm{ii}]}^{\prime }=R_{[\mathrm{N}\text{mathsc}\mathrm{ii}]}+A_R\times \left\{1+\mathrm{erf}\left[\frac{\log (n_{\mathrm{e},\mathrm{H}\text{mathsc}\mathrm{ii}}/n_{\mathrm{e},\mathrm{H}\text{mathsc}\mathrm{ii}}^0)}{{\sigma }_R}\right]\right\},$$

$${L^{\prime }}_{[\mathrm{N}\text{mathsc}\mathrm{ii}]}^{\mathrm{tot}}\simeq L_{[\mathrm{N}\text{mathsc}\mathrm{ii}]}^{\mathrm{tot}}\left(T_{\mathrm{gas},\mathrm{H}\text{mathsc}\mathrm{ii}}\right)\times 10^{0.12\left[R_{[\mathrm{N}\text{mathsc}\mathrm{ii}]}-R_{[\mathrm{N}\text{mathsc}\mathrm{ii}]}^{\prime }(n_{\mathrm{e},\mathrm{H}\text{mathsc}\mathrm{ii}})\right]}.$$

Full text

cosmology: observations – cosmology: theory – galaxies: ISM – infrared: diffused background – large-scale structure of universe

1 Introduction

Line intensity mapping (LIM) is an emerging observational technique developed to statistically measure the intensity field fluctuations of a given spectral line (see Kovetz et al. 2017 for a recent review). While traditional galaxy surveys are restricted by the detection limit of individual sources, LIM is sensitive to the emission from all galaxies, providing a complementary probe of faint objects. Due to its statistical nature, LIM is most effective at constraining how average physical properties, including the star formation rate, ISM conditions, luminosity function, spatial distribution, etc., of the source galaxy population evolve over cosmic time (Serra et al. 2016, hereafter S16; Kovetz et al. 2017; Chang et al. 2019).

LIM was first pioneered with the redshifted H i 21 cm line signal. It serves as a probe of both the matter density distribution as traced by the atomic hydrogen gas in the interstellar medium (ISM), for example, the baryon acoustic oscillation (BAO) feature in galaxy power spectrum (Chang et al., 2010; Switzer et al., 2013), and the structure of neutral intergalactic medium (IGM) at high redshift, in particular during cosmic reionization (Madau et al., 1997; Furlanetto et al., 2004, 2006; Pritchard & Loeb, 2012). Recently, the application of LIM to other emission lines has gained increasing attention, including CO rotational lines (Pullen et al., 2013; Breysse et al., 2014; Mashian et al., 2015; Li et al., 2016), far-infrared (FIR) fine-structure lines of C ii, N ii, O i and others (Gong et al. 2012; Uzgil et al. 2014; Silva et al. 2015; Yue et al. 2015; S16), and bright optical/UV emission lines such as Ly$\alpha$ and H$\alpha$ (Silva et al., 2013; Pullen et al., 2014; Comaschi & Ferrara, 2016; Gong et al., 2017; Silva et al., 2018).

Substantial theoretical and experimental efforts have been devoted to the detection and interpretation of LIM signals of individual lines. However, a simple, physical model that allows multiple line signals, presumably originating from and thus probing different ISM phases, to be modeled in a self-consistent manner is still lacking. The goal of this work is to develop such a self-consistent framework for determining the integrated line intensities of galaxies observed in the intensity mapping regime. This framework is intended to bridge the gap between commonly-used approaches anchored on scaling relations empirically determined from observations (e.g., Visbal & Loeb 2010; Pullen et al. 2013; Silva et al. 2015; Li et al. 2016; S16) and sophisticated simulations of galaxy-scale hydrodynamics and radiative transfer (e.g., Pallottini et al., 2019; Popping et al., 2019). More specifically, it should be sophisticated enough to capture the relevant ISM physics and employ meaningful physical parameters, yet simple enough to interpret the auto/cross-correlation of the intensities of various lines observed in the intensity mapping regime in terms of coarse-grained galactic ISM properties. Some examples include the mass fraction of different ISM phases, as well as quantities like the photoelectric (PE) heating efficiency and the CO-to-$\rm H_{2}$ ratio, which are closely related to exact physical conditions of the ISM (e.g., temperature, density, radiation field) and therefore of particular interest to LIM surveys of the corresponding lines. Furthermore, this analytical framework should also allow mock signal maps to be readily constructed from given information about the position and physical properties of source populations, thereby enabling straightforward implementation in semi-analytic models of the LIM signals.

We build such a formalism using the information from the cosmic infrared background (CIB). The CIB has, on account of sensitive FIR observations from experiments like Planck (Planck Collaboration XXX, 2014) and Herschel (Viero et al., 2013a), been the subject of detailed modeling efforts. In particular, analytic models connecting the infrared (IR) luminosity of galaxies to the mass and redshift of their host halos have been successful in reproducing the statistical properties of the CIB (e.g., Shang et al. 2012; S16; Wu & Doré 2017). In this work, we follow and extend the ideas presented in S16 by employing the CIB model as a starting point for models of both line and continuum emission from galaxies as a function of redshift and halo mass. Taken at face value, the IR luminosities assumed in these models imply a corresponding dust mass, gas mass, and metallicity, which in turn can inform predictions of emission from various interstellar lines, including H i, [C ii], [N ii] and CO(1-0). We work through these consequences, with an eye toward testable predictions from upcoming intensity mapping experiments of these lines.

This paper is organized as follows. In Section 2, we present a simple analytic model that describes a variety of physical properties of dark matter halos hosting the line-emitting galaxies, such as their star formation rate, dust mass, metallicity and so forth. We then discuss in Section 3 how we model the emission of H i, [C ii], [N ii] and CO lines as tracers of different phases of the ISM, based on our model of halo properties. In Section 4, we review the theoretical framework of estimating the power spectrum signal of intensity mapping experiments, as well as the uncertainty associated with the measurements. We compare the predicted strengths of different lines to constraints from the literature in Section 5. We then present four case studies in Section 6 to demonstrate how physical conditions of multi-phase ISM may be probed by and extracted from with intensity mapping experiments. We outline prospects for further improving and extending our simple modeling framework, before briefly concluding in Section 7. Throughout the paper, we assume a flat, $\Lambda$CDM cosmology consistent with the measurement by the Planck Collaboration XIII (2016).

2 A Simple Analytic Model of Mean Halo Properties

An important criterion for choosing first targets for LIM surveys is the overall brightness of the spectral line, which is determined by many different factors, including abundance, excitation potential, critical density, destruction and/or scattering and so forth. In many cases, nevertheless, the line signal either directly traces the star-forming activity (e.g., [C ii], [N ii]) or indirectly probes the gas reservoir closely associated with star formation (e.g., H i, CO). Therefore, it is critical to understand and model the star formation of galaxies well enough in order to properly estimate the production of lines in LIM.

The majority of starlight from young stars at optical/UV wavelengths is absorbed and reprocessed into IR radiation by interstellar dust, naturally giving rise to the connection between IR observations of galaxies and their star formation rate. Because the fraction of spatially resolved galaxies decreases rapidly with increasing wavelength in the IR/sub-millimeter regime, the observed CIB mean intensity and fluctuations provide a useful probe of global star-forming activities. Combining the halo model formalism describing the clustering of galaxies at different angular scales (Cooray & Sheth, 2002) and the observed angular anisotropy of the CIB, Shang et al. (2012) developed a simple parametric form for the infrared luminosity of galaxies as a function of halo mass and redshift, which has been successfully applied to reconstruct the observed angular CIB auto- and cross-power spectra (Planck Collaboration XXX 2014; S16; Wu & Doré 2017). In this section, we extend the discussion in S16 and present a simple, CIB-based model for the mean properties of dark matter halos, such as their infrared luminosity, dust and gas mass, metallicity, etc., which are essential ingredients for the line emission models in this work.

2.1 IR luminosity

We work in the aforementioned framework of the halo model for CIB anisotropies introduced by Shang et al. (2012), which has been exploited in various contexts, including the modeling of high-redshift emission lines (e.g., Planck Collaboration XXX 2014; S16; Wu & Doré 2017; Pullen et al. 2018). In this model, the specific luminosity emitted by a galaxy hosted by a halo of mass $M$ at redshift $z$ at the observed frequency $\nu$ is given by

[TABLE]

where $L_{{\rm IR},0}$ is a normalization constant (see Table LABEL:tb:cib for a summary of fiducial parameter values taken for the CIB model), whereas $\Phi\left(z\right)$, $\Sigma\left(M\right)$, and $\Theta\left[\left(1+z\right)\nu\right]$ are functions to be specified. $\Phi(z)$ governs the evolution of the luminosity–mass relation with redshift, driven, e.g., by an increase in the star formation rate at fixed halo mass with increasing redshift. This is modeled as a power law

[TABLE]

where Wu & Doré (2017) found a best-fit value of $s=3.6$. However, we note that the exact value of $s$ is not well-constrained by the integrated CIB intensity and less steep slopes have indeed been suggested by some other CIB analyses and galaxy evolution models (see discussion in S16).

$\Theta\left[\left(1+z\right)\nu\right]$ describes the frequency dependence of the dust emission as a function of redshift. Over most of the FIR frequency range, the dust emission in a galaxy is modeled as a modified blackbody of temperature $T_{\rm d}$ and spectral index $\beta$,

[TABLE]

where $B_{\nu}\left[T_{\rm d}(z)\right]$ is the Planck function at a dust temperature

[TABLE]

where $T_{0}=24.4$ K is the typical dust temperature in a star-forming galaxy at $z=0$, and the redshift dependence is taken to be $\alpha=0.36$ following Planck Collaboration XXX (2014) and Wu & Doré (2017). The high frequency component is modeled as a power law to account for emission from small, stochastically-heated grains. The full SED is given by

[TABLE]

where the frequency $\nu_{0}$ at any given redshift is determined by having

[TABLE]

satisfied at $\nu=\nu_{0}$. We adopt $\beta=1.75$ and $\gamma=1.7$ (Planck Collaboration XXX, 2014; Wu & Doré, 2017), which yield $\nu_{0}=3.3,2.1,1.6$ and 1.3 THz or wavelength equivalents 92, 143, 185, and 222 $\mu$m at redshifts $z=0,1,2$ and 3, respectively. The redshift-dependent normalization factor $A\left(z\right)$ is defined such that

[TABLE]

for all $z$.

$\Sigma(M)$ links the IR luminosity to the halo mass and is modeled as a lognormal relation:

[TABLE]

where $M_{\rm eff}$ describes the most efficient halo mass at hosting star formation, and $\sigma_{L/M}$ accounts for the range of halo masses mostly contributing to the infrared luminosity. This functional form captures the fact that the star formation efficiency is suppressed for halo masses much lower or much higher than $M_{\rm eff}$ (Mo et al., 2010; Furlanetto et al., 2017; Kravtsov et al., 2018), due to various feedback mechanisms such as input from supernova explosions and active galactic nuclei (AGNs). The total infrared luminosity (8–1000 $\mu$m) is then

[TABLE]

2.2 Star Formation History

From the total infrared luminosity, it is straightforward to derive the star formation rate as a function of halo mass and redshift thanks to the well-established correlation between them (Kennicutt, 1998; Madau & Dickinson, 2014). In this work, we simply take

[TABLE]

where $\mathcal{K}_{\rm IR}=1.73\times 10^{-10}\,M_{\odot}\,\mathrm{yr^{-1}}\,L^{-1}_{\odot}$, consistent with a stellar population with a Salpeter initial mass function (IMF) and solar metallicity. The star formation rate density (SFRD) can consequently be written as

[TABLE]

where $\mathrm{d}N/\mathrm{d}M$ is the dark matter halo mass function defined for the virial mass $M_{\rm vir}$ (Tinker et al., 2008). Figure 1 shows a comparison between cosmic SFRDs predicted by the adopted CIB model and those from the literature. The data points represent estimated SFRDs based on both dust-corrected UV observations (Cucciati et al., 2012) and infrared/sub-millimeter observations of obscured star formation (Gruppioni et al., 2013; Rowan-Robinson et al., 2016). Also shown is the maximum-likelihood model of cosmic SFRD from Robertson et al. (2015) based on extrapolating the galaxy IR and UV luminosity functions down to $10^{-3}L_{\star}$. The agreement between the CIB-derived SFRD and the optical/UV-derived SFRD may be improved with different modeling choices (e.g., Maniyar et al., 2018). However, this comes at the expense of phenomenological parameterizations of the effective bias factor of dusty galaxies, and we therefore do not follow that approach here.

2.3 Dust Mass

Since we have specified both the dust luminosity and the dust temperature from the CIB model, it is possible to estimate the implied dust mass. Assuming that the dust mass is dominated by larger grains whose emission can be described by a modified blackbody with a single dust temperature, the dust luminosity, mass, and temperature can be related via

[TABLE]

where the normalization constant $P_{0}$ is the power emitted per mass of dust at temperature $T_{0}$.

To estimate $P_{0}$, we note that Planck Collaboration Int. XVII (2014) found the Galactic H i-correlated dust emission to be well-described by Equation 3 with $T_{\rm d}\simeq 20$ K and $\beta\simeq 1.6$. Further, they derived an 857 GHz dust emissivity per H of $\epsilon_{857}=4.3\times 10^{-21}$ MJy sr*-1* cm2 H*-1*. Thus,

[TABLE]

where we have assumed a gas-to-dust mass ratio of 100 (Draine et al., 2007). By using this formalism to estimate the dust mass, we are implicitly assuming that physical properties (e.g., composition) of dust grains do not evolve systematically with redshift or metallicity, only their abundance per H atom.

Because $\Theta$ is normalized to unity (see Equation 5), the redshift dependence of $L_{\rm IR}$ is determined entirely by $\Phi(z)$, and thus

[TABLE]

The implied cosmic density of dust, $\Omega_{\rm d}$, as a function of redshift is

[TABLE]

where $\rho_{\rm crit,0}$ denotes the critical density of the universe at the present time. In Figure 2, we plot the redshift evolution of the dust density parameter $\Omega_{\rm d}$, which is compared with a compilation of previous dust abundance measurements by Thacker et al. (2013), including constraints from integrating low-$z$ dust mass functions (Dunne et al., 2011), extinction measurements from the Sloan Digital Sky Survey111The combined data set from Thacker et al. (2013) is adopted here, which assumes that the halo dust content does not evolve significantly with redshift. (Ménard et al., 2010; Ménard & Fukugita, 2012) and 2dF (Driver et al., 2007), and cosmic far-infrared background anisotropy (Thacker et al., 2013).

2.4 Hydrogen Mass

Insofar as gas and dust are well mixed, $\Phi\Sigma$ encodes the total hydrogen mass in the halo. However, the correspondence is not direct since the dust luminosity depends on not only the amount of dust present but also the dust temperature, which is assumed to evolve with redshift (see Equation 5). We therefore introduce the modification

[TABLE]

where $K(z)=\zeta(1+z)^{\xi}$ is a normalization factor that sets the total amplitude of $M_{\rm H}$. The amplitude and redshift dependence of $K$ are determined by approximately matching the hydrogen–halo mass relation over $0<z<3$ predicted by Popping et al. (2015) as shown in Figure 3, while at the same time yielding a gas metallicity of approximately $Z_{\odot}$ at $z=0$ (discussed in next section). For our fiducial model, we take $\zeta=0.005$, $\xi=-1$ for $0<z<1$ and $\zeta=0.0025$, $\xi=0$ otherwise. These values are chosen such that the overall redshift dependence of $M_{\rm H}$ roughly agrees with the product of inferred growth rate of halo mass, which scales as $(1+z)^{1.5}$ at $z\lesssim 1$ and $(1+z)^{2.5}$ at higher redshifts (McBride et al., 2009), and the average star formation efficiency, which may carry an extra factor of $(1+z)^{1-1.5}$ depending on the exact physical mechanisms coupling the stellar feedback (e.g., supernova explosions) to galaxies (Sun & Furlanetto, 2016; Furlanetto et al., 2017). Indeed, the gas-to-stellar mass ratio of $M>10^{11}\,M_{\odot}$ halos of interest has been found to be only weakly dependent on redshift (Popping et al., 2015). The mass dependence, on the other hand, is motivated since the same physical mechanisms preventing star formation at both ends of halo masses also play a role in regulating the hydrogen mass in a galaxy.

The total mass of hydrogen in our model can be written as

[TABLE]

If we express the fractions of molecular and ionized hydrogen as $f_{\rm H_{2}}$ and $f_{{\mathrm{H\,\mathsc{ii}}}}$ respectively, then the masses of hydrogen in three different phases become

[TABLE]

As a fiducial value, we set $f_{\rm H_{2}}=0.2$, typical for most galaxies up to $z\sim 1$ and the most massive ones up to $z\sim 2$ (see, e.g., Popping et al. 2012). Likewise, we adopt $f_{{\mathrm{H\,\mathsc{ii}}}}=0.1$, based on the estimated masses of different ISM phases from Tielens (2005). We note that a factor of 1.36 accounting for the helium abundance is needed to connect the total hydrogen mass to the total gas mass, i.e., $M_{\rm gas}=1.36M_{\rm H}$ (Draine et al., 2007).

Using Eq. 16 and the fiducial molecular gas fraction $f_{\rm H_{2}}=0.2$, we can also obtain the cosmic evolution of the molecular gas density $\rho_{\rm H_{2}}$, whose comparison against the cosmic SFRD (especially the peak of star formation at $z\sim 2$) provides vital information about the fueling and regulation of star formation by cold gas. Constraints on $\rho_{\rm H_{2}}$ have so far been placed primarily by observations of the CO rotational transitions. Figure 4 shows how $\rho_{\rm H_{2}}$ as a function of redshift, computed with our fiducial choice of $f_{\rm H_{2}}=0.2$, compares with constraints derived from various CO LIM experiment and deep galaxy surveys, including COLDz (Riechers et al., 2019), COPSS II (Keating et al., 2016), ASPECS Pilot (Decarli et al., 2016) and ASPECS large program (Decarli et al., 2019). Planned LIM experiments such as COMAP (Li et al., 2016) and TIME (Crites et al., 2014) and next-generation Very Large Array (ngVLA) concepts (Walter et al., 2019) are expected to greatly reduce the substantial uncertainties present in current limits.

2.5 Metallicity

If the dust-to-metals ratio is assumed to be constant, the metallicity $Z$ of interstellar gas in our model can be expressed as a function of the dust mass as

[TABLE]

Recent hydrodynamic galaxy formation simulations have indeed found little variation in the dust-to-metals ratio with redshift or metallicity above 0.5$Z_{\odot}$ (Li et al., 2019).

Our simple model of the gas-phase metallicity gives no halo mass dependence, which is likely an oversimplification given that effects like galactic winds regulating the metallicity of galaxies may evolve with halo mass in a non-trivial way. Nevertheless, as shown in Figure 5, our predicted redshift evolution of metallicity is broadly consistent with that estimated semi-analytically by Fu et al. (2013), using the Millennium-II Simulation (Boylan-Kolchin et al., 2009) combined with an $\rm H_{2}$ prescription specified by the gas surface density, metallicity and a constant clumping factor (Krumholz et al., 2009; McKee & Krumholz, 2010). We note that $Z$ only evolves moderately for $M>10^{11.5}\,M_{\odot}$, a halo mass range that our CIB model is calibrated and most sensitive to. Therefore, in the context of the CIB model a mass-independent gas metallicity is likely a fair approximation. Figure 5 also shows the cosmic metallicity evolution inferred from gamma-ray burst (GRB) observations for comparison. By analogy to damped Ly$\alpha$ (DLA) systems of quasars, Savaglio (2006) uses strong absorption lines due to the intervening neutral gas to estimate the metallicity evolution of GRB-DLA systems and compare it with the average metallicity derived for a sample of GRB hosts at $z<1$.

3 Models of Emission Lines

Based on the mass and redshift dependencies of a wide range of halo properties derived in Section 2, we construct a model of the emission lines that trace star formation and the ISM in galaxies. In this section, we present our line emission models of the H i 21 cm line, [C ii] 158 $\mu$m line, the 122 and 205 $\mu$m [N ii] lines, and the CO(1-0) 2.6 mm line. Each of these lines probes a somewhat different phase of the ISM, ranging from the coldest molecular gas to the warm ionized medium (WIM). As such, their joint analysis can reveal rich information about the multi-phase ISM, as will be illustrated in the following sections.

Figure 6 illustrates how our modeling framework connects the emission from each of these lines to the phases of the ISM. Young stars formed in dense regions of a giant molecular cloud (GMC) are surrounded by H ii regions ionized by UV radiation, whose physical conditions may be probed by FIR [N ii] and [C ii] lines. Photodissociation regions (PDRs) occupy the interface of H ii regions and cold molecular gas traced by CO lines and produce the majority of [C ii] emission, which is the main cooling mechanism balancing the photoelectric heating by dust grains. Together molecular gas clouds compose roughly half of the total ISM mass, whereas warm/cold atomic gas contributing most of the H i mass is responsible for the remaining half.

3.1 H i 21cm Line

The hyperfine structure H i 21 cm line serves as a direct probe of the atomic hydrogen content of galaxies, so its abundance and clustering properties can be straightforwardly modeled with the H i–halo mass relation derived. The H i mass is related to the mean brightness temperature, the relevant observable for H i maps, via (e.g., Bull et al., 2015; Wolz et al., 2017)

[TABLE]

where $\mathcal{C}_{{\mathrm{H\,\mathsc{i}}}}$ is the conversion factor from the mean H i density to the mean brightness temperature and $A_{21}=2.88\times 10^{-15}\,\mathrm{s^{-1}}$ is the Einstein coefficient corresponding to the 21 cm line. The mean H i mass density is expressed as (Padmanabhan et al., 2017)

[TABLE]

3.2 [C ii] 158 $\mu$m Line

The 158 $\mu$m [C ii] line is one of the most important metal cooling lines in the interstellar medium and can alone account for $\sim 0.1$% of the total FIR emission of a galaxy (Stacey et al., 1991; Malhotra et al., 1997). Empirically, the emission in the [C ii] line correlates with both FIR dust emission (Crawford et al., 1985; Wright et al., 1991) and star formation (Stacey et al., 1991; De Looze et al., 2014).

The strong correlation between the [C ii] and IR luminosity can be understood with a model in which the cooling of interstellar gas is dominated by [C ii] emission and the heating is dominated by photoelectric emission from dust grains. If the dust converts a fraction $\epsilon_{\rm PE}\ll 1$ of UV and optical radiation absorbed into photoelectric heating and the remainder into infrared emission, then the total heating rate is proportional to $\epsilon L_{\rm IR}$. We can therefore approximate

[TABLE]

where the factor $\left(1-f_{{\rm H}_{2}}\right)$ accounts for the fact that dust is present and will radiate in molecular clouds where there is little atomic C. $\epsilon_{\rm PE}$ is taken to be a free parameter in the model with a fiducial value of $5\times 10^{-3}$, which yields an $L_{[{\mathrm{C\,\mathsc{ii}}}]}$/$L_{\rm IR}$ ratio consistent with that estimated from observations of the LMC (e.g., Rubin et al., 2009) and nearby galaxies (e.g., De Looze et al., 2014). We note that the observed proportionality between SFR and $L_{[{\mathrm{C\,\mathsc{ii}}}]}$ is reproduced here since SFR is correlated with $L_{\rm IR}$ (Equation 10).

A number of simplifications are inherent in this prescription. For instance, other cooling lines (e.g., [O i]) can be important relative to [C ii] (Tielens & Hollenbach, 1985; Young Owl et al., 2002). Second, the photoelectric efficiency of dust grains is a function of the grain charge. As gas density and radiation intensity increase, $\epsilon_{\rm PE}$ is expected to decrease (Bakes & Tielens, 1994), and so we might expect systematic changes in the $L_{\rm IR}$–$L_{[{\mathrm{C\,\mathsc{ii}}}]}$ relation with galaxy properties just from this effect. Finally, unlike the dust emission, the [C ii] line can saturate at high gas temperatures and radiation intensities, breaking the linear correlation (Muñoz & Oh, 2016; Rybak et al., 2019). These effects are most pronounced in gas of extreme density and temperature and may account for the breakdown of the $L_{\rm IR}$–$L_{[{\mathrm{C\,\mathsc{ii}}}]}$ correlation in luminous and ultraluminous galaxies. We do not incorporate these effects into our model at this time, but we discuss potential implementation in Section 7.

3.3 [N ii] 122 and 205$\mu$m Lines

The emission from singly ionized nitrogen, which has an ionization potential of 14.53 eV, traces H ii regions (see Figure 6). When the density is lower than the critical density, collisional de-excitation can be neglected and the luminosity of the [N ii] 122 and 205 $\mu$m lines can be approximated by the balance between the rates of collisional excitation and radiative de-excitation. For an ionized gas cloud of volume $V$,

[TABLE]

where $q_{\nu}$ denotes the collisional excitation coefficient, with $q_{122}=2.57\times 10^{-8}\,\mathrm{cm^{3}\,s^{-1}}$ and $q_{205}=6.79\times 10^{-8}\,\mathrm{cm^{3}\,s^{-1}}$ (e.g., Herrera-Camus et al., 2016). Meanwhile, the ionization equilibrium of H ii regions gives

[TABLE]

where $Q_{0}$ is the rate of hydrogen photoionization sourced by UV photons from O and B stars and $\alpha_{\rm B}=2.6\times 10^{-13}\left(T_{\rm gas,{\mathrm{H\,\mathsc{ii}}}}/10^{4}\,\mathrm{K}\right)^{-0.76}\,\mathrm{cm^{3}\,s^{-1}}$ is the case B recombination coefficient, a reasonable assumption for typical H ii regions where the mean free path of ionizing photons is small. For Population II stars with a Salpeter IMF, each stellar baryon produces $N_{\rm ion}\simeq 4000$ ionizing photons on average (Loeb & Furlanetto, 2013), in which case $Q_{0}$ can be related to the star formation rate by

[TABLE]

where we take the helium mass fraction to be $Y=0.25$. The ionization rate can then be related to the luminosity of [N ii] lines by

[TABLE]

which gives

[TABLE]

where $n_{\rm N^{+}}/n_{\rm H^{+}}$, under the assumption that the second ionization of nitrogen ($\rm N^{+}\rightarrow N^{++}$) with a potential of 29.6 eV is negligible, can be approximated by the $\rm N/H$ ratio $\mathrm{N/H}=\mathrm{(N/H)_{\odot}}\times\left[Z(z)/Z_{\odot}\right]\simeq 7.4\times 10^{-5}\left[Z(z)/Z_{\odot}\right]$ (Asplund et al., 2009).

In order to model the impact of electron number density $n_{\rm e,{\mathrm{H\,\mathsc{ii}}}}$ on the strength of [N ii] line emissions, we exploit a simple parameterization of the [N ii] 122 $\mu$m/205 $\mu$m line ratio as a function of $n_{\rm e,{\mathrm{H\,\mathsc{ii}}}}$

[TABLE]

where $R_{[{\mathrm{N\,\mathsc{ii}}}]}=L^{122}_{[{\mathrm{N\,\mathsc{ii}}}]}/L^{205}_{[{\mathrm{N\,\mathsc{ii}}}]}\sim 0.65$ is the line ratio in the low-density limit discussed above. We further take the normalization factor $A_{R}$ to be 4.7, the characteristic density $n^{0}_{\rm e,{\mathrm{H\,\mathsc{ii}}}}$ to be $10^{2.5}\,\mathrm{cm^{-3}}$ and the transition width $\sigma_{R}$ to be 1, in order to obtain a good fit to the results from Herrera-Camus et al. (2016) over $1\,\mathrm{cm^{-3}}\lesssim n_{\rm e,{\mathrm{H\,\mathsc{ii}}}}\lesssim 10^{5}\,\mathrm{cm^{-3}}$, as illustrated in Figure 7. Meanwhile, there is also a non-trivial evolution of the total [N ii] luminosity with the electron number density (or effectively $R^{\prime}_{[{\mathrm{N\,\mathsc{ii}}}]}$) owing to the increasingly important collisional de-excitation at higher densities, whose effect can be approximated by

[TABLE]

The resulting [N ii] line luminosities depend on both the temperature and the density of H ii regions:

[TABLE]

and

[TABLE]

In our model of [N ii] emission, we set the gas temperature to be $T_{\rm gas,{\mathrm{H\,\mathsc{ii}}}}\simeq 10^{4}$ K, which is a characteristic temperature of ${\mathrm{H\,\mathsc{ii}}}$ regions where ionized nitrogen is expected to be found (e.g., Goldsmith et al., 2015; Herrera-Camus et al., 2016). Meanwhile, it is important to note that, alternatively to the empirical prescription presented, the dependence of [N ii] line ratio on $n_{\rm e,{\mathrm{H\,\mathsc{ii}}}}$ may also be derived ab initio from the transition rates of collisionally coupled states of [N ii] (see, e.g., Goldsmith et al. 2015).

3.4 CO(1-0) Line

The CO(1-0) rotational transition ($\lambda=2.6$ mm) is a powerful tracer of the molecular gas content of both individual molecular clouds and of galaxies (e.g., Solomon et al., 1987; Dame et al., 2001; Ivison et al., 2011; Saintonge et al., 2011). In molecular clouds, the CO(1-0) line is generally optically thick, and so the line luminosity $L_{\rm CO}$ is independent of the CO abundance. For a virialized molecular cloud, it can be shown that $L_{\rm CO}$ is proportional to the cloud mass (e.g., Draine, 2011; Bolatto et al., 2013), with the constant of proportionality designated $\alpha_{\rm CO}$. Even in this idealized case of a homogeneous cloud, $\alpha_{\rm CO}$ depends on the precise conditions within the cloud. Draine (2011) derives the dependence of $\alpha_{\rm CO}$ on the excitation temperature $T_{\rm exc}$ and molecular gas density $n_{{\rm H}_{2}}$ as

[TABLE]

where we have adopted a factor of 1.36 to convert from hydrogen mass to total gas mass, which accounts for the abundance of He (Bolatto et al., 2013). We note that for a density $n_{\rm H_{2}}=2\times 10^{3}\,\mathrm{cm^{-3}}$, typical of GMCs, $T_{\rm exc}=10$ K implies a CO-to-$\rm H_{2}$ conversion factor of $\alpha_{\rm CO}\approx 4.4\,M_{\odot}(\mathrm{K\,km\,s^{-1}\,pc^{2}})^{-1}$, consistent with the value inferred from observations (Bolatto et al., 2013).

A population of virialized molecular clouds will likewise have a linear relationship between the total molecular gas mass and the integrated CO(1-0) line luminosity provided that the covering factor is low enough just that the CO emission from one cloud is unlikely to be absorbed by another cloud (Dickman et al., 1986; Bolatto et al., 2013).

Under these assumptions, we can write the CO luminosity directly in terms of the molecular gas mass $M_{{\rm H}_{2}}\equiv f_{{\rm H}_{2}}M_{\rm H}$:

[TABLE]

We treat $\alpha_{\rm CO}$ as a parameter to be fit. While there are indications that $\alpha_{\rm CO}$ may vary systematically with other galaxy properties, e.g., metallicity (Genzel et al., 2012; Bolatto et al., 2013; Sandstrom et al., 2013), we do not consider such variations here.

4 Intensity Mapping Framework

4.1 Modeling the Fluctuation Signals

In this section, we introduce a simple, generic halo occupation distribution (HOD) model, which is used to compute the power spectra that describe the spatial fluctuations of various signals emitted from discrete galaxies. Incorporating the correlation of subhalo structure (e.g., satellite galaxies) via such an HOD model is non-trivial, since both observational and theoretical studies have shown that massive dark matter halos tend to host more than one galaxy at low redshifts (e.g., Gao et al., 2011; McCracken et al., 2015), with a peak in subhalo abundance for a given halo mass at $z\sim 2-3$ as found by Wetzel et al. (2009). The original HOD model describes the occupation of halos by central and satellite galaxies (Kravtsov et al., 2004). Here, we generalize it to describe the fluctuations in line signals associated with the clustering of both central and satellite galaxies by weighting the galaxy number counts by a measure of the signal strength $S_{\nu}$ at observed frequency $\nu$ for a given halo mass and redshift, which means slightly differently for different signals (see later text). In particular, we define the number-count-weighted signal strengths of central and satellite galaxies,

[TABLE]

and

[TABLE]

where $N_{\mathrm{cen}}$ is the number of central galaxies in a halo, which is equal to 1 for $M>M_{\mathrm{min}}=10^{10}\,M_{\odot}$ and 0 otherwise (Wu & Doré, 2017), and $\mathrm{d}n/\mathrm{d}m$ represents the subhalo mass function, for which we adopt the fitting function in Tinker & Wetzel (2010). We consequently define the mean radiation strength as

[TABLE]

We note that in our expression, for different signals, $\bar{j}_{\nu}$ represents slightly different physical quantities and thus carries different units by convention. Specifically, $\bar{j}_{\nu}$ denotes the mean volume emissivity, intensity, and brightness temperature for the CIB,222This shows how the HOD formalism is originally defined in the CIB anisotropy model, provided here for completeness and better illustrating our generalization. [C ii]/[N ii]/CO lines and H i 21cm line, respectively. For the signals under consideration, we have

[TABLE]

[TABLE]

[TABLE]

where the units of signal strengths are $\rm erg\ s^{-1}\ Hz^{-1}sr^{-1}$, $\rm cm\ erg\ s^{-1}\ Hz^{-1}sr^{-1}$ and $\rm mK\ cm^{3}$, respectively. The mapping from frequency to line-of-sight distance is given by $y(z)=\mathrm{d}\chi/\mathrm{d}\nu=c(1+z)/\left[\nu H(z)\right]$, where $\chi$ denotes the comoving radial distance.

Generally, the power spectrum of a pair of signals at frequencies $\nu$ and $\nu^{\prime}$ (auto-correlation if $\nu=\nu^{\prime}$ and cross-correlation otherwise) can be expressed as the sum of one-halo, two-halo, and shot-noise components, namely

[TABLE]

The one-halo term characterizes the contribution to the fluctuations from emitters residing in the same halo. Assuming that the occupation statistics of central and satellite galaxies are independent and that the latter is Poissonian, we have

[TABLE]

where $u(k|M,z)$ is the normalized Fourier transform of the halo density profile (Navarro et al., 1997; Bhattacharya et al., 2013). The two-halo component describes the contribution from emitters residing in different halos,

[TABLE]

where $P_{\delta\delta}(k,z)$ is the dark matter power spectrum and

[TABLE]

with $b(M,z)$ being the halo bias factor (Tinker et al., 2008). Finally, the shot-noise component describes the self-correlation of emitters due to their discrete nature,

[TABLE]

which can be considered as the $k\rightarrow 0$ limit of the one-halo term in the absence of satellite galaxies (see, e.g., Wolz et al., 2017).

Finally, following Sun et al. (2018), in order to take into account of the stochasticity of individual galaxies, we introduce a simple parameterization of a lognormal distribution of line brightness below for a given halo mass and redshift. The probability density can be expressed as

[TABLE]

where $\mu_{\nu}=\log[S_{\nu}(M,z)]$ is the aforementioned mean line strength and $\sigma_{\nu}=0.3$ dex is our fiducial level of scatter reflecting the typical galaxy-to-galaxy variation in line production. It is straightforward to show that the power spectrum averaged over the lognormal distribution is essentially a rescaling of the power spectrum without scatter, specified by the additive correction factors in the following relations:

[TABLE]

which applies to the two-halo term of power spectrum scaling as the square of the first luminosity moments, and

[TABLE]

which applies to the one-halo and shot-noise terms of power spectrum scaling as the second luminosity moment. Figure 8 shows how the power spectrum is affected by the above correction factors in the presence of a non-trivial scatter $\sigma_{\nu}$. For comparison, the open squares indicate the average of power spectra directly drawn from 1000 random realizations of the lognormally distributed $S_{\nu}(M,z)$ relation, which agrees well with the one analytically derived using $\langle\mu_{\nu}\rangle$ and $\langle 2\mu_{\nu}\rangle$ as shown by the solid curve (a sum of rescaled one-halo, two-halo, and shot-noise components). In our power spectrum analysis, we include these correction factors to obtain constraints on the lognormal scatter together with physical properties of the ISM. We further assume, for simplicity, that similar physical processes (e.g., regulations of galaxy evolution by star formation, outflows and interactions, variations of stellar population and ISM properties) give rise to the stochasticity for a given halo mass and redshift, and therefore line luminosities considered in this work all share the same lognormal scatter $\sigma$.

4.2 Sensitivity Analyses

In this section, we describe the formalism to forecast the sensitivity to the power spectrum signal, assuming a given experimental setup. For a three-dimensional survey of volume $V_{\mathrm{s}}=L_{x}L_{y}L_{z}$, the observed 3D power spectrum $\mathcal{P}(K)$ for a given mode $K$ in the Fourier space of the observing frame is related to the true, spherically-averaged power spectrum $\Delta^{2}(k)=k^{3}P(k)/2\pi^{2}$ by

[TABLE]

where $W(k,K)$ is a convolution kernel commonly referred to as the “window function,” which is determined by the survey geometry. Here we only consider the simple situation that the survey volume is large enough such that $W(k,K)$ can be well-approximated by a function sharply peaking at $k\sim K$, which yields $\mathcal{P}(K)\approx P(k)$. Following S16, we write the uncertainty of the power spectrum $P(k)$ as the sum of a sample variance (i.e., cosmic variance) term and a thermal noise term. In particular, for the auto power spectrum $P_{\nu\nu}(k)$, we have

[TABLE]

where $G(k)$ denotes a smoothing factor due to finite spatial and spectral resolutions, which attenuates the power spectrum at large $k$ values beyond resolvable scales and is defined as Li et al. (2016)

[TABLE]

where $\mu=\cos\theta$ is the cosine of the angle a given $k$ vector makes with respect to the line of sight. For any given frequency channel width $\delta_{\nu}$, the spatial and spectral resolutions in physical units are given by $\sigma_{\parallel}(z)=k^{-1}_{\parallel,\mathrm{max}}(z)=y(z)\delta_{\nu}$ and $\sigma_{\perp}(z)=k^{-1}_{\perp,\mathrm{max}}(z)=\chi(z)\sqrt{\Omega_{\rm beam}}$, respectively. For the cross power spectrum $P_{\nu\nu^{\prime}}(k)$, we have

[TABLE]

where

[TABLE]

The (averaged) power spectrum of thermal noise is scale-independent and can be expressed as

[TABLE]

Using the radiometer equation, we can compute the on-sky sensitivity from the noise equivalent flux density (NEFD) or system temperature $T_{\rm sys}$, the beam size

[TABLE]

and the observing time per voxel

[TABLE]

as

[TABLE]

where $D_{\rm ap}$ and $N_{\rm feeds}$ represent the instrument’s effective aperture size and number of feeds (i.e., the number of spatial channels or spectrometers simultaneously on sky), respectively; the radio astronomy convention is adopted in the second equality. The voxel size can be derived from the spectral and angular resolutions as

[TABLE]

As long as the survey has proper spectral and angular resolutions to sample the $k$ space in a roughly isotropic manner, the number of (independent) modes $N_{\rm modes}$ can be calculated as (e.g., Furlanetto & Lidz 2007; Li et al. 2016; S16)

[TABLE]

where the factor of $1/2$ comes from the fact that the power spectrum is the Fourier transform of a real-valued function and thus only half of the Fourier plane contains independent information. The total signal-to-noise ratio (S/N) of a measured power spectrum is then defined to be (Gong et al., 2012; Li et al., 2016)

[TABLE]

The values of relevant instrumental parameters, adopted to guarantee significant detections of the LIM signals with comparable total S/N in our analysis, are summarized in Table 3 for each of the four case studies to be discussed in Section 6. We note that certain requirements presented exceed the scope of planned surveys and better resemble future mission concepts, for instance, a 10 m class, FIR telescope in space like the Origins Space Telescope (OST) for measuring [C ii] and [N ii] at intermediate redshifts, as well as the large number of feeds that will be enabled by the successful deployment of broadband, on-chip spectrometers like SuperSpec (Hailey-Dunsheath et al., 2014) and DESHIMA (Endo et al., 2019). Meanwhile, the detector noise levels assumed for some signals (e.g., [N ii]) are substantially more optimistic than what may be achieved from the ground, and therefore require observations in space, in which case an NEFD of order of $\rm 10\,mJy\,s^{1/2}$, corresponding to a noise equivalent power (NEP) of a few times $\rm 10^{-19}\,W\,Hz^{-1/2}$, is achievable (Bradford et al., 2008, 2018). Even though we make no effort to carefully build these case studies on existing or planned experiments, auto-/cross-correlation opportunities based on real experiments in similar contexts will be described.

5 Comparison to Existing Observational Constraints

From LIM observations of the large-scale distribution of neutral hydrogen, constraints have been placed on the H i density parameter, defined as the ratio of the H i density to the critical density of the universe at $z=0$, namely, $\Omega_{{\mathrm{H\,\mathsc{i}}}}=\rho_{{\mathrm{H\,\mathsc{i}}}}/\rho_{\rm c,0}$ or equivalently the mean H i brightness temperature $\bar{T}_{{\mathrm{H\,\mathsc{i}}}}$ as defined in Equation 22. The top two panels of Figure 9 show the product of H i density parameter and bias factor, degenerate when constrained by the large-scale clustering of H i, and the mean 21cm brightness temperature predicted by our model, respectively, which are found to be in good agreement with observed values at $z\sim 0.8$ (Chang et al., 2010; Switzer et al., 2013). The corresponding H i power spectra $\Delta^{2}_{{\mathrm{H\,\mathsc{i}}}}$ derived from our HOD model at $z=0,0.8$ and 1 are shown in the bottom panel of Figure 9, together with the deep-field results (detections only) from Switzer et al. (2013). While the detections shall be interpreted as upper limits since residual, correlated foregrounds are very likely present, predictions by our reference ISM model are still broadly consistent with H i observations available to date.

The top panel of Figure 10 shows the $L_{[{\mathrm{C\,\mathsc{ii}}}]}$–SFR relations derived from our model assuming different photoelectric heating efficiency (from bottom to top, $\epsilon_{\rm PE}=3\times 10^{-4},1\times 10^{-3},3\times 10^{-3},1\times 10^{-2}$ and $3\times 10^{-2}$) and how they compare with the best-fit relation to a large sample of galaxies of various populations (starburst galaxies, dwarfs, ULIRGs, AGNs, high-$z$ galaxies, etc.) taken from De Looze et al. (2014). Recently, Pullen et al. (2018) and Yang et al. (2019) report a tentative detection of excess emission in the 545 GHz Planck map that can be attributed to redshifted [C ii] line emission. From angular cross-power spectra of high-frequency Planck maps with BOSS quasars and CMASS galaxies, a joint constraint on the product of mean [C ii] intensity and bias factor $b_{[{\mathrm{C\,\mathsc{ii}}}]}I_{[{\mathrm{C\,\mathsc{ii}}}]}=2.0^{+1.2}_{-1.1}\times 10^{5}\,\mathrm{Jy\,sr^{-1}}$ is inferred at 95% confidence level. In the bottom panel of Figure 10, we compare our model predictions at the five different $\epsilon_{\rm PE}$ values against the measurement from Yang et al. (2019). We note that a relatively high $\epsilon_{\rm PE}$ is required to match the measured level of $b_{[{\mathrm{C\,\mathsc{ii}}}]}I_{[{\mathrm{C\,\mathsc{ii}}}]}$, which may lead to tension with the observed $L_{[{\mathrm{C\,\mathsc{ii}}}]}$–SFR relation. Such a discrepancy is also observed by Pullen et al. (2018) and Yang et al. (2019) when comparing against phenomenological models (e.g., Gong et al., 2012; Silva et al., 2015) based on local observations. While it is possible that the $L_{\rm IR}$–$L_{[{\mathrm{C\,\mathsc{ii}}}]}$ relation is different at these redshifts, in which case a deviation from the proportionality $L_{[{\mathrm{C\,\mathsc{ii}}}]}\propto\mathrm{SFR}$ may be implied (see, e.g., the data-driven model of [C ii] emission presented by Padmanabhan 2019), the observed excess may also be produced by non-[C ii] factors such as interloper lines or redshift evolution of CIB parameters. Future, high-resolution [C ii] LIM surveys will help clarify this discrepancy.

Measuring CO power spectrum from dedicated LIM experiments, such as COPSS II (Keating et al., 2016), COMAP (Li et al., 2016) and Y. T. Lee Array (Ho et al., 2009), or galaxy surveys (e.g., Uzgil et al., 2019) is an emerging field. In Figure 11, we show our model predictions of the CO(1-0) power spectrum at $z=1$, evaluated for three pairs of excitation temperature $T_{\rm exc}$ and molecular gas density $n_{\rm H_{2}}$ to illustrate how sensitive CO power spectrum is to these gas properties. For comparison, the best estimate from an empirical model fit to a compilation of existing observations, including constraints on CO luminosity function and power spectrum obtained at redshifts $0<z<3$, taken from Padmanabhan (2018) is shown by the shaded band. The prediction of our reference ISM model is in good agreement with the observational constraints.

As there has not been any LIM measurement of [N ii] lines because of their faintness, in Figure 12 we only compare our reference $L_{[{\mathrm{N\,\mathsc{ii}}}]}$-$L_{\rm IR}$ model against results from the literature and then present the [N ii] power spectra it predicts. The top panel of Figure 12 shows the relations between [N ii] line luminosities and the star formation rate predicted by our reference ISM model assuming $T_{{\mathrm{H\,\mathsc{ii}}}}=10^{4}\,$K and $n_{\rm e,{\mathrm{H\,\mathsc{ii}}}}=10^{2}\,\mathrm{cm^{-3}}$. Estimates from previous work are shown for comparison, including scaling relations (Spinoglio et al. 2012; S16)333The scaling relation for [N ii] 205$\mu$m is not provided by Spinoglio et al. (2012), for which we assume a line ratio of $L^{122}_{[{\mathrm{N\,\mathsc{ii}}}]}/L^{205}_{[{\mathrm{N\,\mathsc{ii}}}]}=3$ following S16 (corresponding to $n_{\rm e,{\mathrm{H\,\mathsc{ii}}}}\sim 100\,\mathrm{cm^{-3}}$, as can be seen from Figure 7). based on a sample of local galaxies observed with the ISO-LWS spectrometer (Clegg et al., 1996) and compiled by Brauher et al. (2008), and relations derived by Herrera-Camus et al. (2016) based on an observationally-motivated prescription assuming a uniform $n_{\rm e,{\mathrm{H\,\mathsc{ii}}}}=10^{2}\,\mathrm{cm^{-3}}$. Given the relatively large dispersion that exists in the existing data (see, e.g., Spinoglio et al., 2012), our simple model is deemed satisfactory despite the fact that it may slightly overestimate the local [N ii] luminosities. The bottom panel of Figure 12 shows the power spectra of [N ii] 122 $\mu$m and 205 $\mu$m lines, evaluated at $z=2$ for two different values of the H ii region electron number density to illustrate the density effect on the [N ii] line ratio.

6 Inferring ISM Properties from Auto/Cross-Correlations

Both auto-correlation and cross-correlation analyses serve as a powerful tool to study the ISM physics when LIM data sets of multiple lines are available. The latter, however, has the advantage of avoiding contamination from uncorrelated foregrounds (line and continuum), which are usually a few orders of magnitude brighter than and/or spectrally blended with the signal of interest, therefore presenting a great challenge to reliably measuring the line intensity fluctuations (Lidz et al. 2009; Pullen et al. 2013; Silva et al. 2015; S16; Beane et al. 2019; see also Switzer et al. 2019). In the rest of this section, we present several case studies in order to demonstrate how the population-averaged physical properties of different ISM phases, such as their gas temperature and density, might be reliably extracted by auto/cross-correlating the intensity fields of different tracers.

We adopt a Bayesian analysis framework and fit parameters of ISM properties with the affine-invariant Markov Chain Monte Carlo (MCMC) code emcee (Foreman-Mackey et al., 2013). The likelihood function for fitting the mock power spectra can be expressed as

[TABLE]

where $N_{k}$ is the number of $k$ bins in which auto or cross power spectra are measured and $N_{s}$ represents the number of auto/cross-correlation surveys being included. The probability of the data vector $\hat{x}$ is described by a normal distribution

[TABLE]

where $\sigma_{ij}$ represents the gaussian error associated with the measurement. Broad, uninformative priors on the model parameters $\hat{\theta}$ are used, whose values are to be stated below for each individual case study. Furthermore, for all the following case studies, we adopt the same range and binning scheme for $k$ which yield 15 bins evenly-spaced in $\log k$ over $-1.5<\log[k/(h/\mathrm{Mpc})]<1$. We stress that while all four case studies presented below are evaluated at $z\sim 2$ for a redshift interval of $\Delta z=\pm 0.25$, the same exercise could be repeated at different redshifts in order to study the redshift evolution of different ISM properties, which is one of the most important applications of the modeling framework presented.

It is also important to point out that our model implicitly enforces a linear relation between line luminosity and halo mass, which is likely an oversimplification given the complicated physics involved in line production. The only physical parameter that modifies the shape of line luminosity function (and therefore the shape of power spectrum) is the log scatter $\sigma_{\nu}$. Consequently, whether the constraining power comes from the clustering or shot-noise regime of the power spectrum only makes a moderate difference in our analysis, as will be shown in Case I. The scale dependence of constraining power in a power spectrum analysis without such simplification can be found in recent studies (e.g., Yue & Ferrara, 2019).

6.1 Case I: Multi-Phase Diagnosis with H i and CO

As the first example, we investigate how the multi-phase ISM may be probed by a combination of H i and CO LIM observations, which trace atomic and molecular hydrogen, respectively. Because the total gas mass is constrained implicitly by the CIB, the H i measurement constrains both the atomic and molecular gas fraction. The CO measurement can then in principle break the degeneracy between the total amount of molecular gas and $\alpha_{\rm CO}$.

We consider two independent, mock measurements of H i and CO auto power spectra, generated at $z\sim 2$ assuming the reference ISM model described in Table 2 and experimental setups specified in Table 3, which yield a total S/N of approximately 8 for each signal444Summed over all $k$ bins; see Equation (61). Forthcoming single-dish/interferometric suveys, including FAST (Bigot-Sazy et al., 2016), CHIME (Bandura et al., 2014) and SKA (Dewdney et al., 2013) for H i and COMAP (Li et al., 2016; Chung et al., 2019) and mmIME (Keating et al. in preparation) for CO, will carry out these auto-correlation measurements directly, even though in both cases the signal is expected to be heavily contaminated by line/continuum foregrounds. Broad, flat priors over $0<f_{\rm H_{2}}<0.5$, $0<f_{\rm{\mathrm{H\,\mathsc{ii}}}}<0.5$, $10^{1}<n_{\rm H_{2}}/\mathrm{cm^{-3}}<10^{5}$ and $0<\sigma<1$ are assumed for the MCMC analysis. The MCMC sampling is constructed with 60 walkers, 500 burn-in steps—well above the estimated autocorrelation time ($\sim 50$ steps) emcee returned, and another 500 steps for sampling.

The left panel of Figure 13 shows mock observed power spectra of H i 21cm and CO auto-correlation signals at $z\sim 2$, where the error bars are calculated from the assumed instrument parameters. The joint and marginalized posterior distributions of free parameters constrained by the mock auto power spectra under the MCMC framework are shown in the right panel of Figure 13. Note that we have converted the posterior of $n_{\rm H_{2}}$ into the more commonly seen $\alpha_{\rm CO}$ factor using the assumed $T_{\rm exc}=10\,$K.

From the comparison between posterior distributions and true values (orange plus signs), as well as the fact that none of them are prior dominated, constraining power on all four parameters is observed. Even $f_{\rm H_{2}}$ and $\alpha_{\rm CO}$, though still strongly correlated, are individually constrained in this analysis. However, more precise estimation of the molecular gas content of galaxies from CO power spectrum measurements is conditional on how well $\alpha_{\rm CO}$ can be reliably determined, even if additional information about the atomic hydrogen content from H i LIM is available. In practice, the exact value of $\alpha_{\rm CO}$ could vary in a non-trivial way with physical conditions of molecular gas in galaxies, especially the gas temperature distribution and metallicity. As a result, how LIM might be exploited to better determine its value is an interesting topic to be explored (see Section 7 for further discussion).

In addition to the default scenario using the full power spectrum in all $k$ bins, we consider an alternative scenario, where only the shot-noise power can be measured, while holding the overall S/N fixed. This resembles deep, targeted observations by, e.g., ALMA, from which information about large-scale intensity fluctuations is not available. As indicated by the gray error bars in the left panel of Figure 13, we assume two $\rm S/N\sim 8$ measurements of H i and CO power spectrum at $k\approx 1\,h/\mathrm{Mpc}$ where shot noise is dominant. Due to the implicitly assumed linearity between line luminosity and halo mass, similar constraining power on the parameter space is observed, except that the measured log scatter $\sigma_{\nu}$ becomes biased, which can be easily understood given that in our model it is the only parameter sensitive to the shape of the power spectrum.

6.2 Case II: Multi-Phase Diagnosis with H i, [C ii] and CO

Given the observed degeneracy between $f_{\rm H_{2}}$ and $\alpha_{\rm CO}$ in the previous case study, which introduces ambiguity to the interpretation of CO LIM results in terms of a molecular gas census, we investigate in this case how the inclusion of [C ii] data, an indirect tracer of the molecular hydrogen fraction as indicated by Equation (24), may help alleviate such a degeneracy. Additionally, we investigate how the constraining power on the parameter space may differ between using the three separate auto power spectra and using the $3(3-1)/2=3$ cross-correlation measurements available, which has the advantage of being immune to contamination from uncorrelated foregrounds as suggested in S16.

Mock data sets of LIM observations are again created assuming the reference ISM model parameters and instrument parameters listed in Table 2 and Table 3, respectively. We note that when accounting for the effect of finite beam size in the cross-correlation sensitivity analysis, we conservatively evaluate for the coarser beam throughout our calculations. The overall S/N of cross-correlation data ($\rm SNR_{tot}\sim 5$ for each cross signal) is consequently lower than that of auto-correlation data. At intermediate redshifts, experiments like EXCLAIM (Switzer, 2017) and TIM (Aguirre & STARFIRE Collaboration, 2018) will measure [C ii] in tomography, which, when spatially overlapped, may be combined with the H i and CO surveys mentioned in the previous case to obtain their mutual cross-correlations. Broad, flat priors over $0<f_{\rm H_{2}}<0.5$, $0<f_{\rm{\mathrm{H\,\mathsc{ii}}}}<0.5$, $10^{-4}<\epsilon_{\rm PE}<10^{-1}$, $10^{1}<n_{\rm H_{2}}/\mathrm{cm^{-3}}<10^{5}$ and $0<\sigma<1$ are assumed for the MCMC analysis. The MCMC sampling is constructed with 50 walkers, 500 burn-in steps—sufficiently larger than the estimated autocorrelation time ($\sim 60$ steps)—and another 500 steps for sampling.

In the left panel of Figure 14, we show the mock power spectrum data sets in addition to what has been shown in Figure 13, including auto-correlation of [C ii] and mutual cross-correlations of the three lines considered, all evaluated at $z\sim 2$. The corresponding constraining power on the parameter space of our mock auto-correlation (brown dashed contours) and cross-correlation (gray dashed-dotted contours) and auto/cross-correlation-combined (black solid contours) data sets is presented in the right panel of Figure 14 as joint and marginalized posterior distributions.

From the posterior constrained by auto-correlations, which becomes less biased from the true value after including [C ii] data, it is clear that the degeneracy between $f_{\rm H_{2}}$ and $\alpha_{\rm CO}$ has been substantially reduced, although considerable uncertainty is still associated with $f_{\rm H_{2}}$. Other parameters, including $\alpha_{\rm CO}$, $\epsilon_{\rm PE}$ and $\sigma$, are well constrained by the auto-correlations from their marginalized posteriors, except for $f_{{\mathrm{H\,\mathsc{ii}}}}$ which is not directly traced by any of the lines. The constraining power from cross-correlations, on the other hand, is not as good—particularly for $f_{{\mathrm{H\,\mathsc{ii}}}}$ of which the constraint is prior dominated—yet still significant in general. While formally when $N\geq 3$, perfectly correlated lines are present, the mean line intensities shall be constrained equally well by their mutual cross-correlations; the poorer performance can be largely explained by the lower overall S/N of cross-correlation data. For completeness, we show also the total constraining power combining both auto-correlation and cross-correlation data sets, even though it is only slightly improved compared with the auto-only case. Both kinds of measurement are subject to realistic but different limitations — while analysis based on auto correlations tends to be more foreground-contaminated in general, it has the advantage of not requiring the experiments to be spatially overlapped, as long as the cosmic variance of individual surveys can be properly accounted for.

6.3 Case III: Probing H ii Regions with [N ii] Lines.

Another straightforward application of our line model is to use the two [N ii] lines to constrain the state of ionized ISM, especially its electron number density $n_{\rm e,{\mathrm{H\,\mathsc{ii}}}}$ directly probed by the [N ii] fine-structure line ratio (see, e.g., Goldsmith et al. 2015 and Díaz-Santos et al. 2017 for applications of the $[{\mathrm{N\,\mathsc{ii}}}]_{205}/[{\mathrm{N\,\mathsc{ii}}}]_{122}$ ratio as a diagnostic of $n_{\rm e,{\mathrm{H\,\mathsc{ii}}}}$ to the Galactic plane and local galaxies). Here, we consider two types of measurements, namely, the cross power spectrum of the two [N ii] lines and their respective auto power spectra.

Mock data sets of LIM observations are created assuming the reference ISM model ($n_{\rm e,{\mathrm{H\,\mathsc{ii}}}}=100\,\mathrm{cm^{-3}}$, $T_{\rm gas,{\mathrm{H\,\mathsc{ii}}}}=10^{4}$ K and $\sigma=0.3$ dex), together with experimental setups specified in Table 3. We note that while [N ii] lines tend to be spectrally covered by [C ii]-targeted experiments like EXCLAIM and TIM, at intermediate redshifts the required specifications in this (and the next) case study for a detection may only be achievable for next-generation space missions like OST owing to the faintness of [N ii] emission. Broad, flat priors over $1<n_{\rm e,{\mathrm{H\,\mathsc{ii}}}}/\mathrm{cm^{-3}}<10^{5}$, $10^{3}<T_{\rm gas,{\mathrm{H\,\mathsc{ii}}}}/\mathrm{K}<10^{5}$ and $0<\sigma<1$ are assumed for the MCMC analysis. The MCMC sampling in either case is done with 100 walkers, 1000 burn-in steps—well above the estimated autocorrelation time ($\sim 100$ steps)—and another 1000 steps for sampling.

Figure 15 shows the posterior distributions of the electron number density $n_{\rm e,{\mathrm{H\,\mathsc{ii}}}}$, the gas temperature $T_{\rm gas,{\mathrm{H\,\mathsc{ii}}}}$ and the lognormal scatter in line intensity $\sigma$, as constrained by the two types of observations, respectively. With the assumed model and survey parameters, the auto and cross power spectra $P^{[{\mathrm{N\,\mathsc{ii}}}]}_{122}$, $P^{[{\mathrm{N\,\mathsc{ii}}}]}_{205}$ and $P^{[{\mathrm{N\,\mathsc{ii}}}]}_{122\times 205}$ are measured at a total S/N of $\mathrm{SNR_{tot}}\sim 5.2$, 4.7 and 6.7, respectively. Both methods are able to determine the density and temperature without significant bias. Nevertheless, the auto power spectra of both [N ii] lines together are much more effective than the just the cross power at breaking the degeneracy between the density and temperature as probed by the line ratio (see Figure 7). We note that the difference in the constraining power between these two contrasting cases serves as an example of the importance of determining the amplitude of each individual tracer, through measurements of either the individual auto power spectra, or all mutual cross power spectra when $N\geq 3$ lines are detected, as suggested in S16 and demonstrated in the previous case study.

6.4 Case IV: Dissecting [C ii] Origin with [C ii]–[N ii] 205$\mu$m Line Ratio

While in Section 3.2 we have assumed that the [C ii] line is solely attributed to the atomic gas in the PDRs so as to keep the line model simple, a small yet non-trivial fraction of the observed [C ii] emission may actually originate in ionized gas phases as suggested by several recent studies (Hughes et al., 2015; Croxall et al., 2017; Cormier et al., 2019). Therefore, in this final example we consider a slight extension of the [C ii] model presented: we rewrite the total [C ii] emission observed with LIM as $L^{\rm tot}_{[{\mathrm{C\,\mathsc{ii}}}]}=L_{[{\mathrm{C\,\mathsc{ii}}}]}/f_{[{\mathrm{C\,\mathsc{ii}}}]}^{\rm neutral}$, where $L_{[{\mathrm{C\,\mathsc{ii}}}]}=(1-f_{\rm H_{2}})\epsilon_{\rm PE}L_{\rm IR}$ is the contribution from the neutral ISM (PDRs) defined in Section 3.2 and $f_{[{\mathrm{C\,\mathsc{ii}}}]}^{\rm neutral}$ is an extra parameter introduced here to describe the fraction of $[{\mathrm{C\,\mathsc{ii}}}]$ emission contributed by the neutral ISM. Following Croxall et al. (2017), we use the ratio of [C ii]/[N ii] 205 $\mu$m lines, whose critical densities for electron collisions are very similar ($n^{\rm crit}_{[{\mathrm{C\,\mathsc{ii}}}]}\sim 45\,\mathrm{cm^{-3}}$ and $n^{\rm crit}_{[{\mathrm{N\,\mathsc{ii}}}],205}\sim 32\,\mathrm{cm^{-3}}$), as a diagnostic of $f_{[{\mathrm{C\,\mathsc{ii}}}]}^{\rm neutral}=1-R_{\rm ionized}L_{[{\mathrm{N\,\mathsc{ii}}}],205}/L_{[{\mathrm{C\,\mathsc{ii}}}]}$, where $R_{\rm ionized}\approx 4$ denotes the ionized gas [C ii]/[N ii] ratio implied by their respective collision rates with electrons (Blum & Pradhan, 1992; Tayal, 2008, 2011), assuming Galactic gas-phase abundances.

Figure 16 demonstrates the constraining power on the neutral-phase contribution $f_{[{\mathrm{C\,\mathsc{ii}}}]}^{\rm neutral}$ to the observed [C ii] emission, estimated from the line ratio $L_{[{\mathrm{N\,\mathsc{ii}}}],205}/L_{[{\mathrm{C\,\mathsc{ii}}}]}$ inferred from mock LIM observations at $z\sim 2$ (left panel). For the MCMC analysis, we assume a gaussian prior $\mathcal{N}(4,0.4)$ for $R_{\rm ionized}$, whereas a broad, flat prior is used for $f_{[{\mathrm{C\,\mathsc{ii}}}]}^{\rm neutral}$. The sampling is done with 100 walkers, 500 burn-in steps—sufficiently large compared with the estimated autocorrelation time ($\sim 30$ steps)—and another 500 steps for sampling.

Our reference ISM model assumes a high $f_{[{\mathrm{C\,\mathsc{ii}}}]}^{\rm neutral}\sim 0.9$ (orange plus sign), consistent with the finding that [C ii] emission arises mostly from the neutral ISM. For the survey specifications given in Table 3, we find that a population-averaged, neutral-phase contribution $f_{[{\mathrm{C\,\mathsc{ii}}}]}^{\rm neutral}$ can be robustly determined by simultaneously observing [C ii] and [N ii] 205 $\mu$m lines with LIM. We note that, in principle, the simple diagnostic described may be subject to density effects on the [C ii]/[N ii] line ratio and $R_{\rm ionized}$, although in both cases the dependence on $n_{\rm e}$ is found to be weak (Croxall et al., 2017).

7 Discussion and Conclusion

We have presented a simple analytical framework to self-consistently model the production of emission lines in the multi-phase ISM, based on the mean dark matter halo properties derived from a model fit to the observed CIB anisotropy. The redshift evolution of cosmic star formation, dust mass, gas (total and molecular) mass, gas-phase metallicity, and the strengths of H i, [C ii], [N ii] and CO lines predicted by our model have been compared with observations, showing that our model, despite its simplicity, can describe the production of lines in the ISM in a physically-motivated way. We have illustrated how this modeling framework can be used to reconstruct the average properties of different ISM phases, such as the mass fractions and densities of neutral and ionized gas, the photoelectric heating efficiency in the PDRs, and so forth, over a wide range of redshifts from multi-tracer LIM observations.

Our analysis underscores the importance of cross-correlation analyses. While equivalent information may be obtained from the auto-correlation of respective tracers, cross-correlation analysis in the same cosmological volume is minimally susceptible to foreground contamination. With the large number of upcoming LIM experiments targeting lines produced in different ISM phases, e.g., CCAT-Prime (Stacey et al., 2018), CHIME (Bandura et al., 2014), COMAP (Li et al., 2016), CONCERTO (Lagache, 2018), HIRAX (Newburgh et al., 2016), SKA (Santos et al., 2015), SPHEREx (Doré et al., 2014), Tianlai (Xu et al., 2015), TIM (Aguirre & STARFIRE Collaboration, 2018) and TIME (Crites et al., 2014), our understanding of the ISM evolution and physical processes dominating line emission over cosmic time is expected to be greatly deepened by the coarse-grained view built up from LIM surveys with multiple tracers.

The simplicity and modularity of the model presented here lends itself to straightforward improvements and extensions to incorporate more sophisticated treatments of both galaxy evolution and ISM physics motivated by observational and theoretical studies. For instance, to more reliably apply this framework to galaxies at higher redshifts, including the reionization era, it would be valuable to introduce additional calibrations and constraints from data sets at other wavelengths. Currently the star formation history is anchored only to FIR emission constrained by the CIB anisotropy, which is sensitive mostly to galaxies at redshift $z\lesssim 3$ (e.g., Viero et al., 2013b). However, much of the information about the galaxy–halo connection, feedback-regulated galaxy and ISM evolution, and so forth, is encoded in data at shorter wavelengths, e.g., the galaxy UV luminosity function (UVLF). We therefore expect the exact mass and redshift dependence of mean halo properties (see Section 2) to be better constrained out to the epoch of reionization by combining IR and UV data, which will be explored in future work.

Given the necessarily coarse-grained picture of galactic ISM properties painted by LIM, we have employed physically-motivated but ultimately simple prescriptions for the line emission physics. Particularly as new observations yield more model constraints, the line physics can be refined. Notably, our prescription for the [C ii] emission does not account for the deficit relative to $L_{\rm IR}$ observed in luminous and ultraluminous galaxies (e.g., Malhotra et al., 1997). It should be explored whether this can be recovered within the modeling framework by introducing the effect of dust charging on $\epsilon_{\rm PE}$ (Bakes & Tielens, 1994) and the saturation of [C ii] at high gas temperatures (Muñoz & Oh, 2016). Further, given that this effect appears to be a strong function of galaxy luminosity, this may have important testable implications for the predicted [C ii] power spectra.

Another important caveat to our simple prescription lies in the interpretation of LIM signals in terms of globally-averaged ISM properties. In reality, the ensemble of gas clouds within a galaxy, while all contributing to the same line emission, may have a wide distribution of physical properties (e.g., $\rm H_{2}$ gas temperature). Likewise, these distributions may vary significantly among different galaxy populations. As a result, interpreting LIM data in terms of a single “mean” property is an oversimplification. A more robust extraction of ISM physics from LIM data sets could be achieved by modeling these distributions directly, perhaps incorporating prior information about ISM conditions that are known to vary systematically in different galaxy populations, as well as how the line production is coupled to these distributions through radiative processes. Such modeling is beyond the scope of this paper and will be the subject of future investigation.

Finally, in this work we selected a small subset of available lines to illustrate the power of LIM to probe the multi-phase ISM. However, the model is readily extensible to other lines. For instance, the [O i] 63 $\mu$m and [O iii] 88 $\mu$m lines are also important cooling lines, and so the sum of emission from these lines and [C ii] may yield a more robust correlation with $L_{\rm IR}$ (De Looze et al., 2014), with their relative importance as a function of redshift and galaxy properties providing constraints on the physical state of the emitting gas. Simultaneous measurements of multiple CO rotational lines are both a powerful probe of the physics of molecular gas and a means of validation since the CO lines should be spatially correlated. In addition to CO, H2 rotational lines can be used to constrain the molecular gas content of galaxies, meanwhile shedding light on the gas temperature distribution (Pereira-Santaella et al., 2014; Togi & Smith, 2016). Optical and UV lines of hydrogen such as Ly$\alpha$, H$\alpha$ and H$\beta$ are also being actively pursued by LIM experiments and should be incorporated, particularly given their potential to probe metal-poor environments in the very high redshift universe.

We would like to thank the anonymous referee for comments that helped improve this paper. We would like to thank Hao-Yi (Heidi) Wu for helpful discussion on the CIB model, as well as Garrett (Karto) Keating and Ryan Keenan for compiling and sharing the constraints on cosmic molecular gas content. We are also grateful to Jamie Bock, Matt Bradford, Patrick Breysse, Paul Goldsmith, Adam Lidz, Lunjun Liu, Lluis Mas-Ribas and Anthony Pullen for constructive discussion and comments on this work. Part of the research described in this paper was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. Software: corner Foreman-Mackey (2016), emcee Foreman-Mackey et al. (2013), hmf Murray et al. (2013), matplotlib Hunter (2007), numpy Walt et al. (2011) and scipy Jones et al. (2001).