{'fraction_non_alphanumeric': 0.040660284475354545, 'fraction_numerical': 0.025996606405932505, 'mean_word_length': 4.568412457716085, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dissecting the high-z interstellar medium through intensity mapping cross-correlations

We explore the detection, with upcoming spectroscopic surveys, of three-dimensional power spectra of emission line fluctuations produced in different phases of the Interstellar Medium (ISM) by ionized carbon, ionized nitrogen and neutral oxygen at redshift z>4. The emission line [CII] from ionized carbon at 157.7 micron, and multiple emission lines from carbon monoxide, are the main targets of planned ground-based surveys, and an important foreground for future space-based surveys like the Primordial Inflation Explorer (PIXIE). However, the oxygen [OI] (145.5 micron) line, and the nitrogen [NII] (121.9 micron and 205.2 micron) lines, might be detected in correlation with [CII] with reasonable signal-to-noise ratio (SNR). These lines are important coolants of both the neutral and the ionized medium, and probe multiple phases of the ISM. We compute predictions of the three-dimensional power spectra for two surveys designed to target the [CII] line, showing that they have the required sensitivity to detect cross-power spectra with the [OI] line, and the [NII] lines with sufficient SNR. The importance of cross-correlating multiple lines is twofold. On the one hand, we will have multiple probes of the different phases of the ISM, which is key to understand the interplay between energetic sources, the gas and dust at high redshift. This kind of studies will be useful for a next-generation space observatory such as the NASA Far-IR Surveyor. On the other end, emission lines from external galaxies are an important foreground when measuring spectral distortions of the Cosmic Microwave Background spectrum with future space-based experiments like PIXIE; measuring fluctuations in the intensity mapping regime will help constraining the mean amplitude of these lines, and will allow us to better handle this important foreground.


INTRODUCTION
Intensity mapping, introduced in Madau et al. (1997); Suginohara et al. (1999); Shaver et al. (1999), is an observational technique for measuring brightness fluctuations of emission lines produced by sources below the detection limit. Atomic and molecular emission lines, produced at a given redshift, are observed as fluctuations redshifted at a certain frequency, enabling us to map the three-dimensional structure of the Universe and compute, for each redshift slice, statistical quantities of interest such as the power spectrum. Intensity mapping, by measuring the aggregate radiation emitted by all galaxies in a given redshift slice, does not suffer from the incompleteness problem, while traditional galaxy surveys, being fluxlimited, do not detect the faintest galaxies. This can be a serious disadvantage if the galaxy luminosity function has a sufficiently steep end, as shown in Uzgil et al. (2014). One of the first and main targets of intensity mapping is the 21 cm neutral hydrogen line (Battye et al. 2004;Chang et al. 2010;Bull et al. 2015), which, in principle, opens a new window on both the formation of structures at high redshift and the history of reionization (Furlanetto et al. 2006). However, lines from other atoms and molecules can be used to constrain the physics of the ISM in a broad redshift range. The carbon [CII] fine-structure line at 157.7 µm, arising from the 2 P 3/2 → 2 P 1/2 fine-structure transition, is one of the most promising lines not only to understand star-formation in galaxies (Boselli et al. 2002;De Looze et al. 2011Herrera-Camus et al. 2015), but also to constrain the epoch of reionization and the physics of the ISM (Gong et al. 2011(Gong et al. , 2012Uzgil et al. 2014;Silva et al. 2015;Lidz & Taylor 2016;Cheng et al. 2016). Both theory and observations indicate that the atomic [CII] fine-structure is the dominant coolant of the neutral ISM (Hollenbach & Tielens 1999;Bernard-Salas et al. 2012), and one of the brightest lines in the Spectral Energy Distribution (SED) of a typical star-forming galaxy, with luminosities ranging from 0.01% to 1% of the total infrared luminosity (Stacey et al. 1991;Maiolino et al. 2005;Iono et al. 2006;Maiolino et al. 2009;Stacey et al. 2010;Ivison et al. 2010;Wagg et al. 2010;De Breuck et al. 2011). In fact, carbon is the forth most abundant element in the Universe. It has a low ionization potential, only 11.26 eV (see Table 1), below the 13.6 eV of hydrogen ionization; this ensures it is present both in the ionized and in the neutral medium. Moreover, the [CII] fine-structure transition of ionized carbon is characterized by a low temperature (91 K), and low critical density for collisions with hydrogen 4 . Intensity mapping from the rotational transitions of carbon monoxide, and, in particular, the lowest order transition CO(1-0) at 115 GHz, have also received increased attention in the past few years. Carbon monoxide emission lines at a given redshift act as a foreground contamination both for CMB observations (Righi et al. 2008;De Zotti et al. 2016), and for [CII] intensity mapping surveys targeting background galaxies at higher redshifts (Gong et al. 2012;Lidz & Taylor 2016;Cheng et al. 2016). Carbon monoxide molecules are easily produced from carbon and oxygen in star-forming regions, and CO intensity mapping provides information on the spatial distribution and redshift evolution of star formation in the Universe (Visbal & Loeb 2010;Carilli 2011;Lidz et al. 2011;Gong et al. 2011;Pullen et al. 2013;Breysse et al. 2014). At far-infrared (FIR) frequencies, many other lines can in principle be targeted by intensity mapping surveys, such as [OI] (63 µm and 145 µm), [NII] (122 µm and 205 µm), [OIII] (52 µm and 88 µm), and [CI] (610 µm and 371 µm), while proposed lines in other frequency bands include measurements of HeII (0.164 µm) to constrain properties of Population III stars (Visbal et al. 2015), Lyα (0.1216 µm) to probe reionization and star formation (Pullen et al. 2014), and OII (0.3737 µm) and Hα (0.6563 µm) to study the large scale clustering at redshifts 1 < z < 4 (Fonseca et al. 2016). As emphasized in Lidz & Taylor (2016), the sensitivity of intensity mapping measurements will rapidly increase in the near future, thanks to advances in detector technology, and some surveys are already in progress, or have been planned, to perform intensity mapping of one or more emission lines from sources at multiple redshifts. The CO Power Spectrum Survey (COPPS) (Keating et al. 2015) recently published measurement of the CO abundance and power spectrum from the CO(1-0) transition in the redshift range 2.3 < z < 3.3 (Keating et al. 2016), and the Carbon Monoxide Mapping Array Pathfinder (COMAP see Li et al. (2016)) has been proposed to study the CO emission at similar redshifts. Experiments targeting the [CII] emission line include the Tomographic Ionized-Carbon Mapping Experiment (TIME-Pilot, Crites et al. (2014)), and CONCERTO (CarbON CII line in post-rEionization and ReionizaTiOn epoch, Lagache et al., in preparation), while the Spectrophotometer for the History of the Universe, Epoch of Reionization, and Ice Explorer (SPHEREx) will focus on Lyα, Lyβ and [OIII] . The Cryogenic-Aperture Large Infrared-Submillimeter Telescope Observatory (CALISTO) (Bradford et al. 2015) has been proposed to measure, among other things, multiple FIR fine-structure transitions such as [NeII], [OI], [OIII] and, for z < 2, [CII]. Foregrounds are an important concern for intensity mapping surveys. Apart from the continuum emission from our Galaxy, a survey targeting an emission line observed at a given frequency ν obs will also detect the sum of emissions of N atoms or molecules α i coming from redshifts z i , whose lines are redshifted to the same observed frequency, so that the measured intensity I ν obs can be written as: (1) Different methods to overcome this difficulty have been proposed so far. Some authors (Visbal et al. 2011;Breysse et al. 2015;Silva et al. 2015) explore the possibility of mitigating this contamination by progressively masking the brightest pixels in the observed map. However, when dealing with [CII] maps at very high redshift (e.g. z = 7), a percentage of the signal will be masked in the process, and such a loss of information translates in a underestimation of the amplitude of the measured power spectrum (Breysse et al. 2015). This is unfortunate because, while the cosmological information content of the measured power spectrum is mainly encoded in its shape (primordial non-Gaussianity, neutrino masses, modified gravity can all be tested by looking at the shape of the clustering power spectrum), most of the meaningful astrophysical processes are constrained by the amplitude of the spectrum. Another method, recently discussed in Lidz & Taylor (2016); Cheng et al. (2016), exploits the fact that the interloper lines, being emitted at different redshifts respect to the targeted line, will introduce an anisotropic component in the power spectra, due to the incorrect redshift projection. A third method to mitigate contamination from different lines has been proposed by Visbal & Loeb (2010); Visbal et al. (2011), and involves the cross-correlation between maps measured at different frequencies, whose emission comes from atoms and molecules at the same redshift. Since all contaminant lines in each map will generally come from different redshifts, they will not contribute to the signal in the cross-correlation, but they will only add noise to the measurement. While cross-correlation measurements are generally more complicated to be carried out, most surveys proposed so far work in a broad frequency range, and multiple crosscorrelations produced at the same redshift among lines from different atoms and molecules might be attempted, at least in the non-linear regime. If the amplitudes of the lines to be cross-correlated is large enough, the information content from these measurements will be vast, and it will enable us to constrain various physical processes of the ISM. In this paper we propose the use of cross-correlation measurements among various emission lines from carbon, oxygen, and nitrogen to constrain the mean amplitude of each emission line at redshift z > 4. Using measurements of the Cosmic Infrared Background (CIB) angular power spectra from Herschel/SPIRE (Viero et al. 2013) and Planck (Planck Collaboration et al. 2014c), coupled to a compilation of star formation rate density (SFRD) measurements from Madau & Dickinson (2014), we constrain the galaxy FIR luminosity as a function of the halo mass at all relevant redshifts. By using scaling relations from Spinoglio et al. (2012) to link the intensity of emission lines to the constrained galaxy infrared luminosity, we compute 3D emission line power spectra for all relevant lines. Focusing on two experimental setups, corresponding to present and future ground-based surveys, we show that multiple cross-correlations with the [CII] line can constrain the mean amplitudes of all lines. This is important not only to constrain average properties of the ISM of galaxies at high redshift, but also because, as shown in Mashian et al. (2016);De Zotti et al. (2016); Carilli et al. (2016), especially the CO and [CII] line emission from galaxies across cosmic time distort the CMB spectrum at a level that must be taken into account by future space-based surveys aiming at measuring the tiny spectral distortions of the CMB, such as PIXIE. Intensity mapping, by constraining the mean amplitude of the signal, will healp dealing with this important foreground. In Sect. 2 we will derive the formalism used to compute emission line power spectra from the Halo model. We will then discuss in Sect. 3 the physics of the ISM in the context of emission lines from carbon, oxygen, and nitrogen, with particular focus on all possible cross-correlations to be performed using the experimental setups discussed in Sect. 5. Finally we will discuss our main results in Sect. 6. Throughout this paper, we adopt the standard flat ΛCDM model as our fiducial background cosmology, with parameter values derived from the best-fit model of the CMB power spectrum as measured by Planck Collaboration et al. (2014a).

A HALO MODEL FOR EMISSION LINE AMPLITUDES
The computation of 3D auto-and cross-power spectra of intensity line emission is performed in the context of a Halo model developed in Shang et al. (2012), where the galaxy luminosity is linked to the mass of the host dark matter halo with a simple parameteric form. It has been successfully applied to the interpretation of the latest measurements of angular CIB power spectra from Herschel/SPIRE (Viero et al. 2013) and Planck (Planck Collaboration et al. 2014c).
Using the latest measurements of CIB auto-and cross-power spectra at 250, 350 and 500 µm from Viero et al. (2013), together with a compilation of measurements of SFRD in the redshift range 0 < z < 6 (Madau & Dickinson 2014), we are able to constrain the galaxy infrared luminosity as a function of halo mass and redshift. We then use known scaling relations from Spinoglio et al. (2012) to compute the amplitudes of emission lines from carbon, oxygen, and nitrogen with respect to the constrained galaxy infrared luminosity. This allows us to compute the amplitudes of 3D power spectra for all relevant emission lines at all redshifts. This approach is very similar to that discussed in Cheng et al. (2016).
2.1. The Halo model for CIB anisotropies The halo model is a phenomenological description of the galaxy clustering at all angular scales (Cooray & Sheth 2002). Assuming that all galaxies live in virialized dark matter structures, called halos, and using a recipe to populate halos with galaxies, the clustering power spectrum results from the sum of two components: a 1-halo term, related to correlations between galaxies in the same halo, and responsible for the clustering at small angular scales, and a 2-halo term, which describes the power spectrum at large angular scales, and is due to correlations between galaxies belonging to separated dark matter halos. The angular power spectrum of CIB anisotropies, observed at frequencies ν and ν , is defined as: δ I lm,ν δ I l m ,ν = C l,νν δ νν δ mm where I ν is the specific intensity at that frequency, given by: here χ(z) denotes the comoving distance at redshift z, a(z) is the scale factor, and j(ν, z) is the comoving emission coefficient. In Limber approximation (Limber 1954), Eqs. 2 and 3 can be combined to give the clustering angular power spectrum as: where P νν (k, z) is the 3D power spectrum of the emission coefficient, expressed as: This term is composed by the mentioned 1-halo and 2-halo components. Thus, together with a scale independent shotnoise power spectrum, describing the contribution from random fluctuations due to the Poisson distribution of sources, the total CIB angular power spectrum is: This quantity will be computed and fit to Herschel/SPIRE measurements of CIB angular power spectra in order to constrain the galaxy infrared luminosity.
Below we show how to compute the two clustering terms. This formalism will be useful in Sect. 2.3, when computing 3D power spectra of emission lines. The mean emissivityj ν (z) from all galaxies is computed from the infrared galaxy luminosity function dn/dL as: where the galaxy luminosity L (1+z)ν is observed at the frequency ν with a flux given by: Neglecting any scatter between galaxy luminosity and dark matter halo mass, the luminosity of central and satellite galaxies can be expressed as L cen,(1+z)ν (M H , z) and L sat,(1+z)ν (m SH , z), where M H and m SH denote the halo and sub-halo masses, respectively. We can thus rewrite Eq. 7 as the sum of the contributions from central and satellite galaxies as:j here dN/dm (Tinker et al. 2008) and dn/dm (Tinker et al. 2010) denote the halo and sub-halo mass function respectively, while N cen is the number of central galaxies in a halo, which will be assumed equal to zero if the mass of the host halo is lower than M min = 10 11 M , and one otherwise. Introducing f cen ν and f sat ν as the number of central and satellite galaxies weighted by their luminosity, as: and the power spectrum coefficient of CIB anisotropies at the observed frequencies ν and ν can be written as the sum of a 1-halo term and 2-halo term as, respectively: where and u(k, M, z) is the Fourier transform of the Navarro-Frenk-White (NFW) density profile (Navarro et al. 1997), with concentration parameter from Duffy et al. (2010). The term b(M, z) denotes the halo bias (Tinker et al. 2010). The linear dark matter power spectrum P lin (k) is computed using CAMB (http://camb.info/). The final ingredient to be specified is the link between galaxy luminosity and host dark matter halo mass. Following Shang et al. (2012), we assume a parametric function, where the dependence of the galaxy luminosity on frequency, redshift, and halo mass is factorized in three terms as: The parameter L 0 is a free normalization parameter whose value is set by the amplitude of both the CIB power spectra and the SFRD. It has no physical meaning, and it will not be discussed further in the rest of the paper. A very simple functional form (see Blain et al. 2003, and reference therein) is assumed for the galaxy SED: where T d is the dust temperature averaged over the redshift range considered, and β is the emissivity of the Planck function B ν (T d ). We note that we discarded a redshift dependence of the dust temperature, because it is not very well constrained by the data. The power-law function at frequencies ν ≥ ν 0 has been found more in agreement with observations than the exponential Wien tail (see also Hall et As explained in Shang et al. (2012), a power law is motivated by the study of the star formation rate (SFR) per unit stellar mass, or specific star formation rate (sSFR). Assuming that the stellar mass to halo mass ratio does not evolve substantially with redshift, the ratio of galaxy infrared luminosity L IR to halo mass has an evolution similar to the sSFR, thanks to the correlation between SFR and infrared luminosity (Kennicutt 1998).
where M eff describes the most efficient halo mass at hosting star formation, while σ L/m accounts for the range of halo masses mostly contributing to the infrared luminosity. Such a functional form captures the fact that, for halo masses much lower and much higher than M eff , various mechanisms prevent an efficient star formation (Benson et al. 2003;Silk 2003;Bertone et al. 2005;Croton et al. 2006;Dekel & Birnboim 2006;Béthermin et al. 2012a;Behroozi et al. 2013).

Analysis
We perform a Monte Carlo Markov Chain (MCMC) analysis of the parameter space, using a modification of the publicly available code CosmoMC (Lewis & Bridle 2002), and fitting to six CIB auto-and cross-power spectra from Viero et al. (2013) in the multipole range 200 < l < 23000. We also add a dataset for the SFRD as a function of redshift by averaging multiple measurements, discussed in Madau & Dickinson (2014), in eleven redshift bins in the range 0 < z < 6. We vary the following set of parameters: and we add six free parameters A i=1,...6 to model the amplitudes of the CIB shot-noise power spectra. All parameters have a uniform prior, and we fix the emissivity index to β = 1.5 (Planck Collaboration 2014), and σ 2 L/M = 0.5 (Shang et al. 2012;Planck Collaboration et al. 2014c). With a total χ 2 value of 104.9 for 97 degrees of freedom, we obtain a very good fit to the data. In Table 3, we quote mean values and marginalized limits for all free parameters used in the fit, while in Fig. 1 we plot the Herschel/SPIRE measurements of the CIB power spectra, together with our best estimates of the 1-halo, 2-halo, shot-noise, and total power spectrum. It is important to note that there is a relevant uncertainty associated to measurements of the SFRD, especially at the high redshifts considered in this work. The compilation of measurements extrapolated from Madau & Dickinson (2014) (plotted in Fig. 2), is based on galaxy counts, and there is a number of uncertain steps in the conversion from galaxy counts and luminosities to star formation rates, mainly related to assumptions on conversion factors and dust attenuation. When considering clustering meausurements, the Planck Collaboration, using a Halo model similar to the one presented in this paper, and fitting to CIB power spectra between 217 GHz (1381 µm) and 857 GHz (350 µm) in the multipole range 50 < l < 2000, infer a much higher SFRD at high redshifts (Planck Collaboration et al. 2014c), respect to the values found here by fitting Herschel-SPIRE data and star formation rate density data from Madau & Dickinson (2014) 2014c)). The reason for this discrepancy is mainly due to the different values inferred for the parameter δ in Eq. 17. Planck Collaboration et al. (2014c) found δ = 3.6 ± 0.2 (see Table 9 of Planck Collaboration et al. (2014c)), while we find δ = 2.6 ± 0.2, compatible with (Viero et al. 2013). We checked that the fitting to star formation rate density data from Madau & Dickinson (2014) is not responsible for such a divergence, by performing an MCMC run with only one measurement of the local SFRD at z = 0.07 from Madau & Dickinson (2014) (thus being compatible with Planck's analysis, since they use a prior on the local SFRD from Vaccari et al. (2010)). As it is clear from Fig. 2 (2014c) can be explained by a combination of multiple factors involving our ignorance of the exact values of some key parameters, such as the amplitudes of the shot noise power spectra and the redshift evolution of the galaxy luminosity, coupled to differences in the datasets considered. CIB anisotropies are mostly sourced by galaxies at redshift 1 < z < 4 and, in this range, a simple power law might not be a good description of the redshift evolution of the galaxy luminosity/halo mass relation. Some semianalytic models of galaxy formation and evolution find a power law slope of ∼ 2. Shot noise for 350x500 µm 2659 ± 80 S 500x500 Shot noise for 500x500 µm 1600 ± 61  2008), but also a more gradual evolution, with different slopes for low redshift and high redshift sources (Wu et al. 2016). On the other end, observations are more in agreement with a steep evolution with redshift (Oliver et al. 2010), or with a steep evolution followed by a plateau for z ∼ 2 (Bouché et al. 2010;Weinmann et al. 2011), which is also not easily explained by theoretical arguments (Bouché et al. 2010;Weinmann et al. 2011). Planck Collaboration et al. (2014c is indeed able to find lower values for the star formation rate density at early times, more in agreement with this work, but only when they  2014c)). The differences between the two datasets in terms of angular scales and related uncertainties can also be responsible for the difference values inferred for δ. Planck data probe CIB anisotropies at large scales with very high precision. However, because of its angular resolution, Planck is not able to access multipoles higher than l ∼ 2000, where the 1-halo term and the shot-noise dominate the clustering, and are degenerate. Uncertainties in the contribution of these two terms to the small-scale clustering (Planck Collaboration et al. (2014c) used free amplitudes for the shot-noise power spectra, with flat priors based on current measurements, such as, e.g., Béthermin et al. (2012b)) translates in an uncertainty in the inferred constraints on the halo model parameters. On the other end, Herschel/SPIRE data probe both large and small scales, but while adding information at small scales helps disentangling the relative contributions to the total power from the 1halo term and the shot-noise, the largest scales are measured with much larger uncertainty than Planck. Finally, Planck and Herschel probe a different frequency range, which might affects results. Thus, it is possible that the differences in the datasets used, coupled with uncertainties regarding the level of the shot-noises, and a poor description of the redshift evolution of the sources, determine different values for the parameter δ. It is clear that the higher the value of the star formation rate density, the greater the value for the mean emission from all atoms and molecules. This would translate in large amplitude for the emission line power spectra. In order to be as independent as possible on the particular values of the Halo model parameters used to constrain the galaxy infrared luminosity, we compute predictions for the 3D power spectra of emission lines using both the mean values found by fitting Herschel/SPIRE data (quoted in Table 3) and the mean values quoted in Table 9 of Planck Collaboration et al. (2014c). The geometric average of these two estimates will be our best estimate of the power spectrum of the emission lines. In the rest of the paper we will focus on predictions based on these average estimates of the power spectra. In Fig. 3

Intensity mapping power spectrum from the Halo model
The analysis presented in the previous section has been necessary to constrain the main parameters describing the galaxy SED and its dependence on halo mass and redshift. The galaxy infrared luminosity is: where the extremes of integration correspond to the wavelength range 8 < λ < 1000 µm. We can use scaling relations provided in Spinoglio et al. (2012), to express the emission line luminosity I α (where α denotes emission lines from the atoms and molecules considered, such as carbon, oxygen, nitrogen) as a function of the constrained infrared luminosity, as: where all luminosities are in units of 10 41 erg s −1 . These scaling relations are obtained from a sample of local galaxies compiled by Brauher et al. (2008) using all observations collected by the LWS spectrometer (Clegg et al. 1996) onboard ISO (Kessler et al. 1996). Regarding the [NII] 205 µm emission line, whose luminosity is not found in Spinoglio et al.  .
In Table 1 we summarize the values used for slopes, intercepts and their uncertainties, together with their associated transition, and transition temperatures from Kaufman et al. (1999); Cormier et al. (2015). The emission line luminosity at each redshift for each halo mass can now be expressed as previosly done for the galaxy luminosity (see Eq. 15) as: where the term F(M, z) contains the global dependence on redshift and halo mass as and we use the parameter values from Table 3 to compute the term F(M, z). This functional form allows us to link the emission line luminosity of a galaxy to its host halo mass, and to evolve the amplitude of all emission lines with redshifts.
We note that this model assumes that the redshift evolution of all emission lines is the same, since it follows the evolution of the galaxy infrared luminosity (through the parameter δ). Different emission lines might have different a evolution with redshift, and more sophisticated models could incorporate redshift-dependent scaling relations for each line. However, current data do not allow us to constrain the exact dependence on redshift of each emission lines. Thus, to keep the analysis as simple as possible, we do not consider such a scenario.
It is easy to see that, assuming that each halo hosts only one galaxy (a good approximation because, at high redshift, halos are not very massive, see also Lidz et al. (2011)), and in the limit of sufficiently large scales (so that the NFW profile approaches unity), the clustering auto-power spectrum of emission line α can be written as: where Introducing an effective, scale independent, bias term as: the clustering power spectrum of emission line α can be expressed as: where the average specific intensityĪ α (z) is: and z α denotes the redshift of emission of the atom or molecule α. Analogously, the shot-noise power spectrum can be expressed as:

THE PHYSICS OF THE ISM WITH EMISSION LINES AND EMISSION LINE RATIOS
Understanding the main heating and cooling processes of the ISM is a key goal of astronomy, because they play a fundamental role in the formation of stars, and thus in the galaxy evolution. Space missions such as Planck and Herschel, together with the Stratospheric Observatory for Infrared Astronomy (SOFIA) and the Atacama Large Millimeter Array (ALMA), are now giving new insights on these physical processes, providing spatially resolved maps of the interstellar dust in our Galaxy, and measuring atomic and molecular emission lines from the main phases of the ISM both in the Milky Way (Pineda et al. 2013(Pineda et al. , 2014Goicoechea et al. 2015), and in external galaxies (see e.g.  2016)). The gas in the ISM of galaxies is observed in three main phases; a cold and dense neutral medium (T≥ 50K) is in rough pressure equilibrium (with P/k ∼ 10 3 − 10 4 Kcm −3 ) with a hot (T≥ 10 6 K), ionized phase, and an intermediate, warm (T≥ 8000K) phase, which can be either neutral or ionized, depending on the gas density (Wolfire et al. 1995).
Various mechanisms contribute to the heating and cooling of the ISM. For a gas with hydrogen density n, temperature T, cooling rate per unit volume of Λ(T), and heating rate per unit volume of Γ(T), the thermal balance between heating and cooling is expressed in terms of a Generalized Loss Function L: For a gas at constant thermal pressure nT, equilibrium occurs when L = 0 and the explicit form for Λ and Γ depends on the heating and cooling process considered, as explained below. The investigation of the thermal balance and stability conditions of the neutral ISM started with Field et al. (1969), who first presented a model of the ISM based on two thermally stable neutral phases, cold and warm, heated by cosmic-rays. Subsequent analyses by many authors focused on the heating provided by the photoelectric ejection of electrons from dust grains by the interstellar radiation field (Draine 1978;Wolfire et al. 1995;Kaufman et al. 1999;Wolfire et al. 2003).
Most of the Far-Ultraviolet (FUV) starlight impinging on the cold neutral medium is absorbed by dust and large molecules of polycyclic aromatic hydrocarbons (PAH), and then reradiated as PAH infrared lines and infrared continuum radiation. However, as pointed out by Tielens & Hollenbach (1985b), in photodissociation regions (PDRs), the photoelectric heating of dust grains provides an efficient mechanism (0.1% − 1%) at converting the FUV heating into atomic and molecular gaseous line emission. The physics of heating processes in PDRs can be understood in terms of a limited set of parameters, namely the density of Hydrogen nuclei density n and the incident FUV (6eV < hν < 13.6 eV) parametrized in units of the local interstellar field, G 0 (Tielens & Hollenbach 1985a,b;Kaufman et al. 1999), in units of the Habing field (1.6 · 10 −3 ergs cm −2 s −1 ). The basic mechanism for gas heating and cooling is the following: about 10% of incident FUV photons eject photoelectrons from dust grains and PAH molecules, which cool by continuum infrared emission. The photoelectrons (with energy of about 1 eV) heat the gas by collisions, and the gas subsequently cools via FIR fine-structure line emission. The entire process thus results in the conversion of FUV photons to FIR continuum emission plus spectral line emis-sion from various atoms and molecules. As an example, the computation of the heating due to small grains is given by (Bakes & Tielens 1994): the radiation field G 0 quantifies the starlight intensity, and is the fraction of FUV photons absorbed by grains which is converted to gas heating (i.e. heating efficiency), and it depends on G 0 T 1/2 n e , where n e denotes the electron density Wolfire et al. (1995). A detailed calculation of the main heating processes in the ISM, including the effect from photoelectric heating, cosmic rays, soft X-rays, and photoionization of CI is presented in Wolfire et al. (1995); Meijerink & Spaans (2005). The cooling rate Λ of each atom/molecule depends on both the number density and the equivalent temperature of each species. A recent estimate of the cooling rate of the [CII] line for temperatures between 20 K and 400 K is (Wiesenfeld & Goldsmith 2014): where n(C + ) denotes the carbon number density, and T kin the kinetic temperature of the gas. Numerical codes compute a simultaneous solution for the chemistry, radiative transfer, and thermal balance of PDRs, providing a phenomenological description of the interplay among three main parameters n, G 0 and T (see e.g. Kaufman et al. (1999)) for all emission lines. The observed intensity of line emissions can thus be compared with models to constrain these parameters. with characteristics similar to PIXIE will be able to detect most of the main cooling lines from both PDRs and from the ionized medium of high redshift galaxies. Below we summarize some useful diagnostics of the ISM provided by these important lines (see also Cormier et al. (2015)).
• [CII] emission line: It is hard to overestimate the importance of the [CII] emission line in constraining physical properties of the interstellar medium. Because of its low ionization potential, the [CII] line arises both from ionized and neutral gas. In PDRs, the low gas critical density for collisions with Hydrogen and the low excitation temperature for the [CII] 2 P 3/2 − 2 P 1/2 transition (only 92 K, see Table 1), make C + one of the major coolant of the neutral ISM. Moreover, since the [CII] line is generally one of the brightest lines in starforming galaxies, it is potentially a very strong indicator of star formation rate (SFR) (Boselli et al. 2002 Hudson & Bell (2004). Being in the same ionization stage, their ratio directly determines the electron density of the ionized gas in HII regions. For electron densities n e larger than 10 cm −3 , the 122/205 µm line ratio R 122/205 increases as a function of n e , starting from R 122/205 ∼ 0.6 for n e ∼ 10 cm −3 , and reaching the value R 122/205 ∼ 3 (the value used in this paper) for n e ∼ 100 cm −3 (Tayal 2011;Goldsmith et al. 2015).

Moreover, combined measurements of line emission from [NII] and [CII]
can be used to estimate the amount of [CII] emission coming from the ionized medium (Malhotra et al. 2001;Oberst et al. 2006;Decarli et al. 2014;Hughes et al. 2016). Recently Goldsmith et al. (2015), using data from the PACS and HIFI instruments onboard Herschel, estimated that between 1/3 and 1/2 of the [CII] emission from sources in the Galactic plane arise from the ionized gas. The [NII]/[CII] ratio is also useful to estimate the metallicity of a galaxy (Nagao et al. 2012). Finally, the [NII] emission lines, arising from gas ionized by O and B type stars, directly constrains the ionizing photon rate, and thus the star formation rate (Bennett et al. 1994;McKee & Williams 1997).
• Oxygen 63 µm and 145 µm lines: Oxygen has a ionization potential of 13.62 eV, just above that of hydrogen. 145 µm respectively) they contribute significantly to the cooling of the ISM only for high FUV fields and/or high densities. The measurement of the mean amplitude of the [OI] lines with intensity mapping would give us clues regarding the mean value of the G 0 field and the mean density of PDRs at high redshifts (Meijerink et al. 2007).
The intensity mapping technique would constrain the mean amplitude of multiple emission lines, together with their ratio, thus probing mean properties (such as mean radiation field, mean electron density in HII regions, mean density of various atoms, molecules) at high redshifts.

MULTIPLE CROSS-CORRELATIONS CONSTRAIN THE PHYSICS OF THE ISM
As previously stated, the cross-correlation signal between different emission lines coming from the same redshift is important not only to avoid contamination from foreground lines (assuming that, at the frequencies considered in the crosscorrelation measurements, foregrounds are not correlated), but also to help constraining the mean amplitude of each signal. This is particularly true at sufficiently small scales, where the SNR is larger. If we assume that all lines are emitted by the same objects (a reasonable assumption, especially if the emission lines are not distant from each other, as in the case of the FIR lines such as [CII], [NII] and [OI]), it will be possible to constrain the mean amplitudes of emission lines I i=1,...N , just by looking at all cross-correlation power spectra 5 . For a survey working in a given frequency range where N lines are detected, there are N(N − 1)/2 cross-correlation measurements to be performed and, assuming there is perfect correlation among lines, it is sufficient that N ≥ 3 to be able to constrain the mean emission from all lines. The chances of detecting auto-and cross-power spectra strongly depend on the amplitude of the spectra, which, as already seen, is very uncertain. In the following we will consider predicted measurements of multiple combinations of emission line power spectra for two different surveys. The first one corresponds to a survey of the [CII] emission line similar to the proposed CONCERTO. The second one, referred in literature as CII-Stage II, and described in Silva et al. (2015); Lidz & Taylor (2016), is more sensitive, and corresponds to an evolution of currently planned [CII] surveys. As already emphasized, emission line power spectra are strongly contaminated by interlopers lines emitted by molecules at different redshifts. In case of [CII], the main confusion results from foreground emission of CO molecules undergoing rotational transitions between states J and J-1.
As an example, [CII] emission from z = 6 is observed at frequency ν obs = 271.6 GHz, and it is mainly contaminated by CO rotational transitions J = 3 → 2 (z = 0.27), J = 4 → 3 (z = 0.70), J = 5 → 4 (z = 1.12), J = 6 → 5 (z = 1.54), and J = 7 → 6 (z = 1.97). Emission lines beyond this transition have a negligible contribution to the total foreground due to CO molecules, and we will not consider them in the rest of the paper. Using linear scaling relations from Visbal & Loeb (2010) to express the amplitude of the various CO emission lines as a function of the infrared luminosity, it is possible to estimate the contamination due to the main CO rotational lines. In the following, when plotting the [CII] auto-power spectra at various redshits, we will also plot the CO auto-power spectrum computed as the sum of the the main CO rotational transitions involved (from 3 → 2 to 7 → 6), in order to highlight the amplitude of this foreground.

EXPERIMENTAL SETUPS AND PREDICTIONS
In order to measure high-redshift fluctuations with sufficient SNR at the scales of interest, it is important to optimize the survey area. All predictions considered in this section are based on measurements spanning a redshift range ∆ z ∼ 0.6 which corresponds to a frequency range of B ν ∼ 20 GHz at z = 7 for the [CII] line. We follow Gong et al. (2012) to compute uncertainties on the power spectra. The primary goal of the first survey considered, called CON-CERTO, is to detect [CII] fluctuations in the redshift range 4.5 < z < 8.5. It is based on a spectrometer working in the frequency range 200 < ν < 360 GHz, with spectral resolution δ ν ∼ 1.5 GHz. Such a frequency window imposes the use of a so-called "sub-millimetre" telescope, with primary aperture size D = 12 m, and moderate angular resolution. The instrumental noise is thus computed for a total observing time of t survey = 1500 hours, and a number of spectrometers N sp = 1500. The survey area considered here is two square degrees, and is optimized to ensure high SNR in the wavenumber range of 0.1 < k < 1 h / Mpc. Accounting for realistic observational conditions and the total atmospheric transmission, the Noise Equivalent Flux density (NEFD), computed as the sensitivity per single pixel divided by the square root of the number of spectrometers, is equal to NEFD = 155 mJy sec 1/2 , for a spectral resolution of δ ν = 1.5 GHz. The on-sky sensitivity σ N can be expressed as: where is the beam area (in steradians), and the beam FWHM is given by: where λ obs is the observed wavelength. Values for σ N at z = 5, z = 6, and z = 7 are 15, 11, and 8.3 MJy/sr √ (sec) respectively. The observing time per pixel is given by: where ∆Ω survey is the total survey area covered. Assuming a spherically averaged power spectrum measurement, and a directionally independent on sky sensitivity σ N , the variance of the power spectrum is: where N m (k, z) denotes the number of modes at each wavenumber: the term ∆ k is the Fourier bin size, and V s (z) is the survey volume, expressed as: The averaged noise power spectrum in Eq. 37 is: where the volume surveyed by each pixel is:  withȳ and λ α is the wavelength of the line α is the rest frame. In Fig. 4  [CII] auto-power spectrum will be detected with high significance (SNR > 50), while the [CII] cross-correlations with oxygen and nitrogen at these scales will not be very significant (SNR∼ 3 and SNR∼ 0.5 respectively). However, considering smaller scales (larger wavenumbers) the SNR increases significantly, and it will enable us to constrain the mean quantities I [CII] , I [OI] , and I [NII] . Given the CONCERTO frequency coverage, at z = 6.0 it is possible to add the cross-correlation   with [NII] (122 µm). As shown in Fig. 5, the cross-correlation of carbon with oxygen and nitrogen seems to be barely detectable at linear scales. However, in the non-linear regime, it might still be possible to measure these cross-correlations, and thus constrain the mean amplitude of these emission lines. As already described in Sect. 3, by looking at the cross-power spectra [CII]x[NII] (121.9 µ)m, and [CII]x[NII] (205.2 µm), we would be able to measure the mean ratio [NII] (205.2 µm) / [NII] (121.9 µm), which is useful not only to constrain the electron density of the low-ionized gas in HII regions, but also to infer the mean emission of [CII] from PDRs, and to constrain the global star formation rate. The mean ratio between [OI] (145.5 µm) and CII is also a useful diagnostic of mean properties of properties of PDRs, such as the hydrogen density and the strength of the radiation field. The second experimental setup, called CII-Stage II, has been introduced in Silva et al. (2015) as an appropriate baseline to ensure detection of [CII] spectra in case of a pessimistic [CII] amplitude (see also Lidz & Taylor (2016)). It consists of a dish with diameter D = 10 m, with 16000 bolometers and N sp = 64 beam spectrometers, observing in the frequency range 200 < ν < 300 GHz, with a frequency resolution of 0.4 GHz. The total survey area is 100 deg 2 for a total observing time of t survey = 2000 hours, and a NEFD of 5 mJy sec 1/2 . As it appears from Fig. 6, the cross-correlation of carbon with oxygen and nitrogen is now detectable with high SNR at z = 6. A space-based survey, being not limited by the atmosphere, would be able to perform measurements on a still wider frequency range, and thus perform measurements of high-redshift correlations with other interesting lines such as [

DISCUSSION
We have developed a consistent framework to compute predictions of 3D power spectra of multiple FIR cooling lines of the ISM. Using measurements of CIB power spectra, together with measurements of star formation rate density from Madau & Dickinson (2014), it is possible to constrain the galaxy FIR luminosity at all redshift, which can be directly linked to emission line amplitudes through scaling relation from Spinoglio et al. (2012). Present and upcoming groundbased surveys aiming at measuring the power spectrum of the bright [CII] line, should be able to detect also the crosscorrelation between the [CII] line and other lines produced in all phases of the ISM, such as [NII] (122 µm and 205 µm), and [OI] (145.5 µm). Multiple measurements of cross-power spectra between [CII] and other emission lines will allow us to constrain the mean amplitude of each signal, and they will be key to gain insight into the mean properties of the ISM. Future surveys, such as PIXIE (Kogut et al. 2011(Kogut et al. , 2014, working in a broad frequency range, will detect many more atomic and molecular lines emitted from moderate to high redshift with high SNR, allowing us to obtain multiple probes of all phases of the ISM. Moreover, the cross-correlation of the target line with galaxy number densities from future surveys such as, e.g., LSST (LSST Science Collaboration et al. 2009), will be a powerful method to eliminate line foregrounds. Line emissions from multiple atoms/molecules at multiple redshifts are also an important foreground for future surveys aiming at constraining CMB spectral distortions. In Fig. 7 we plot µ-type and y-type spectral distortions with µ = 5 · 10 −8 and y = 1 · 10 −8 , corresponding to the current PIXIE 5σ sensitivity limits, together with the sum of the spectra from carbon monoxide emission lines (from J = 1 → 0 to J = 7 → 6), and the spectra from all emission lines considered in this work. The CO spectra have been computed using scaling relations from Visbal & Loeb (2010) to link the CO line emission to the star formation rate, and the Kennicutt relation to express the star formation rate in terms of the galaxy infrared luminosity (Kennicutt 1998). The amplitude of the global signal from CO lines is similar to what found by Mashian et al. (2016) using a radiative transfer modeling technique, even if the shape is slightly different. We note that, even if foreground lines do not have a simple spectral dependence, unlike other foregrounds that can be modeled with power law such as synchrotron or thermal dust, their shape is still monotonic in frequency, and thus very different with respect to the CMB spectral distortions. However, foreground subtraction will require a very good knowledge of the amplitude and shape of the total signal provided by the sum of these lines. The intensity mapping technique, by constraining the mean amplitude of the signal in multiple redshift bins, will help constraining the global contamination signal. Finally, it is clear that an aggressive program to model the amplitude of all emission lines at all redshifts is necessary to have a detailed interpretation of upcoming measurements. Scaling relations are useful to work with, but they provide little information on the main physical mechanisms governing the line emission. Moreover, they are based on few observations performed at some given redshift, and their redshift evolution is not very well known. Different physical conditions can dominate the line emission at different epochs, strongly affecting the amplitude of the signal. As an example, at high redshift, the CMB strongly suppresses the [CII] emission from the cold neutral medium, leaving only the emission from PDRs (Vallini et al. 2015). The redshift evolution of the galaxy infrared luminosity (which governs the evolution of the line emission in our model) is determined by the power law parameter δ (see Eq. 17), which, as stated earlier, is quite uncertain, especially at high redshifts. On the other hand, semi-analytic models of galaxy formation and evolution often involve a large number of assumptions and free parameters, and such a complexity makes them difficult to use. A third approach, intermediate between the two, and based on present and upcoming measurements from, e.g. ALMA and SOFIA, should be developed to model the line intensity of all relevant emission lines, together with their redshift evolution. Such a model, possibly based on the physics of photodissociation regions, ionized medium, and molecular clouds, will offer an important guidance in interpreting upcoming and future intensity mapping observations, and thus constrain the mean properties of high-redshift galaxies.