EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 609 (January 2018) A&A, 609 (2018) A130 Full HTML
Free Access
Issue
A&A
Volume 609, January 2018
Article Number A130
Number of page(s) 19
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/201732019
Published online 05 February 2018

© ESO, 2018

1. Introduction

One of the final frontiers in piecing together a coherent picture of cosmic history relates to the period 300–900 million years after the Big Bang (redshifts 6 <z< 15). During this time, the Universe underwent two major changes. Firstly, the earliest stars and galaxies began to shine, bathing the Universe in starlight. Secondly, the intergalactic medium transitioned from a neutral to a fully ionized gas, a timespan known as the epoch of reionization (EoR). Connecting these two changes is highly desirable and after years of effort, recent breakthroughs showed that reionization occured at 6 <z< 10 (Planck Collaboration XLVII 2016) and that UV-selected star-forming galaxies likely dominated the reionization process (e.g. Robertson et al. 2015). Active galactic nuclei can also potentially contribute to reionization (Giallongo et al. 2015); the exact role of the two populations is still unclear.

Another remarkable result of cosmology in the last decade is the realization that the star formation rate (SFR) density at redshifts z> 1 is higher than at present by about an order of magnitude and that half of the energy produced since the surface of last scattering has been absorbed and reemitted by dust (Dole et al. 2006), in dusty star-forming galaxies (DSFG). Most of the light produced at high redshift thus reaches us in the wavelength range 100 μm–1 mm (Lagache et al. 2005). Contribution of DSFG to the global star formation history is roughly known up to z = 3 (Madau & Dickinson 2014). But at higher redshifts and in the EoR, it is an uncharted territory. At such early epochs (z> 5) dust is surely present even if in small amounts (Riechers et al. 2013; Watson et al. 2015) and can strongly affect SFR measurements based on UV-luminosity.

With the advent of the Atacama Large Millimetre Array (ALMA) and NOEMA, it is now possible to measure the dust content of very high redshift galaxies, but also to use far-infrared fine-structure lines (as [OIII] or [CII]) to study the physical conditions of their interstellar medium (ISM). The [OIII] line, originating from diffuse and highly ionized regions near young O stars, is a promising line (Inoue et al. 2016) that might gain in importance in low-metallicity environments where photo-dominated regions (PDRs) may occupy only a limited volume of the ISM. The [CII] line, predominantly originating from PDRs at high redshift (Stacey et al. 2010; Gullberg et al. 2015), can provide SFR estimates that are not biased by dust extinction, although it has been found to depend strongly on the metallicity (Vallini et al. 2015; Olsen et al. 2017). This line can also be used to measure the systemic redshift of the galaxies (e.g. Pentericci et al. 2016). In addition, the [CII]-line ALMA surveys will derive the line luminosity functions, thus measuring the abundance and intensity distributions of [CII] emitters (Aravena et al. 2016).

Due to its relatively low ionization potential, [CII] is the dominant form of the element under a large variety of conditions. The C+ ion has only two fine structure levels in the ground electronic state. The lower J = 1 /2 level has statistical weight gl = 2. The upper J = 3/2 level has statistical weight gu = 4, and lies at equivalent temperature T∗ = ΔE/k = 91.25 K above the ground state. The measured transition frequency is 1900.537 GHz (Cooksy et al. 1986) corresponding to a transition wavelength of 157.74 μm, making the [CII] line easily accessible from the ground for 4.5 ≲ z ≲ 8.5. These redshifts marks an important epoch when the ISM in typical galaxies matures from a nearly primordial, dust-free state at z ~ 8, during the EoR, to the dust- and metallicity-enriched state observed at z ~ 4.

Consequently, we investigate in this paper the correlation between SFR, [CII] luminosity and metallicity, and predict the luminosity function of [CII] line emitters at z ≥ 4. We use the semi-analytical model (SAM) described in Cousin et al. (2015b), that we combine with the CLOUDY photoionisation code (Ferland et al. 2013, 2017). For each galaxy in the SAM (that has its own mass, SFR, metallicity, size, etc) we define an equivalent PDR characterised by its own properties (i.e. interstellar radiation field, gas metallicity, mean hydrogen density) and run CLOUDY to derive its [CII] emission, taking into account the CMB (heating and attenuation). We are well aware that using global galaxy characteristics to predict the [CII] line emission ignores the complex properties of galaxies at very high redshift in which differential dust extinction, excitation and metal enrichment levels may be associated with different subsystems assembling the galaxies (e.g., Carniani et al. 2017; Katz et al. 2016; Pallottini et al. 2017b). Complex hydrodynamical simulations are being undertaken (e.g. Olsen et al. 2017) but future developments and more statistics are needed to make detailed comparisons with observations (see the discussion in Katz et al. 2016). In the meantime, the low computational cost of SAMs makes them a powerful tool to model large volumes of the sky and to sample a large diversity of galaxy properties.

The paper is organized as followed: in Sect 2, we present briefly our SAM and validate its use for predicting the [CII] emission at very high redshift. Then, we describe our model for [CII]-line emission, and we quantify the effects of cosmic microwave background (CMB) and galaxy properties (as gas metallicity) on [CII] luminosity (Sect. 3). We discuss in Sect. 4 the L[CII]–SFR relation, and compare it with recent observations. Section 5 is dedicated to the [CII] deficit. In Sect. 6 we present the [CII] luminosity function from z = 4 to 8 and discuss its evolution. Finally, we conclude in Sect. 7. Throughout the paper, we use Chabrier (2003) initial mass function.

2. Galaxy formation in the early Universe

2.1. Brief description of the semi-analytical model

We used the SAM presented in Cousin et al. (2015b,a, 2016). In addition to the original prescriptions detailed in Cousin et al. (2015b), the extension of the model described in Cousin et al. (2016) tracks the metal enrichment in both the gas phase and stellar populations, which is essential to predict the [CII] emission. The chemodynamical model is applicable from metal-free primordial accretion to very enriched interstellar gas contents.

The SAM is combined with dark-matter merger trees extracted from a pure N-body simulation. The simulation is based on a WMAP-5yr cosmology (Ωm = 0.28, ΩΛ = 0.72, fb = 0.16, and h = 0.70) and covers a volume of [100 /h] 3 Mpc with 10243 particles. Each particle has a mass mp = 1.025 × 108M⊙. Haloes and sub-structures (satellites) are identified using HaloMaker (Tweed et al. 2009).

Dark matter haloes grow following a smooth accretion, with a dark-matter accretion rate Ṁdm derived from particles that are newly detected in the halo and that have never been identified in an other halo. Baryons are then progressively accreted following (1)where is the effective baryonic fraction depending on the virial halo mass and redshift. This fraction is computed following Gnedin (2000) and Kravtsov et al. (2004) photoionization models but with an effective filtering mass as defined in Okamoto et al. (2008).

Our SAM assumes a bimodal accretion (Khochfar & Silk 2009; Benson & Bower 2011), based on a cold and a hot reservoirs that are both fed by the metal-free cosmological accretion. In addition, the hot reservoir receives the galactic metal-rich ejecta. As the metallicity of the wind phase depends on the galaxy metal enrichment process, the metallicity of the hot reservoir evolves with time. Metals are initially formed by stars in the galaxies. The enriched gas is then ejected by supernova and active galactic nuclei (AGN) feedback (see Cousin et al. 2015b, for the detailed implementation of the supernovae and AGN feedback).

The chemodynamical model (Cousin et al. 2016) can follow the 1H, 4He, 12C, 14N, 16O, and 56Fe elements in the gas phase. Their production in stars and re-injection in the ISM are taken into account for stars with initial mass between 0.1 M⊙ and 100 M⊙ and for metal-free to super-solar metal fraction. It is assumed that stars are formed following a Chabrier (2003) initial mass function.

One of the particularity of the Cousin et al. (2015b,a) model is that the freshly accreted gas is assumed to be no-star-forming. It is progressively converted into star-forming gas and then into stars. The main idea behind the existence of the no-star-forming gas reservoir is that only a fraction of the total gas mass in a galaxy is available to form stars. The reservoir generates a delay between the accretion of the gas and the star formation. In the present paper, we use the conversion between the no-star-forming and star-forming gas as described in Cousin et al. (in prep.), thus assuming an inertial turbulent cascade in the gas. This updated version of Cousin et al. (2015b) SAM is called G.A.S., for Galaxy Assembly from dark-matter Simulations.

Cousin et al. (in prep.) present the model for dust extinction and emission. This modelling is not used for the [CII]-line emission prediction but is mandatory to compute the UV and IR luminosities of G.A.S. galaxies. Stellar spectra are based on Bruzual & Charlot (2003) library. Extinction curves and dust spectral energy distribution are computed using DustEM (Compiègne et al. 2011) and are self consistently applied to the disc and the bulge of the galaxy. A standard slab geometry for old stars in the disc is used (Guiderdoni & Rocca-Volmerange 1987). Additional extinction from burst clouds is applied for young stars in disc using a screen geometry (Charlot & Fall 2000). For the bulge, a standard Dwek geometry is used (Devriendt et al. 1999, and references therein). Effective extinctions predicted by this model are in excellent agreement with Calzetti et al. (2000) extinction law.

2.2. High-redshift stellar-mass and UV luminosity functions

The G.A.S. model is quite successful in predicting a vast number of observations from z = 4 to z = 0 (see Cousin et al. 2015a), including the stellar mass function, stellar-to-dark-matter halo mass relation, SFR density, stellar mass density, and specific SFR. Also, it reproduces well the stellar mass to gas-phase metallicity relation observed in the local universe and the shape of the average stellar mass to stellar metallicity relations (Cousin et al. 2016).

thumbnail Fig. 1

Observed and predicted stellar mass functions from z = 4 to z = 7. Observations are from Song et al. (2016), Duncan et al. (2014), Grazian et al. (2015), and Caputi et al. (2015). Our SAM predictions are shown in grey.

In this section, we extend the comparison between the model and observations to z> 4 to check the model validity at very high redshift z ~ 4–9. At such high z, stellar mass functions and UV-luminosity functions are the only observables that can be used.

We first compare the model prediction with the stellar mass functions (SMF). Stellar mass assembly is one of the most fundamental property of galaxy evolution, that does not depend in a SAM on for example complex metal-dependent extinction curve. SMF has been measured up to z ≃ 7, although with a quite large dispersion in the data points at z ≃ 7. The comparison between model predictions and observations is shown in Fig. 1. We have an overall excellent agreement between the two.

We show on Fig. 2 the model prediction for the UV luminosity function together with the most recent observations at z = 4 (the comparison at higher z is similar). The model has no constraint at the faint end due to our mass resolution (our model contains the contribution of all galaxies only for M⋆ ≥ 107M⊙). At the bright end, we limit the comparison when the number of galaxies in the simulation is >5 in the given luminosity bin. Observations are not corrected for extinction so we show the model prediction with and without extinction corrections. We can see that such corrections are important only for MUV ≤ −19. We have a very good agreement between the model and observations up to z ≃ 6. At z ≃ 7–8, our model slightly overestimates the number of MUV ≤ −21 objects. This may be caused by an underestimate of extinction corrections, which are very large for bright-UV galaxies.

These comparisons between model predictions and observations at z ≳ 4 give us confidence in using our SAM as a reference model to predict the [CII]-line emission.

thumbnail Fig. 2

Observed and predicted UV luminosity functions at z ~ 4. Observations are from Bouwens et al. (2015), Finkelstein et al. (2015). Our SAM predictions are shown in grey (with and without extinction correction, solid and dotted lines, respectively).

3. Modelling [CII] emission

[CII] line from high-z galaxies has been computed both through numerical simulations (e.g. Nagamine et al. 2006; Vallini et al. 2013, 2015; Olsen et al. 2017) and semi-analytical models (e.g. Gong et al. 2012; Muñoz & Furlanetto 2014; Popping et al. 2016). Here we take advantage of the excellent agreement of our SAM predictions at z> 4 with current constraints to revisit the expected [CII] signal from high-z galaxies.

3.1. Origin of [CII] emission in distant galaxies

The single [CII] fine structure transition is a very important coolant of the atomic ISM and of PDRs in which carbon is partially or completely in ionized form. Carbon has an ionization potential of 11.3 eV (compared to 13.6 eV for hydrogen), implying that line emission can originate from a variety of phases of the ISM: cold atomic clouds (CNM), diffuse warm neutral and ionized medium (WNM and WIM) and HII regions. Excitation of the [CII] fine structure transition can be via collisions with hydrogen molecules, atoms, and electrons. For example, for the WNM and WIM conditions (Tk = 8000 K; Wolfire et al. 2003) the critical density for excitation of [CII] by H atoms is ~1300 cm-3, and for electrons ~45 cm-3.

thumbnail Fig. 3

[CII] excitation (or spin) temperature, Tex, as a function of total gas density n and kinetic temperature on interstellar gas, Tkin (computed for Tbg = 2.726 K). The upper (lower) panel is dedicated to optically thin (thick) medium. The black solid and dot-dash lines correspond to Tex = 16 K and Tex = 21 K, that are the CMB temperature at z = 5.0 and z = 6.7, respectively. At these redshifts, [CII] emission is suppressed for kinetic temperature below these lines, due to the CMB. For optically thin emission, this suppression affects mostly the cold neutral medium (Tkin ~ 50–120 K, n ~ 20–200 cm-3).

Observationally, it is tremendously difficult to separate the contribution of [CII] emission from all different components. Analysis of [CII] observations is also complicated by the fact that it is difficult to determine the optical depth of the line (e.g. Neri et al. 2014). In the ISM of our Galaxy, because of the density contrast between the CNM and WNM, the [CII] emission associated with the WNM is expected to be a factor of ~20 weaker than that associated with the CNM for a given HI column density (Pineda et al. 2013). In the Galactic plane, Pineda et al. (2014) estimate that 80% of the [CII] comes from atomic and molecular regions, and 20% from ionized gas. In local star-forming galaxies, Croxall et al. (2017) show that 60–80% of [CII] emission originates from neutral gas. This fraction has a weak dependence on the dust temperature and surface density of star formation, and a stronger dependence on the gas-phase metallicity. For metallicities corresponding to the bulk of our galaxies at high redshift (see Fig. 6), the fraction of [CII] emission originating in the neutral phase approaches 90%. At higher redshift, in the interacting system BR1202-0725 at z = 4.7, while [CII] emission arises primarily in the neutral gas for the sub-millimetre galaxy and the quasar, [CII] emission seems to be associated with the ionized medium (H II regions) for one Lyman-α emitter of the system (Decarli et al. 2014). Studying 20 dusty star-forming galaxies from the SPT sample at 2.1 <z< 5.7, Gullberg et al. (2015) found that [CII] emission is consistent with PDRs. Similarly, Stacey et al. (2010) found that the bulk of the [CII] emission line (70%) is originating from PDRs in twelve z ~ 1–2 galaxies.

Theoretically, Olsen et al. (2015) computed the [CII] emission from cosmological smoothed particle hydrodynamics simulations in seven main sequence galaxies at z = 2 and found the ionized gas to have a negligible contribution (<3%). Most of [CII] emission (≳70%) originates from the molecular gas phase in the central ≤1 kpc of their galaxies, whereas the atomic/PDR gas dominates (>90%) further out (>2 kpc). In two zoom-in high-resolution (30 pc) simulations of prototypical M⋆ ~ 1010M⊙ galaxies at z = 6, representative of typical lyman break galaxies at this redshift, 95% of [CII] emission comes from dense gas located in the H2 disc (Pallottini et al. 2017b,a). In their simulations of galaxy formation during the EoR, Katz et al. (2016) found that the majority of [CII] mass is associated with cold neutral clumps and that the [CII] emission (although not computed) is likely to originate in cold, neutral gas, or in PDRs close to young stars.

Thus, it is reasonable to assume that [CII] at high redshift originates mainly from the CNM and PDRs.

However, at the redshifts of interest, one has to consider that the CMB represents a strong background against which the line flux is detected. Indeed, the fraction of the intrinsic line flux observed against the CMB radiation approaches to zero when the excitation temperature (Tex) is close to the CMB temperature. One has thus to check in which physical conditions the [CII] line is being attenuated. To get a first hint, we follow Goldsmith et al. (2012) to compute the [CII] excitation temperature and transpose part of their results to the case of distant unresolved galaxies. The computation is detailed in Appendix A. The deexcitation collision rate coefficients (valid for the range of temperature probed here) are extracted from Barinovs et al. (2005), Goldsmith et al. (2012), Wiesenfeld & Goldsmith (2014): Unlike the situation for molecular tracers, for [CII] emission we do not have the possibility of multiple transitions and many isotopologues to allow determination of the volume density, temperature, and optical depth, and thus obtain a reasonably accurate determination of the column density. While Gullberg et al. (2015) found low to moderate [CII] optical depth, τ[CII] ≲ 1, in a sample of lensed dusty star forming galaxies covering the redshift range z = 2.1–5.7, optically thick [CII] emission was proposed by Neri et al. (2014) for the high-z sub-millimetre source HDF850.1. As the situation is not clear, we use a range of N(C+) and δv such as to cover the optically thin and thick case (Eqs. (A.11) and (A.12)). We consider Tbg = TCMB(z = 0) and calculate the excitation (or spin) temperature, Tex, of the [CII] transition. We show in Fig. 3 the relation between the total gas density n, the kinetic temperature Tkin and Tex for the optically thin and thick cases. We also show the curve corresponding to the CMB temperature at z = 5 and z = 6.7. In any case, [CII] emission from the warm (≃104 K), low density (≲0.1 cm-3) component of the ISM is suppressed at high redshift by the CMB. For the optically thin case, [CII] emission from gas density n ≲ 100 cm-3 (the CNM) will be mostly completely attenuated for z> 6.5 (as also noticed by Vallini et al. 2015), the CMB effect becoming negligible only for galaxies at z ≤ 4.5. In that case, only [CII] from PDRs will reach the instrument. In the optically thick case, only the very cold, low density neutral medium will be affected by CMB attenuation.

Based on these arguments and observational and theoretical constraints detailed above, we can assume that [CII] emission in high-z galaxies arises predominantly from PDRs.

We caution the readers that the above simple calculations do not account for the temperature structure of the clouds. Indeed, as shown in the next section (Sect. 3.2), large CMB attenuation is seen at higher cloud depths, where the temperatures are below ~100 K. At the deepest point of the cloud, where the temperature is only ~50 K, radiative and collisional excitation rates are comparable, but deexcitations are primarily spontaneous. Thus the Aul term dominates the other terms by a few dex in Eq. (A.9).

PDRs are well-studied structures with intense characteristic emissions. Theoretical models addressing the structure of PDRs have been available for approximately 40 yr (Hollenbach et al. 1971; Jura 1974; Glassgold & Langer 1975; Black & Dalgarno 1977). The PDR gas mass fraction in star forming galaxies ranges from a few percent for quiescent systems like the Milky Way up to more than 50% for starburst galaxies like M82 making PDRs important on galactic scales (e.g. Stacey et al. 1991; Malhotra et al. 2001).

Consequently, we assume that the [CII]-line emission originates from PDRs and use CLOUDY to compute its luminosity. Similarly, Popping et al. (2016) only account for the contribution by PDRs to the [CII] luminosity of galaxies when coupling their semi-analytic model of galaxy formation with a radiative transfer code.

3.2. CMB effects on [CII] emission

The structure of a PDR is well established (Hollenbach & Tielens 1999). The outermost layer, which is exposed to the ambient radiation is ionized, and its thickness determined by the ionization parameter. This is followed by a neutral layer, and yet deeper lies a molecular layer. [CII] is present at varying degrees across the PDR. In Fig. 4 we show the structure of a PDR of modest density (log nH = 2.4), exposed to an intense interstellar radiation field, ISRF (ISRF = 3.2 × 103G0, where G0 is the Habing Field, see Sect. 3.3.1). These conditions are close to those predicted in high-z galaxies (see Fig. 6). In this example, the ionized skin of the cloud has a thickness of ~0.1% of the total, and is at a temperature of 10 000 K. The temperature drops sharply in the neutral layer of the cloud, to below 100 K, and the [CII] 158 μm emissivity is reduced. This is illustrated by the shallower slope between 1018.5 and 1020 cm in the top panel of the figure, which presents the emergent intensity in the line (integrating inward). The ionization fraction of [CII] (which is the singly ionized carbon to total carbon mass fraction) drops to below 1% in the molecular core of the cloud, and the intensity flattens at these depths.

Apart from the ISRF, the [CII] emission is influenced by the CMB at high redshift, which peaks at ~130–210 μm at redshifts 4–7. The photon occupation number (Eq. (A.8)) of the CMB dominates the corresponding number in the ISRF for intensities less than ~103–104G0. In our CLOUDY models, both fields are isotropic and subject to removal when corrected line intensities are computed.

Table 1

The ten metallicities considered in our grids of models and the abundances of the five main elements of the ISM we are including in CLOUDY.

On account of the high temperature, the level populations of the transition in the ionized layer are set primarily by collisions, leading to a very small correction for isotropic radiation. The insignificance of radiation in this layer is illustrated by the coincidence of the line intensity at z = 4 and 7 as a function of depth, shown in the top panel of Fig. 4. As the temperature drops in the neutral zone, radiative pumping becomes more important, as does the correction for isotropic radiation. The correction is more important at z = 7, because of the higher photon occupation number. The net effect is that the emergent line intensity corrected for isotropic radiation does not vary significantly between z = 4 and 7.

thumbnail Fig. 4

Example [CII] emission (top) and kinetic temperature (bottom) computed by CLOUDY. Two redshifts are shown (blue and red lines). On the top figure one can see the effect of CMB attenuation on the [CII] line luminosity. In this example, the PDR is homogeneous with log nH = 2.4, and exposed to a radiation field with log ISRF = 3.5, in addition to the CMB radiation.

3.3. [CII] emission from photon dominated regions

3.3.1. The equivalent PDR model

For each galaxy in our SAM, we need to define an equivalent PDR characterised by three parameters: the mean hydrogen density (nH), gas metallicity (Zg), and interstellar radiation field (ISRF).

The mean Hydrogen density (nH) is computed from the mean hydrogen surface density (ΣH) and disc scale height (hd) following . The average disc scale height (hd) is derived at half-mass radius (which is 1.68rd where rd is the exponential disc radius) by assuming a vertical equilibrium in the disc between the gas and stars, and an homogeneous gas velocity dispersion depending on physical properties of galaxies, which is ~15–25 km s-1 (Cousin et al., in prep.). In our simulated galaxies, hd ~ 100 pc at z = 5. Hydrogen densities for the PDRs are computed using a height 5 times smaller, to take into account the fact that the gas is more concentrated to the equatorial plane. Taking hd/ 10 or hd/ 2.5, rather than hd/ 5, modifies the [CII] luminosities by less than 0.1 dex. The distribution of the PDR scale height hpdr is shown in Fig. 5.

ΣH is computed inside a critical radius rc which sets the limit outside of which the gas is not dense enough to form stars (we consider ΣH(r>rc) < 50 M⊙ pc2). By assuming an exponential gaseous disc, the average hydrogen surface density is given by: (6)rd and MH are the exponential disc radius and mass of hydrogen in the disc, respectively.

Our equivalent PDR model depends also on gas metallicity. Oxygen is the most abundant element formed in stars. It is therefore commonly used as a tracer of the gas-phase metallicity. We define the gas metallicity Zg as the number of oxygen atoms to hydrogen atoms with a logarithmic scale: Zg = 12 + log (O/H). We adopt Z⊙ = 8.94 (Karakas 2010). Correspondance between Zg and metals mass fraction is given in Table 1. We assume that the gas is homogeneously distributed in a given galaxy and consequently in a given equivalent PDR. We thus consider average metallicities.

Finally we need to compute the ISRF produced by young stars. It is is defined as the flux of stellar radiation integrated between 6 and 13.6 eV. We assume a mean distance between gas and stars of D = 50 pc. The exact choice of this mean distance has a very small impact on the predicted L[CII]. As commonly used in the literature, the ISRF is normalised in units of the Habing Field (Habing 1968), G0 = 1, which corresponds to f0 = 1.6 × 10-3 erg s-1 cm-2.

thumbnail Fig. 5

Physical sizes of effective PDRs in our simulated galaxies at z = 4.7. Orange and blue histograms are associated to equivalent PDR scale-height (heq = hpdr) and equivalent PDR radius (req), respectively.

3.3.2. CLOUDY: model grids, parameters and outputs

thumbnail Fig. 6

[CII] luminosities predicted by CLOUDY (for a CMB temperature at z = 5). Each panel is dedicated to a given metallicity bin and shows the [CII] luminosity per unit area as a function of both ISRF and hydrogen density. In each panel we show the location of G.A.S. galaxies (extracted from the snapshot at z = 4.7). Rectangles represent the 0.15 and 0.85 quantiles of the galaxy distribution. The black point marks the median, with its coordinates given in the top left corners. We also give the number of galaxies Nobj present in the considered metallicity bin.

Predictions of [CII] emission are computed with the plasma simulation code CLOUDY (Ferland et al. 2013). We use the C17 version of the code (Ferland et al. 2017) as it incorporates a diminution factor due to an external isotropic radiation field (both the CMB and the ISRF, in our case) in its line intensity estimates. This factor was derived as an extended radiative transfer theorem (Chatzikos et al. 2013), and applied to predictions for hyperfine structure line intensities in Chatzikos et al. (2014).

We have built grids of models based on 560 distinct model parameters. A given grid is divided in ten bins of metallicity Zg. In each metallicity bins, ISRF and hydrogen density are sampled following log(ISRF) = [ −0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5] and log(nH) = [ 0, 1, 2, 3, 4, 5, 6].

In a given metallicity bin, abundances of the five main ISM elements are fixed (helium, carbon, nitrogen, oxygen and iron, see Table 1). They are fixed to the median abundances of all simulated galaxies that have gas metallicity in the given bin. Abundances of all other elements are set to zero (we checked that this had no impact on the [CII] luminosity). A galaxy, in a given metallicity bin, will have the median metallicity of the bin, but individual values for hydrogen density and ISRF. The CLOUDY grids are interpolated to find the corresponding [CII] luminosity and the carbon fraction f[CII].

PDRs are modelled assuming a plan-parallel geometry. The shape of the ISRF is that of Black (1987), as given through the CLOUDY option “Table ISM”. Cosmic rays background is fixed to the fiducial value of 2 × 10-16 s-1. For a given PDR model the computation is stopped at AV = 10 and gas temperature can decrease until T = 10 K.

We have generated five grids of models using five background temperatures (Tbg) corresponding to CMB temperatures at z = 4, 5, 6, and 7. We have also computed the grids for both intrinsic and emergent [CII] luminosities. No significant differences have been observed between these two luminosities. We can therefore assumed that [CII] emission is weakly affected by extinction in our galaxies.

For each model associated to a given set of parameters (ISRF−nH−Zg) we extract from CLOUDY the [CII] luminosity per unit area, l[CII ] (in L⊙ cm-2 sr-1), and the [CII] column density, N[ CII ]. We then compute the carbon fraction (in number) in [CII], f[ CII], by computing the ratio between the [CII] column density and the sum of the column density of all species containing carbon atom(s).

From the [CII] luminosity and column density (l[CII ], N[ CII]) we then define the equivalent surface of the PDR as: (7)where Mc is the carbon mass in the galaxy and mc is the mass of individual carbon atom. This formulation implies that we have an uniform [CII] column density in the PDR.

Combined with the [CII] luminosity per unit area (l[CII]), this PDR equivalent surface leads to the following [CII] luminosity, (8)We show on Fig. 5 the distribution of PDR sizes in our simulated galaxies, which is computed using . The median size is around 450 pc. Sizes range from 345 to 745 pc for 70% of the objects. These sizes are in line with the PDR region of M82 (ranging from 300 pc for Joy et al. 1987 to 600 pc Carlstrom & Kronberg 1991). They are in general smaller by a factor ~2 than the estimated sizes of the lensed DSFGs found by SPT and covering the redshift range z = 2.1–5.7 (see Table 2 of Gullberg et al. 2015). This is not surprising as those SPT DSFGs have de-magnified far-IR luminosities (42 <λ< 500 μm) of LFIR ~ 5 × 1012L⊙, thus mean SFR much larger than the average population (see Fig. 7).

thumbnail Fig. 7

L[CII]–SFR relation. Predictions from our model are shown for a set of redshifts from z = 4 to z = 7.6. In each panel the whole sample of G.A.S. galaxies is shown in grey scale. The average relation is plotted with a solid black line. The black dashed line shows the relation given in Eq. (10). Yellow to red coloured points mark the gas metallicity of a randomly selected sample of simulated galaxies (note that the observed tendency of high-metallicity galaxies to fall either above or below the mean trend, that is to make an envelope, is only a trick of the eye caused by the plotting; galaxies with high metallicities (Zg> 8.8) are spread over the whole area, with a higher density of objects at high SFR). Our predictions are compared to a large sample of observational data that are detailed in Table B.1. Amplification corrections on luminosity and SFR, when available, are applied. For dusty star forming galaxies, SFR are converted directly from LIR using the Kennicutt (1998) conversion factor assuming a Chabrier (2003) IMF where SFR (M⊙ yr-1) = 1.0 × 10-10LIR(L⊙). The blue solid line shows the De Looze et al. (2014) relation for the local dwarf galaxy sample.

thumbnail Fig. 8

Same as Fig. 7 but with the SFR determined from the UV (observed, so attenuated) and IR emission, following Kennicutt (1998) standard conversions between LUV and LIR and SFR. This mimics what is being doing when computing the SFR from the UV and IR emission of galaxies. Taking the “observed” SFR (this figure) rather than the “true” SFR from the model (Fig. 7) decreases the dispersion and removes a large fraction of the outliers.

3.3.3. Effect of metallicity, ISRF and density on [CII] emission

We show on Fig. 6 the [CII] luminosity variation as a function of ISRF, hydrogen density and gas metallicity. The grids of models are shown for a background temperature corresponding to the CMB at z = 5. We also show the location of our simulated galaxies extracted from the SAM at z = 4.7. The majority is found in regions where ISRF and hydrogen densities are in the ranges (log ISRF, log nH) = ( [ 3.0, 4.5], [2.0, 3.5]). These ranges of ISRF and hydrogen densities are mainly associated with low gas metallicities (>98% of galaxies with Zg ≤ Z⊙). In the small fraction of galaxies (~2%) with higher gas metallicities (Zg>Z⊙), the physical conditions are different: both the radiation field and hydrogen density are higher, with (log ISRF, log nH) = ( [ 3.5, 5.5], [3.0, 6.0]). These galaxies have discs that are smaller (rd = 0.1 ± 0.2 kpc vs. rd = 0.4 ± 0.2 kpc) and flatter (hd = 0.01 ± 0.04 kpc vs. hd = 0.12 ± 0.10 kpc). Smaller sizes and scale heights lead to a higher gas density. The star formation is thus higher than in the other discs of the sample (SFR = 156 ± 441 M⊙ yr-1 vs. SFR = 5 ± 56 M⊙ yr-1). The strong ISRF associated to the high-metallicity galaxies is therefore explained by the high star formation activity. We finally note that the high gas metallicity of these galaxies is also associated to a high stellar metallicity (Z⋆ = 0.79 ± 0.48 Z⊙ vs. Z⋆ = 0.05 ± 0.12 Z⊙).

Comparing the SPT data with PDR model tracks from Kaufman et al. (1999), Gullberg et al. (2015) obtained a rough estimate of the radiation field strength and gas density of 103<G0< 104 and 102<n< 105 cm-3 for the z> 4 objects, which is in line with our model.

4. L[CII] – SFR relation

The [CII] transition has great potential as a SFR tracer at high redshift. In this section we examine the correlation between L[CII] and SFR at 4 <z< 8 obtained from our model.

De Looze et al. (2014) analyze the applicability of the [CII] line to reliably trace the SFR in a sample of low-metallicity dwarf galaxies from the Herschel Dwarf Galaxy Survey and, furthermore, extend the analysis to a broad sample of galaxies of various types and metallicities in the literature (see also Herrera-Camus et al. 2015; Sargsyan et al. 2014). They found that the L[CII ]–SFR relation has a quite high dispersion, with 1σ = 0.38 dex. Including all the samples from the literature (ex: AGNs, ULIRGS, high-z galaxies) the dispersion increases to 0.42 dex. The scatter in the L[ CII ]–SFR relation increases towards low metallicities, warm dust temperatures, and large filling factors of diffuse, highly ionized gas. At high redshift (z ≃ 7) and using numerical simulations, Vallini et al. (2015) find that the L[ CII ]–SFR relation holds (and is consistent with observations of local dwarf galaxies), with eventual displacements due to extremely low metallicities or a modified Kennicutt-Schmidt relation. The results from their models (obtained assuming a constant metallicity and ΣSFR ∝ ΣH2) are well described by the following best-fitting formula: (9)where L[ CII] is expressed in solar units, and the SFR in M⊙ yr-1.

We show on Fig. 7 the L[CII ]–SFR relation derived from the coupling between G.A.S. and CLOUDY. The predictions are compared to a large set of observational measurements mostly based on UV- or sub-millimetre-selected galaxies where SFRs are either derived from UV flux, deduced from SED-fitting analysis, or computed from LIR. All observational data are compiled in Table B.11. In Fig. 7, we also compare our predictions with the local L[ CII]–SFR relation measured by De Looze et al. (2014) for local dwarf galaxies. Most of these galaxies have metallicities comparable with the bulk of our simulated galaxies at high redshift.

As expected, at each redshift, there is a relation between SFR and the [CII] luminosity, albeit with a very large scatter (0.62 dex at z = 4.0 and 0.51 dex at z = 7.6). To investigate the origin of the scatter, and which one of the three parameters (ISRF, nH or Zg) contributes the most, we compute for each galaxy its [CII] luminosity, fixing two parameters to their median value while keeping the third one to its original value. We find that ISRF, nH, and Zg contribute roughly equally to the scatter.

As shown on Fig. 7, predictions are in remarkably good agreement with the majority of observational data points. A source of dispersion in the observed L[CII ]–SFR relation may come from the fact that the [CII] emission may not overlap with the bulk of UV emission that is used to determine the SFR (Maiolino et al. 2015; Carniani et al. 2017). Nevertheless, the current observations seem to be less scattered than the model. This may be explained by the different timescales that are implicitly assumed when measuring the SFR in galaxies. Our SAM is based on the progressive structuring of the diffuse accreted gas, following an inertial cascade which depends on the fraction of gas already structured in the galaxies (i.e. the gas mass fraction in giant star-forming clouds). This leads to brief and intense star-forming episodes, separated by phases in which the gas starts again its cascade structure (Cousin et al., in prep.). Most of these star-forming episodes occurs on timescales shorter than those assumed when converting UV and far-IR luminosities in SFR using Kennicutt (1998) relations. To investigate the impact on timescales used to determine the SFR, we compute the SFR of galaxies in our model using the UV and IR luminosities, following what is being done from the observations (SFR = SFRUV−obs + SFRIR using Kennicutt (1998) relations to convert luminosities into SFR). We show on Fig. 8 the L[ CII]–SFR relation using this “observed”   SFR. We see that the scatter in the relation decreases; we also see that using the “observed”   SFR removes a large fraction of outliers. This clearly illustrates the importance of timescales when using instantaneous quantities as SFR, and shows that using an average conversion that is the same for all galaxies smoothes the variations.

Two kinds of sample are not falling into our L[CII ]–SFR relation: (i) galaxies with very high SFR and bright [CII] emission (e.g. Gullberg et al. 2015; Oteo et al. 2016; Strandet et al. 2017). They correspond to SFR and L[ CII ] that are not probed in our simulation. Indeed, such objects are rare and are not present in our simulated volume; (ii) galaxies with a strong [CII]-excess emission as compared to their SFR, as tentatively detected in the blind ASPECS survey (Aravena et al. 2016). These galaxies cannot be produced using our assumption that [CII] emission only arises from PDR. The galaxy named “CR7” at z = 6.6 (Matthee et al. 2017) is also a particular case, as it lies close to the average relation but in a region where our simulated galaxies have IR luminosities higher than the upper limit LIR< 3.14 × 1010 L⊙ reported by Matthee et al. (2017). The minimum LIR of our simulated galaxies around the location of CR7 in the L[ CII]–SFR diagram is 4.2 × 1010 L⊙.

thumbnail Fig. 9

L[CII]–SFR relation for galaxies that follow the selection of the simulated galaxy sample by Olsen et al. (2017) at z ~ 6. Our 575 selected galaxies are shown with coloured points (with colour scale reflecting the gas metallicity). Olsen et al. (2017) sample is shown with grey points. The black and the grey solid lines show the mean trend of our whole galaxy sample (Eq. (10)) and of Olsen et al. (2017) sample (their Eq. (6)), respectively. The black dashed line marks the mean trend of our galaxy selection.

We do not observe a strong variation of the mean L[CII ]–SFR relation with redshift. The average trend can be represented by the following law: (10)valid for all redshifts explored here. According to our ranges of available values of SFR and L[ CII ] this average relation is limited to L[ CII]> 107 L⊙ and SFR< 1000 M⊙ yr-1.

We see from Fig. 7 that the De Looze et al. (2014)L[CII ]–SFR relation for the local dwarf galaxies does not really apply to our simulated galaxies but at low [CII] luminosities (LCII ~ 107 L⊙) and z = 4. This mismatch cannot be accounted for only by CMB effects. We show in Fig. C.1, the L[ CII ]–SFR relation when both heating and attenuation by the CMB are ignored. While the L[ CII]–SFR relation becomes more compatible with the local dwarf galaxy sample relation up to z ≃ 6, we still predict a shallower slope at higher redshift.

We show in Fig. 9 the L[CII ]–SFR relation associated to a sub-sample of our simulated galaxies. Following Olsen et al. (2017), this sub-sample of 575 objects has been extracted at z = 5.9 based on three criteria: i) stellar mass between 0.7 and 8 × 109M⊙; ii) SFR between 3 et 23 M⊙ yr-1; and iii) average gas metallicities between 0.15 and 0.45. Compared to the mean trend of the whole galaxy sample at this redshift (Eq. (10)), this sub-sample is biased towards lower L[ CII ] luminosities (by 0.1 dex to 0.5 dex as LCII increases). The L[ CII ]–SFR relation based on our simulated galaxy sub-sample and Olsen et al. (2017) sample are in very good agreement, even if our galaxy sample shows a higher dispersion. This difference could be explained by the different number of objects in the two samples (575 vs. 30). Even if Olsen et al. (2017) galaxies lie mostly in the middle our selection, the mean trend of our sub-sample is shifted to higher L[ CII] by a factor of 0.15 dex. In the light of the very different approaches used between the two models (Olsen et al. 2017 modelled the multi-phased ISM using numerical simulations), the agreement between the two samples is noteworthy. A still better agreement would be obtained by reducing in Olsen et al. (2017) simulations the dust-to-metals ratio by a factor of 2 (consistent with the fact that dust production is less efficient at low metallicities), or by decreasing the slope of the giant molecular cloud mass spectrum from 1.8 (which corresponds to the Galactic mass spectrum) to 1.5 (see the discussion in Sect. 5.1 of Olsen et al. 2017).

thumbnail Fig. 10

Distribution of the simulated galaxies in [CII] luminosities and metallicities at z = 5.9. Coloured points mark the SFR of a randomly selected sample. Blue squares are from the local dwarf galaxy sample (Cormier et al. 2015).

To study the variation of the [CII] luminosities with gas metallicities, we plot on Fig. 10 the distribution of our simulated galaxies in the L[CII] – Zg plane. There is a broad correlation; galaxies with higher [CII] luminosities tend to have higher metallicities. However, the scatter is quite large, with 0.84 and 0.72 dex and z = 5 and 7, respectively. We also show the data points from the local dwarf galaxy sample (Cormier et al. 2015). They are well sampled by the distribution of our simulated galaxies.

Finally, we compare our prediction with Vallini et al. (2015) model (Eq. (9)) in Fig. 11. For this comparison and to be consistent with Vallini et al. (2015), we consider our model at z = 5.9 and with CMB effects. We observe a systematic trend with a decreasing L[CII ]/ ratio with both SFR and Zg. The bulk of our galaxies has higher L[ CII] than that predicted by Eq. (9) (by factors of about 1.5–5). This equation holds when the molecular mass is scaled with the SFR (the Kennicutt-Schmidt relation), adopting N = 1. As discussed in Vallini et al. (2015), the range in power law index (N) depends on a variety of factors, among which the most important ones are the observed angular scales, and the calibration of SFRs. As shown in Fig. 11, we have a better agreement considering a fit of Vallini et al. (2015) model with N = 2: (11)This is expected as the Kennicutt-Schmidt relation predicted by G.A.S. has 1.4 ≤ N ≤ 2.

thumbnail Fig. 11

Ratio of our predicted [CII] luminosity at z ~ 6 (L[CII]) and that predicted by Vallini et al. (2015) model (). The top panel considers from Eq. (9), and thus N = 1, while the bottom panel considers from Eq. (11), and thus N = 2. N is the power law index of the Kennicutt-Schmidt relation, .

thumbnail Fig. 12

L[CII ]/LIR vs. LIR for our sample of simulated galaxies (grey shaded areas) at z = 4.7. Local galaxies (orange circles) are the GOALS luminous infrared galaxy sample (Díaz-Santos et al. 2017); a mean correction has been applied to their L[ CII] to mimic the effects of CMB heating and attenuation. High-resdhift galaxies (blue circles) are extracted from Table  B.1. We give on the bottom of the figure the mean densities, ISRF, and metallicities (in log), for bins of 0.5 dex in log  LIR. The [CII] deficit naturally arises in our model. It is well correlated with the intensity of the ISRF.

5. [CII] deficit

In the early days of [CII] observations of low-redshift galaxies with the Infrared Space Observatory, it was observed that very luminous infrared galaxies such as ULIRGs appear to have a deficit in [CII] emission compared to their FIR luminosities (Luhman et al. 1998; Malhotra et al. 2001). This deficit has been later confirmed with Herschel and extended to lower infrared luminosities, LIRGs. For example, Díaz-Santos et al. (2013) find that LIRGs show a tight correlation of [CII]/FIR with infrared luminosity, with a strong negative trend spanning from ~10-2 to 10-4, as the infrared luminosity increases. The result from high-redshift objects is more mixed, with a large scatter (2 orders of magnitude) at high luminosities. Different explanations for this measured decline have been proposed (e.g. Casey et al. 2014), including the compactness of the starburst, the AGN activity, optically thick [CII] emission, varying IMF or [CII] saturation at high temperature.

As shown on Fig. 12, the [CII] deficit naturally arises in our model, with a decrease of L[CII ]/LIR by about 2 orders of magnitudes from LIR = 109 to 1012 L⊙. The dispersion in the L[ CII ]/LIR increases with LIR. The large dispersion at high LIR is also observed for high-redshift objects (e.g. Fig. 4 of Gullberg et al. 2015). Compared to the GOALS sample of local galaxies, the L[ CII ]/LIR decrease is stronger in the model and simulated galaxies have on average a lower L[ CII]/LIR. This might be due to selection effects (often linked to a limited sensitivity in the observations).

We investigate the origin of the deficit using our model parameters. First, we compute the [CII] transition upper level loading to test the hypothesis of Muñoz & Oh (2016), in which the [CII] deficit observed in the highest IR surface-brightness systems is a natural consequence of saturating the upper fine-structure transition state at gas temperatures above 91 K. We find that from z = 4 to 7, the transition upper level loading is not saturated in the region where the bulk of the [CII] intensity is emitted (see Fig. 4), and that 0.01 <nu/nl< 0.05. The crucial difference between our model and the analytical work of Muñoz & Oh (2016) is that we have strong ISRF for our galaxies. Muñoz & Oh (2016) ignore the effects of the local (isotropic) radiation field, under the assumption of densities in excess of the critical density for that transition. Therefore, this saturation effect cannot be responsible for the [CII] deficit observed in our model. We searched for correlations of the deficit with the different parameters of our model and found that it is strongly correlated with the intensity of the ISRF (see Fig. 12). This is consistent with the analysis of Luhman et al. (2003) which suggests that a high ISRF incident on a moderate density PDR could explain the deficit in their observed ULIRGs. We extend this analysis to lower luminosities, and show that the deficit still holds at very high redshift. For 1010<LIR< 3 × 1011 L⊙, we have a weak correlation of the deficit with the metallicity, with slightly increased metallicity associated with deeper deficits, as observed in Smith et al. (2017) but for higher metallicities (Zg between 8.54 and 8.86).

6. [CII] luminosity function

thumbnail Fig. 13

[CII] luminosity function predicted by the G.A.S.+CLOUDY model from z = 4.0 to z = 7.6. The blue solid curve shows the prediction that accounts for the attenuation of [CII] emission due to the CMB. The blue dotted line would be the luminosity function ignoring the attenuation. At z ≃ 5, we show the observational constraints from Capak et al. (2015). At z ≃ 6, the black squares indicate the observational results of Yamaguchi et al. (2017). We also add the local [CII] luminosity function published by Hemmati et al. (2017; orange dotted line) and model predictions of Popping et al. (2016; grey solid line).

Figure 13 shows our [CII] luminosity function predicted at 4.0 ≲ z ≲ 8. We present two different predictions, with and without CMB attenuation. We see a systematic deviation between the two, which is almost constant with redshift (as expected, see Sect. 3.2). The attenuation induced by the CMB increases slowly with the [CII] luminosity, from 25% at L[CII] ~ 107 L⊙ to 35% at L[CII] ~ 1010 L⊙. This trend is similar at all redshifts.

We also show on Fig. 13Popping et al. (2016) luminosity function predictions, that are also based on a semi-analytical model. Compared to Popping et al. (2016), we predict a smaller (larger) number of [CII]-emitting galaxies in the faint (bright)-end part, with a crossing point at L[CII] ≃ 2 × 108 L⊙.

We can also compare our [CII] luminosity function with that obtained by using the SFR function from Smit et al. (2012) and our mean L[CII ]–SFR relation. We found that such a combination overestimates the [CII] luminosity function, by factor ~6 at z = 4 for L[ CII] = 108 L⊙ for example.

6.1. The functional form of the luminosity function

Our predicted [CII] luminosity function has a power law shape for the whole range of L[CII] probed in our simulation. This shape is quite different from the [CII] luminosity function measured at z = 0 (Hemmati et al. 2017), which agrees well with the form of the IR luminosity function. This IR luminosity function is better fitted either by a double power law (Magnelli et al. 2011), following, (12)or alternatively by a double-exponential function (Saunders et al. 1990; Caputi et al. 2007; Gruppioni et al. 2013),

which is a modified-Schechter function behaving as a power law for L ≪ L⋆ and as a Gaussian in log L for L ≫ L⋆: (13)In these equations, L⋆ is the characteristic luminosity where the transition between the faint and bright regimes occurs, and Φ⋆ is the normalization factor2. Unfortunately, the IR luminosity function has not been measured at z ≥ 4. At lower redshift, it is found that log  increases with redshift, from 10.48 (z = 0) to 12.35 (z ~ 2) in Magnelli et al. (2013), assuming a double power law, and from 10.12 (z = 0) to 11.9 (z ~ 4) in Gruppioni et al. (2013) and 11.40 (z ~ 1) and 11.80 (z ~ 2) in Caputi et al. (2007), assuming a double-exponential function. Thus, a typical  L⊙ is expected at high redshift. This IR luminosity converts to SFR = 100 M⊙ yr-1 (using Kennicutt (1998) and assuming a Chabrier (2003) IMF). Using our L[CII]–SFR relation (Eq. (10)), we obtain log  to 9.1, from z = 7 to 4, respectively. These characteristic luminosities are quite high and difficult to probe with our model. They fall in a regime where we have less than 50 objects in our simulation.

A power law shape is not completely unexpected at these very high redshifts. A single power law provides an equally good fit to the UV luminosity function at z = 8, while at z = 6 and 7, an exponential cutoff at the bright end is moderately preferred (Finkelstein et al. 2015). For the stellar mass function (SMF), the knee is sharpened as time goes by, and a progressive flattening of the low-mass end is observed from z ~ 6 to zero (Davidzon et al. 2017). At z> 4.5 the SMF is best fit by a single power law function with a cut-off at 3 × 1011M⊙. Measurements of the SMF extend down to ~1010M⊙. In this range of masses, the [CII] luminosities range from ~2 × 108 to 5 × 109L⊙. Accordingly, a cut-off may be expected in the [CII] luminosity function but outside of the range of luminosities probed by our simulation. This may explain why we do not see a break in our [CII] luminosity function.

We compare on Fig. 14 the observed UV luminosity function at z = 5 with that predicted from the [CII] luminosity function, applying a UV to [CII] luminosity ratio. For that comparison, it is primordial to correct the observed UV luminosity function for attenuation as dust strongly affects the shape of the observed UV luminosity function. We use our model to derive a mean attenuation per UV magnitude bin, and correct the observed UV data points. We have a quite good agreement between the corrected UV luminosity function and that predicted from [CII] using a constant luminosity ratio, /L[CII ] = 1.6 × 103. This constant ratio is a crude approximation but it is not worth searching for any variation with luminosity given the large and uncertain corrections for dust attenuation. This ratio is much larger than the LUV/L[ CII] ratios of ~100 to 650 obtained for UV-selected galaxies at z = 5 (Barisic et al. 2017). Part of the discrepancy may be attributed to dust attenuation of the observed UV light.

thumbnail Fig. 14

UV luminosity function derived from the [CII] luminosity function at z ~ 5. Small and large symbols are observational measurements without and with extinction correction, respectively. The correction of extinction of observed measurements has been done using G.A.S. (using the two grey cuves of Fig. 2). The solid grey line shows our predicted UV luminosity function from G.A.S., which is in good agreement with the observed data points coming from Bouwens et al. (2015) and Finkelstein et al. (2015). The UV luminosity function predicted from our [CII] luminosity function assuming a fixed [CII] to UV luminosity ratio is shown in blue.

6.2. Redshift evolution

The slope of the power law that fits the luminosity function is not evolving strongly with redshift for 4 ≲ z ≲ 8, and is ≃–1. At such high redshift, the slope of the IR luminosity function for is not known; at lower redshift, it is usually fixed to –0.6 (e.g. Magnelli et al. 2011), but a shallower faint-end slope (α1 = −0.4) has been recently measured at 1.5 <z< 2.5 (Koprowski et al. 2017). At z = 0, the slope of the [CII] luminosity function is equal to –0.42 for , where . The steepening of the faint-end slope of the [CII] luminosity function between z = 0 and z> 4 may reflect the fact that the galaxy population is richer in faint [CII] emitting galaxies with increasing redshift, which is the natural consequence of the hierarchical formation of galaxies. In the UV, the faint-end slopes varies from –1.5 at z = 4 to –2 at z = 7 (Bowler et al. 2015; Finkelstein et al. 2015) but part of the steepening in that case may be explained by a changing impact of dust attenuation with redshift. For [CII], only the CMB is attenuating the luminosity and this does not cause any change in the slope of the power law (see Fig. 13).

At a given [CII] luminosity the density of object decreases with redshift following, (14)valid for redshifts 4.7 ≤ z ≤ 8. The slope of the decrease in density is equal to –0.4; it is close to the slope of –0.31 seen in UV (Finkelstein et al. 2015). The density evolves to higher values by a factor of 20× from z = 7.6 to z = 4.7. At the characteristic luminosity of the local [CII] luminosity function ( L⊙), the number density is 3.8 times lower at z = 4 than at z = 0.

6.3. Comparison with observational constraints

thumbnail Fig. 15

Cumulative [CII] luminosity functions predicted at z ≃ 4 and z ≃ 6 (thick continuous blue and orange lines, respectively, with the thin lines representing ), compared to the current observational constraints (data points with error bars). The dot-dashed lines are derived by relaxing the criterium req<rg which concerns only <0.8% of the galaxies. The dashed lines show the [CII] luminosity functions derived from the carbon mass function of our model.

To date, observational constraints are very sparse, with upper limits from Yamaguchi et al. (2017) at z ≃ 6 and estimates at z ≃ 5 derived from Capak et al. (2015) in Hemmati et al. (2017). As explained in Hemmati et al. (2017), the z ≃ 5 estimates are very rough. They are based on observations of nine Lyman-break galaxies in [CII] using ALMA. To see where these measurements sit compared to the luminosity function, Hemmati et al. (2017) measure the volume for each observation using the area and the redshift width of each ALMA pointing and correct the volume using the number density of Lyman-break galaxies. These factors and the low number statistics makes these estimates very sensitive on choice of bins and therefore uncertain. Moreover, there might exist classes of galaxies that are not selected as Lyman-break galaxies at high redshifts but that are bright in [CII] and can contribute to the luminosity function.

As we can see from Fig. 13, upper limits are not giving very stringent constraints and our predicted luminosity function is well below. Our predicted [CII] luminosity function at z = 4.7 is compatible with the two points estimated from Capak et al. (2015; at 1 and 3σ, respectively).

More observational constraints are available on the cumulative [CII] luminosity functions. At z ≃ 4, there exists both lower (Swinbank et al. 2012) and upper (Matsuda et al. 2015) limits. At 6 <z< 8, the ASPECS blind survey gives only upper limits to the bright end of the [CII] luminosity function, as their detections are candidate [CII]-line emitters (Aravena et al. 2016). We also consider their number density assuming that only the brightest candidate is real. Hayatsu et al. (2017) found two [CII] emitter candidates at z = 6.0 and 6.5 in the ALMA 1.1-mm survey of the SSA22 field. They estimate the luminosity function at z = 6.2 from blind detection on the assumption that one of the two unconfirmed lines is [CII] at z ~ 6. We show on Fig. 15 the comparison between our predictions and those constraints. We do not consider the measurements from Miller et al. (2016) as their ALMA sample is biased to fields of extreme objects at z> 6 and cannot be used to directly constrain the field luminosity function. At z ≃ 4, our model falls between the observational upper and lower limits. At z ≃ 6, it is well below the upper limits. It is also 1.3σ below the two estimates from Hayatsu et al. (2017) and Aravena et al. (2016). These estimates are given for L[CII] ≳ 5.4 × 108 L⊙ and ≳9.1 × 108L⊙, where we have few objects in our simulated volume (110 and 60, respectively – i.e. <1.5%). These galaxies are characterized by a high metallicity, Zg ≥ 8.0, a SFR ≥ 100 M⊙/ yr and therefore a strong ISRF. A 1σ agreement with the measurements would need an increase of the number density of such [CII]-emitting galaxies by a factor of 2.5.

To investigate this small discrepancy, we have built a [CII] luminosity function derived from the carbon mass function. We assume a constant carbon mass to [CII] luminosity ratio (R). By using R ≃ 25 L⊙/M⊙, the mass-derived [CII] luminosity function lies very close to that built with our G.A.S.+CLOUDY model, for L[CII ] ≤ 108 L⊙, as shown on Fig. 15. But we clearly see that this simple model predicts an excess of bright objects (for L[ CII] ≳ 1.5 × 108 L⊙). Some galaxies can reach [CII] luminosities ten times higher than in our fiducial model. The predictions of this simple model is then in agreement with the current observational constraints. This result indicates that the carbon content in our simulated galaxies is sufficient to produce an excess on the bright-end of the luminosity function.

As explained in Sect. 3.3 our fiducial model is based on an equivalent PDR structure. For each galaxy, we assumed a metallicity, hydrogen density and ISRF. These three parameters allowed us to compute, for each PDR, a [CII] luminosity per unit of surface area and an effective surface area of emission. In some rare cases (<0.8%) the equivalent radius of the PDR is larger than the radius encompassing the whole mass of the galaxies (rg = 11rd, containing 99.9% of the mass). In these cases, we artificially limited the PDR equivalent radius req to rg (while keeping ISRF, nH and Zg to their original values). This led to a reduction of the surface area of the emission, SPDR (Eq. (7)), and therefore of the [CII] luminosity. If this size criterion is relaxed, the [CII] luminosity function becomes very close to that obtained with our simple model. Galaxies that violate the size criterion are very small (disc size <1 kpc), extremely dense (nH > 105 cm-3) and have high metallicities (Zg > 8.5). These objects (which represent a very small fraction of the sample, <0.8%) are probably not a realistic population of galaxies and cannot account for the difference seen between the predicted and observed cumulative luminosity functions. With our model, we cannot produce much more high L[CII] objects by simply changing the parameters (i.e. ISRF, nH and Zg) in reasonable proportions. Higher [CII] luminosities could be obtained by considering an additional excitation of the [CII] line by other processes, as AGN emission.

7. Conclusions

We have used our semi-analytical model of galaxy formation (G.A.S.) combined with the photoionisation code CLOUDY to compute the [CII] luminosity for a large number of galaxies at z ≥ 4 (~28 000 at z = 5). With such a large statistical sample, we can investigate the dispersion in the L[CII]–SFR relation as well as derive the [CII] luminosity function. Our model takes into account the effects of CMB heating and attenuation that are important at such high redshifts.

We showed that our model is able to reproduce the L[CII ]–SFR relation observed for ~50 star-forming galaxies at z ≥ 4. However, our model does not contain any galaxy with a very strong [CII]-excess emission as compared to their SFR, as found in the blind ASPECS survey (Aravena et al. 2016). More generally, we found that the L[ CII]–SFR relation is very dispersed (0.51 to 0.62 dex from z = 7.6 to z = 4), the large dispersion being due to combined effects of different metallicities, ISRF and gas contents in the simulated high-redshift galaxies. The high dispersion provides an explanation to the upper limits obtained on a number galaxies at z ≥ 6 (e.g. González-López et al. 2014; Ota et al. 2014). We found that the dispersion and the fraction of outliers are reduced when the SFR of galaxies is derived from the UV and IR luminosities, following what is being done from the observations (SFR = SFRUV−obs + SFRIR). This demonstrates the importance of timescales when using instantaneous quantities as SFR (the timescales being shorter in the model than those assumed in the luminosity-SFR conversions), and the effect of using average conversions. CMB attenuation and heating (which becomes important in the cold gas) also contribute to the dispersion, because its effects depend on the properties of each galaxies (e.g., kinetic temperature and density of the gas). It will be very difficult to correct individual [CII] observations from CMB effects because this would require to know the physical properties of the [CII]-emitting gas.

We observed a small evolution of the L[CII ]–SFR relation with redshift, with a decrease of the [CII] luminosity of only ~30% from z = 4 to z = 7.6 at a given SFR. Our L[ CII]–SFR relation at z ≥ 5 is not compatible with the relation for the local dwarf galaxy sample. Finally, we also showed that there is a broad correlation, with a scatter ~0.8 dex, between the [CII] luminosity and gas metallicity.

We found that our model naturally predicts the [CII] deficit, with a decrease of L[CII]/LIR by about 2 orders of magnitudes from LIR = 109 to 1012 L⊙. We investigated the origin of the deficit and found that it is strongly correlated with the intensity of the ISRF.

We then presented the predictions for the [CII] luminosity function for 4 ≤ z ≤ 8 and log  L[CII ] ≥ 7, which is our completeness limit. On the bright end, our simulations contain less than 10 objects with log  LCII higher than 9.9, 9.9, 9.4, 9.2, and 8.8 at z ≃ 4.0, 4.7, 5.9, 6.7, and 7.6, and these values are thus the upper bounds of our predictions. In these ranges of LCII, the luminosity function has a power law shape with α = −1. This may be explained by a redshift evolution characterized by continued positive luminosity evolution, as seen for the IR luminosity function (Koprowski et al. 2017). The characteristic luminosity  L⊙, measured at z = 0, should then increase by a factor of about 5 at z = 4 and recover the local value at z = 7.6. In the mean time, the number density decreases by a factor of 20× from z = 4.7 to z = 7.6. At those redshifts, we have a reasonable agreement between the UV and [CII] luminosity functions considering a constant luminosity ratio, /L[ CII ] = 1.6 × 103, and assuming a correction for attenuation of UV luminosities derived from our model. Finally we compared our predictions with the few observational constraints. We found that our differential luminosity function is in reasonable agreement with the observational estimates, but our cumulative luminosity function is 1.3σ below the estimates at z = 6 and L[ CII] ≃ 5–9 × 108 L⊙. By relaxing a parameter in the model that constrains the size of the effective PDR and that affects only <0.8% of simulated galaxies, we can increase the number density of bright [CII]-emitters and better match the estimates on the cumulative luminosity function at z = 6.

Our model relies on the assumption that the [CII] line is originating exclusively from PDRs, with one effective PDR defined for each galaxy. It does not take into account the whole complexity of the ISM in galaxies, as the structure of giant clouds or inhomogeneous ISRF, that can affect the [CII] luminosities. So, it is remarquable how this simplified model can reproduce the observations at high redshift. Our G.A.S.+CLOUDY predictions are also in good agreement with those obtained from cosmological zoom simulations of galaxies combined with a multi-phased ISM modelling (Olsen et al. 2017). However, the main limitations of all current models is that they miss the contribution from [CII] that can be excited (i) on a large scale by the dissipation of mechanical energy (turbulence and shocks) in the early stages of the building of galaxy discs (Appleton et al. 2013), and (ii) by the AGN.

1

We did not add the upper limits coming from the [CII] search in bright Lyman-alpha emitters (LAEs), but the locations of such galaxies in the L[CII ]–SFR plot are not unexpected. For example, at the SFR of the three LAEs of González-López et al. (2014), most galaxies in our simulation have much fainter L[ CII] than their reported upper limits. And the galaxies that lie above the upper limits are those with high metallicities (thus they are not dust-poor, bright LAEs).

Acknowledgments

We acknowledge financial support from the “Programme National de Cosmologie and Galaxies” (PNCG) funded by CNRS/INSU-IN2P3-INP, CEA and CNES, France, from the ANR under the contract ANR-15-CE31-0017 and from the OCEVU Labex (ANR-11-LABX-0060) and the A*MIDEX project (ANR-11-IDEX-0001-02) funded by the “Investissements d’Avenir” French government programme managed by the ANR. MC acknowledges the support from the CNES. We warmly thank Maryvonne Gérin for enlightening discussions, and Maxime Ruaud for providing us his code of [CII] excitation. MC thank Jacques Le Bourlot, Franck Le Petit and Benjamin Godard for their help in using PDR Meudon code (which has been used for crosschecks). We thank Gary Ferland for pointing us the new C17 version of CLOUDY, which incorporates the treatment of isotropic backgrounds. We thank Lin Yan and Shoubaneh Hemmati for providing us the data from their paper on the local [CII] luminosity function, Tanio Díaz-Santos for providing us the data for the [CII] deficit of local galaxies, and Steven Finkelstein and Rychard Bouwens for providing us the data points of their high-redshift UV luminosity functions. Finally, we thank the anonymous referee for his/her very useful comments and responsiveness.

References

  1. Appleton, P. N., Guillard, P., Boulanger, F., et al. 2013, ApJ, 777, 66 [NASA ADS] [CrossRef] [Google Scholar]
  2. Aravena, M., Decarli, R., Walter, F., et al. 2016, ApJ, 833, 71 [NASA ADS] [CrossRef] [Google Scholar]
  3. Barinovs, Ğ., van Hemert, M. C., Krems, R., & Dalgarno, A. 2005, ApJ, 620, 537 [NASA ADS] [CrossRef] [Google Scholar]
  4. Barisic, I., Faisst, A. L., Capak, P. L., et al. 2017, ApJ, 845, 41 [NASA ADS] [CrossRef] [Google Scholar]
  5. Benson, A. J., & Bower, R. 2011, MNRAS, 410, 2653 [NASA ADS] [CrossRef] [Google Scholar]
  6. Béthermin, M., Daddi, E., Magdis, G., et al. 2015, A&A, 573, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Black, J. H. 1987, in Interstellar Processes, eds. D. J. Hollenbach, & H. A. Thronson, Jr., Astrophy. Space Sci. Lib., 134, 731 [Google Scholar]
  8. Black, J. H., & Dalgarno, A. 1977, ApJS, 34, 405 [NASA ADS] [CrossRef] [Google Scholar]
  9. Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34 [Google Scholar]
  10. Bowler, R. A. A., Dunlop, J. S., McLure, R. J., et al. 2015, MNRAS, 452, 1817 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bradač, M., Garcia-Appadoo, D., Huang, K.-H., et al. 2017, ApJ, 836, L2 [NASA ADS] [CrossRef] [Google Scholar]
  12. Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000 [NASA ADS] [CrossRef] [Google Scholar]
  13. Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682 [NASA ADS] [CrossRef] [Google Scholar]
  14. Capak, P. L., Carilli, C., Jones, G., et al. 2015, Nature, 522, 455 [NASA ADS] [CrossRef] [Google Scholar]
  15. Caputi, K. I., Lagache, G., Yan, L., et al. 2007, ApJ, 660, 97 [NASA ADS] [CrossRef] [Google Scholar]
  16. Caputi, K. I., Ilbert, O., Laigle, C., et al. 2015, ApJ, 810, 73 [NASA ADS] [CrossRef] [Google Scholar]
  17. Carlstrom, J. E., & Kronberg, P. P. 1991, ApJ, 366, 422 [NASA ADS] [CrossRef] [Google Scholar]
  18. Carniani, S., Maiolino, R., Pallottini, A., et al. 2017, A&A, 605, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Casey, C. M., Narayanan, D., & Cooray, A. 2014, Phys. Rep., 541, 45 [Google Scholar]
  20. Chabrier, G. 2003, ApJ, 586, L133 [Google Scholar]
  21. Charlot, S., & Fall, S. M. 2000, ApJ, 539, 718 [Google Scholar]
  22. Chatzikos, M., Ferland, G. J., Williams, R. J. R., Porter, R., & van Hoof, P. A. M. 2013, ApJ, 779, 122 [NASA ADS] [CrossRef] [Google Scholar]
  23. Chatzikos, M., Ferland, G. J., Williams, R. J. R., & Fabian, A. C. 2014, ApJ, 787, 96 [NASA ADS] [CrossRef] [Google Scholar]
  24. Compiègne, M., Verstraete, L., Jones, A., et al. 2011, A&A, 525, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Cooksy, A. L., Saykally, R. J., Brown, J. M., & Evenson, K. M. 1986, ApJ, 309, 828 [NASA ADS] [CrossRef] [Google Scholar]
  26. Cormier, D., Madden, S. C., Lebouteiller, V., et al. 2015, A&A, 578, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Cousin, M., Lagache, G., Bethermin, M., Blaizot, J., & Guiderdoni, B. 2015a, A&A, 575, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Cousin, M., Lagache, G., Bethermin, M., & Guiderdoni, B. 2015b, A&A, 575, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Cousin, M., Buat, V., Boissier, S., et al. 2016, A&A, 589, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Cox, P., Krips, M., Neri, R., et al. 2011, ApJ, 740, 63 [Google Scholar]
  31. Croxall, K., Smith, J. D. T., Pellegrini, E., et al. 2017, ApJ, 845, 96 [NASA ADS] [CrossRef] [Google Scholar]
  32. Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Decarli, R., Walter, F., Carilli, C., et al. 2014, ApJ, 782, L17 [NASA ADS] [CrossRef] [Google Scholar]
  34. Decarli, R., Walter, F., Venemans, B. P., et al. 2017, Nature, 545, 457 [NASA ADS] [CrossRef] [Google Scholar]
  35. De Looze, I., Cormier, D., Lebouteiller, V., et al. 2014, A&A, 568, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Devriendt, J. E. G., Guiderdoni, B., & Sadat, R. 1999, A&A, 350, 381 [NASA ADS] [Google Scholar]
  37. Díaz-Santos, T., Armus, L., Charmandaris, V., et al. 2013, ApJ, 774, 68 [NASA ADS] [CrossRef] [Google Scholar]
  38. Díaz-Santos, T., Armus, L., Charmandaris, V., et al. 2017, ApJ, 846, 32 [Google Scholar]
  39. Dole, H., Lagache, G., Puget, J.-L., et al. 2006, A&A, 451, 417 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
  40. Duncan, K., Conselice, C. J., Mortlock, A., et al. 2014, MNRAS, 444, 2960 [NASA ADS] [CrossRef] [Google Scholar]
  41. Ferland, G. J., Porter, R. L., van Hoof, P. A. M., et al. 2013, Rev. Mexicana Astron. Astrofis., 49, 137 [Google Scholar]
  42. Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, Rev. Mex. Astron. Astrofis, 53, 385 [Google Scholar]
  43. Finkelstein, S. L., Ryan, Jr., R. E., Papovich, C., et al. 2015, ApJ, 810, 71 [NASA ADS] [CrossRef] [Google Scholar]
  44. Giallongo, E., Grazian, A., Fiore, F., et al. 2015, A&A, 578, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Glassgold, A. E., & Langer, W. D. 1975, ApJ, 197, 347 [NASA ADS] [CrossRef] [Google Scholar]
  46. Gnedin, N. Y. 2000, ApJ, 542, 535 [NASA ADS] [CrossRef] [Google Scholar]
  47. Goldsmith, P. F., Langer, W. D., Pineda, J. L., & Velusamy, T. 2012, ApJS, 203, 13 [Google Scholar]
  48. Gong, Y., Cooray, A., Silva, M., et al. 2012, ApJ, 745, 49 [NASA ADS] [CrossRef] [Google Scholar]
  49. González-López, J., Riechers, D. A., Decarli, R., et al. 2014, ApJ, 784, 99 [NASA ADS] [CrossRef] [Google Scholar]
  50. Grazian, A., Fontana, A., Santini, P., et al. 2015, A&A, 575, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 432, 23 [NASA ADS] [CrossRef] [Google Scholar]
  52. Guiderdoni, B., & Rocca-Volmerange, B. 1987, A&A, 186, 1 [NASA ADS] [Google Scholar]
  53. Gullberg, B., De Breuck, C., Vieira, J. D., et al. 2015, MNRAS, 449, 2883 [NASA ADS] [CrossRef] [Google Scholar]
  54. Habing, H. J. 1968, Bull. Astron. Inst. Netherlands, 19, 421 [Google Scholar]
  55. Hayatsu, N. H., Matsuda, Y., Umehata, H., et al. 2017, PASJ, 69, 45 [Google Scholar]
  56. Hemmati, S., Yan, L., Diaz-Santos, T., et al. 2017, ApJ, 834, 36 [NASA ADS] [CrossRef] [Google Scholar]
  57. Herrera-Camus, R., Bolatto, A. D., Wolfire, M. G., et al. 2015, ApJ, 800, 1 [NASA ADS] [CrossRef] [Google Scholar]
  58. Hollenbach, D. J., & Tielens, A. G. G. M. 1999, Rev. Mod. Phys., 71, 173 [Google Scholar]
  59. Hollenbach, D. J., Werner, M. W., & Salpeter, E. E. 1971, ApJ, 163, 165 [NASA ADS] [CrossRef] [Google Scholar]
  60. Inoue, A. K., Tamura, Y., Matsuo, H., et al. 2016, Science, 352, 1559 [NASA ADS] [CrossRef] [Google Scholar]
  61. Joy, M., Lester, D. F., & Harvey, P. M. 1987, ApJ, 319, 314 [NASA ADS] [CrossRef] [Google Scholar]
  62. Jura, M. 1974, ApJ, 191, 375 [NASA ADS] [CrossRef] [Google Scholar]
  63. Karakas, A. I. 2010, MNRAS, 403, 1413 [NASA ADS] [CrossRef] [Google Scholar]
  64. Katz, H., Kimm, T., Sijacki, D., & Haehnelt, M. 2016, ArXiv e-prints [arXiv:1612.01786] [Google Scholar]
  65. Kaufman, M. J., Wolfire, M. G., Hollenbach, D. J., & Luhman, M. L. 1999, ApJ, 527, 795 [Google Scholar]
  66. Kennicutt, Jr., R. C. 1998, ARA&A, 36, 189 [Google Scholar]
  67. Khochfar, S., & Silk, J. 2009, ApJ, 700, L21 [NASA ADS] [CrossRef] [Google Scholar]
  68. Knudsen, K. K., Richard, J., Kneib, J.-P., et al. 2016, MNRAS, 462, L6 [NASA ADS] [CrossRef] [Google Scholar]
  69. Knudsen, K. K., Watson, D., Frayer, D., et al. 2017, MNRAS, 466, 138 [NASA ADS] [CrossRef] [Google Scholar]
  70. Koprowski, M. P., Dunlop, J. S., Michałowski, M. J., et al. 2017, MNRAS, 471, 4155 [NASA ADS] [CrossRef] [Google Scholar]
  71. Kravtsov, A. V., Gnedin, O. Y., & Klypin, A. A. 2004, ApJ, 609, 482 [NASA ADS] [CrossRef] [Google Scholar]
  72. Lagache, G., Puget, J.-L., & Dole, H. 2005, ARA&A, 43, 727 [NASA ADS] [CrossRef] [Google Scholar]
  73. Luhman, M. L., Satyapal, S., Fischer, J., et al. 1998, ApJ, 504, L11 [NASA ADS] [CrossRef] [Google Scholar]
  74. Luhman, M. L., Satyapal, S., Fischer, J., et al. 2003, ApJ, 594, 758 [NASA ADS] [CrossRef] [Google Scholar]
  75. Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [NASA ADS] [CrossRef] [Google Scholar]
  76. Magnelli, B., Elbaz, D., Chary, R. R., et al. 2011, A&A, 528, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Magnelli, B., Popesso, P., Berta, S., et al. 2013, A&A, 553, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Maiolino, R., Carniani, S., Fontana, A., et al. 2015, MNRAS, 452, 54 [Google Scholar]
  79. Malhotra, S., Kaufman, M. J., Hollenbach, D., et al. 2001, ApJ, 561, 766 [Google Scholar]
  80. Matsuda, Y., Nagao, T., Iono, D., et al. 2015, MNRAS, 451, 1141 [NASA ADS] [CrossRef] [Google Scholar]
  81. Matthee, J., Sobral, D., Boone, F., et al. 2017, ApJ, 851, 145 [NASA ADS] [CrossRef] [Google Scholar]
  82. Miller, T. B., Chapman, S. C., Hayward, C. C., et al. 2016, MNRAS, submitted [arXiv:1611.08552] [Google Scholar]
  83. Muñoz, J. A., & Furlanetto, S. R. 2014, MNRAS, 438, 2483 [NASA ADS] [CrossRef] [Google Scholar]
  84. Muñoz, J. A., & Oh, S. P. 2016, MNRAS, 463, 2085 [NASA ADS] [CrossRef] [Google Scholar]
  85. Nagamine, K., Wolfe, A. M., & Hernquist, L. 2006, ApJ, 647, 60 [NASA ADS] [CrossRef] [Google Scholar]
  86. Neeleman, M., Kanekar, N., Prochaska, J. X., et al. 2017, Science, 355, 1285 [NASA ADS] [CrossRef] [Google Scholar]
  87. Neri, R., Downes, D., Cox, P., & Walter, F. 2014, A&A, 562, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  88. Okamoto, T., Gao, L., & Theuns, T. 2008, MNRAS, 390, 920 [NASA ADS] [CrossRef] [Google Scholar]
  89. Olsen, K. P., Greve, T. R., Narayanan, D., et al. 2015, ApJ, 814, 76 [CrossRef] [Google Scholar]
  90. Olsen, K. P., Greve, T. R., Narayanan, D., et al. 2017, 846, 105 [Google Scholar]
  91. Ota, K., Walter, F., Ohta, K., et al. 2014, ApJ, 792, 34 [NASA ADS] [CrossRef] [Google Scholar]
  92. Oteo, I., Ivison, R. J., Dunne, L., et al. 2016, ApJ, 827, 34 [NASA ADS] [CrossRef] [Google Scholar]
  93. Pallottini, A., Ferrara, A., Bovino, S., et al. 2017a, MNRAS, 471, 4728 [Google Scholar]
  94. Pallottini, A., Ferrara, A., Gallerani, S., et al. 2017b, MNRAS, 465, 2540 [NASA ADS] [CrossRef] [Google Scholar]
  95. Pentericci, L., Carniani, S., Castellano, M., et al. 2016, ApJ, 829, L11 [NASA ADS] [CrossRef] [Google Scholar]
  96. Pineda, J. L., Langer, W. D., Velusamy, T., & Goldsmith, P. F. 2013, A&A, 554, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. Pineda, J. L., Langer, W. D., & Goldsmith, P. F. 2014, A&A, 570, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  98. Planck Collaboration XLVII. 2016, A&A, 596, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  99. Popping, G., van Kampen, E., Decarli, R., et al. 2016, MNRAS, 461, 93 [NASA ADS] [CrossRef] [Google Scholar]
  100. Rawle, T. D., Egami, E., Bussmann, R. S., et al. 2014, ApJ, 783, 59 [NASA ADS] [CrossRef] [Google Scholar]
  101. Riechers, D. A., Bradford, C. M., Clements, D. L., et al. 2013, Nature, 496, 329 [Google Scholar]
  102. Riechers, D. A., Carilli, C. L., Capak, P. L., et al. 2014, ApJ, 796, 84 [NASA ADS] [CrossRef] [Google Scholar]
  103. Robertson, B. E., Ellis, R. S., Furlanetto, S. R., & Dunlop, J. S. 2015, ApJ, 802, L19 [Google Scholar]
  104. Sargsyan, L., Samsonyan, A., Lebouteiller, V., et al. 2014, ApJ, 790, 15 [NASA ADS] [CrossRef] [Google Scholar]
  105. Saunders, W., Rowan-Robinson, M., Lawrence, A., et al. 1990, MNRAS, 242, 318 [NASA ADS] [CrossRef] [Google Scholar]
  106. Smit, R., Bouwens, R. J., Franx, M., et al. 2012, ApJ, 756, 14 [NASA ADS] [CrossRef] [Google Scholar]
  107. Smit, R., Bouwens, R. J., Carniani, S., et al. 2017, ArXiv e-prints [arXiv:1706.04614] [Google Scholar]
  108. Smith, J. D. T., Croxall, K., Draine, B., et al. 2017, ApJ, 834, 5 [NASA ADS] [CrossRef] [Google Scholar]
  109. Song, M., Finkelstein, S. L., Ashby, M. L. N., et al. 2016, ApJ, 825, 5 [NASA ADS] [CrossRef] [Google Scholar]
  110. Spilker, J. S., Marrone, D. P., Aravena, M., et al. 2016, ApJ, 826, 112 [Google Scholar]
  111. Stacey, G. J., Geis, N., Genzel, R., et al. 1991, ApJ, 373, 423 [NASA ADS] [CrossRef] [Google Scholar]
  112. Stacey, G. J., Hailey-Dunsheath, S., Ferkinhoff, C., et al. 2010, ApJ, 724, 957 [NASA ADS] [CrossRef] [Google Scholar]
  113. Strandet, M. L., Weiß, A., De Breuck, C., et al. 2017, ApJ, 842, L15 [NASA ADS] [CrossRef] [Google Scholar]
  114. Swinbank, A. M., Karim, A., Smail, I., et al. 2012, MNRAS, 427, 1066 [Google Scholar]
  115. Tweed, D., Devriendt, J., Blaizot, J., Colombi, S., & Slyz, A. 2009, A&A, 506, 647 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  116. Vallini, L., Gallerani, S., Ferrara, A., & Baek, S. 2013, MNRAS, 433, 1567 [Google Scholar]
  117. Vallini, L., Gallerani, S., Ferrara, A., Pallottini, A., & Yue, B. 2015, ApJ, 813, 36 [NASA ADS] [CrossRef] [Google Scholar]
  118. Wagg, J., Wiklind, T., Carilli, C. L., et al. 2012, ApJ, 752, L30 [NASA ADS] [CrossRef] [Google Scholar]
  119. Walter, F., Decarli, R., Carilli, C., et al. 2012, Nature, 486, 233 [Google Scholar]
  120. Watson, D., Christensen, L., Knudsen, K. K., et al. 2015, Nature, 519, 327 [NASA ADS] [CrossRef] [Google Scholar]
  121. Wiesenfeld, L., & Goldsmith, P. F. 2014, ApJ, 780, 183 [NASA ADS] [CrossRef] [Google Scholar]
  122. Willott, C. J., Carilli, C. L., Wagg, J., & Wang, R. 2015, ApJ, 807, 180 [NASA ADS] [CrossRef] [Google Scholar]
  123. Wolfire, M. G., McKee, C. F., Hollenbach, D., & Tielens, A. G. G. M. 2003, ApJ, 587, 278 [NASA ADS] [CrossRef] [Google Scholar]
  124. Yamaguchi, Y., Kohno, K., Tamura, Y., et al. 2017, ApJ, 845, 108 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: [CII] excitation temperature

The excitation temperature of the [CII] transition is defined by the relative populations of the upper and lower levels, nu and nl, respectively, through the standard equation (A.1)where T∗ is the equivalent temperature (=hν/k), and gu (gl) is the statistical weight of upper (lower) level. The upwards and downwards rate coefficients are related by detailed balance (A.2)where Tkin is the kinetic temperature. Due to the wide range of conditions under which it is the dominant form of carbon, collisional excitation of [CII] by electrons, H, and H2 can all be important. For a single collision partner, the collision rates are equal to the rate coefficients times the density n of that collision partner, thus (A.3)For a region with multiple collision partners, the upwards and downwards rates are the sum of the rates produced by each.

The energy density in the cloud at the frequency of the [CII] transition is given by (A.4)where β is the photon escape probability and Tbg is the temperature of the isotropic CMB radiation field.

For a spherical cloud with a large velocity gradient of the form v ∝ r, the escape probability is given by (A.5)where τ is the peak optical depth of the transition.

Radiative processes include spontaneous emission (rate Aul s-1), stimulated emission (rate BulU), and stimulated absorption (rate BluU). The stimulated rate coefficients are again related by detailed balance through (A.6)\newpageFrom the relationship between the stimulated and spontaneous downwards rates, (A.7)Following Goldsmith et al. (2012), for convenience in dealing with the background, we define (A.8)The rate equation that determines the level populations includes collisional and radiative processes, and is (A.9)The expression for the excitation temperature finally becomes (A.10)The optical depth can be written (A.11)with τ0 being the optical depth which would occur if there were no excitation (i.e. Tex = 0) and (A.12)In this equation, the line profile function at line centre is approximated by δv-1, and N(C+) is the total column density [CII].

Appendix B: Measured L[CII] and SFR for z > 4 star-forming galaxies

We give in Table B.1 a compilation of measured L[CII] and SFR for high-redshift star-forming galaxies.

Table B.1

Compilation of z > 4 star-forming galaxies (i.e. excluding known QSO or AGN) with both [CII] and SFR (or IR luminosity) measurements (so excluding upper limits).

Table B.1

continued.

Appendix C: CMB effect on the L[CII]–SFR relation

We show on Fig. C.1 the L[CII]–SFR relation obtained when ignoring the heating and attenuation of the [CII] line emission by the CMB.

thumbnail Fig. C.1

Same as Fig. 6 but without taking into account CMB effects (both heating and attenuation).

All Tables

Table 1

The ten metallicities considered in our grids of models and the abundances of the five main elements of the ISM we are including in CLOUDY.

In the text
Table B.1

Compilation of z > 4 star-forming galaxies (i.e. excluding known QSO or AGN) with both [CII] and SFR (or IR luminosity) measurements (so excluding upper limits).

In the text
Table B.1

continued.

In the text

All Figures

thumbnail Fig. 1

Observed and predicted stellar mass functions from z = 4 to z = 7. Observations are from Song et al. (2016), Duncan et al. (2014), Grazian et al. (2015), and Caputi et al. (2015). Our SAM predictions are shown in grey.
In the text
thumbnail Fig. 2

Observed and predicted UV luminosity functions at z ~ 4. Observations are from Bouwens et al. (2015), Finkelstein et al. (2015). Our SAM predictions are shown in grey (with and without extinction correction, solid and dotted lines, respectively).
In the text
thumbnail Fig. 3

[CII] excitation (or spin) temperature, Tex, as a function of total gas density n and kinetic temperature on interstellar gas, Tkin (computed for Tbg = 2.726 K). The upper (lower) panel is dedicated to optically thin (thick) medium. The black solid and dot-dash lines correspond to Tex = 16 K and Tex = 21 K, that are the CMB temperature at z = 5.0 and z = 6.7, respectively. At these redshifts, [CII] emission is suppressed for kinetic temperature below these lines, due to the CMB. For optically thin emission, this suppression affects mostly the cold neutral medium (Tkin ~ 50–120 K, n ~ 20–200 cm-3).
In the text
thumbnail Fig. 4

Example [CII] emission (top) and kinetic temperature (bottom) computed by CLOUDY. Two redshifts are shown (blue and red lines). On the top figure one can see the effect of CMB attenuation on the [CII] line luminosity. In this example, the PDR is homogeneous with log nH = 2.4, and exposed to a radiation field with log ISRF = 3.5, in addition to the CMB radiation.

In the text
thumbnail Fig. 5

Physical sizes of effective PDRs in our simulated galaxies at z = 4.7. Orange and blue histograms are associated to equivalent PDR scale-height (heq = hpdr) and equivalent PDR radius (req), respectively.
In the text
thumbnail Fig. 6

[CII] luminosities predicted by CLOUDY (for a CMB temperature at z = 5). Each panel is dedicated to a given metallicity bin and shows the [CII] luminosity per unit area as a function of both ISRF and hydrogen density. In each panel we show the location of G.A.S. galaxies (extracted from the snapshot at z = 4.7). Rectangles represent the 0.15 and 0.85 quantiles of the galaxy distribution. The black point marks the median, with its coordinates given in the top left corners. We also give the number of galaxies Nobj present in the considered metallicity bin.
In the text
thumbnail Fig. 7

L[CII]–SFR relation. Predictions from our model are shown for a set of redshifts from z = 4 to z = 7.6. In each panel the whole sample of G.A.S. galaxies is shown in grey scale. The average relation is plotted with a solid black line. The black dashed line shows the relation given in Eq. (10). Yellow to red coloured points mark the gas metallicity of a randomly selected sample of simulated galaxies (note that the observed tendency of high-metallicity galaxies to fall either above or below the mean trend, that is to make an envelope, is only a trick of the eye caused by the plotting; galaxies with high metallicities (Zg> 8.8) are spread over the whole area, with a higher density of objects at high SFR). Our predictions are compared to a large sample of observational data that are detailed in Table B.1. Amplification corrections on luminosity and SFR, when available, are applied. For dusty star forming galaxies, SFR are converted directly from LIR using the Kennicutt (1998) conversion factor assuming a Chabrier (2003) IMF where SFR (M⊙ yr-1) = 1.0 × 10-10LIR(L⊙). The blue solid line shows the De Looze et al. (2014) relation for the local dwarf galaxy sample.

In the text
thumbnail Fig. 8

Same as Fig. 7 but with the SFR determined from the UV (observed, so attenuated) and IR emission, following Kennicutt (1998) standard conversions between LUV and LIR and SFR. This mimics what is being doing when computing the SFR from the UV and IR emission of galaxies. Taking the “observed” SFR (this figure) rather than the “true” SFR from the model (Fig. 7) decreases the dispersion and removes a large fraction of the outliers.

In the text
thumbnail Fig. 9

L[CII]–SFR relation for galaxies that follow the selection of the simulated galaxy sample by Olsen et al. (2017) at z ~ 6. Our 575 selected galaxies are shown with coloured points (with colour scale reflecting the gas metallicity). Olsen et al. (2017) sample is shown with grey points. The black and the grey solid lines show the mean trend of our whole galaxy sample (Eq. (10)) and of Olsen et al. (2017) sample (their Eq. (6)), respectively. The black dashed line marks the mean trend of our galaxy selection.

In the text
thumbnail Fig. 10

Distribution of the simulated galaxies in [CII] luminosities and metallicities at z = 5.9. Coloured points mark the SFR of a randomly selected sample. Blue squares are from the local dwarf galaxy sample (Cormier et al. 2015).
In the text
thumbnail Fig. 11

Ratio of our predicted [CII] luminosity at z ~ 6 (L[CII]) and that predicted by Vallini et al. (2015) model (). The top panel considers from Eq. (9), and thus N = 1, while the bottom panel considers from Eq. (11), and thus N = 2. N is the power law index of the Kennicutt-Schmidt relation, .

In the text
thumbnail Fig. 12

L[CII ]/LIR vs. LIR for our sample of simulated galaxies (grey shaded areas) at z = 4.7. Local galaxies (orange circles) are the GOALS luminous infrared galaxy sample (Díaz-Santos et al. 2017); a mean correction has been applied to their L[ CII] to mimic the effects of CMB heating and attenuation. High-resdhift galaxies (blue circles) are extracted from Table  B.1. We give on the bottom of the figure the mean densities, ISRF, and metallicities (in log), for bins of 0.5 dex in log  LIR. The [CII] deficit naturally arises in our model. It is well correlated with the intensity of the ISRF.

In the text
thumbnail Fig. 13

[CII] luminosity function predicted by the G.A.S.+CLOUDY model from z = 4.0 to z = 7.6. The blue solid curve shows the prediction that accounts for the attenuation of [CII] emission due to the CMB. The blue dotted line would be the luminosity function ignoring the attenuation. At z ≃ 5, we show the observational constraints from Capak et al. (2015). At z ≃ 6, the black squares indicate the observational results of Yamaguchi et al. (2017). We also add the local [CII] luminosity function published by Hemmati et al. (2017; orange dotted line) and model predictions of Popping et al. (2016; grey solid line).

In the text
thumbnail Fig. 14

UV luminosity function derived from the [CII] luminosity function at z ~ 5. Small and large symbols are observational measurements without and with extinction correction, respectively. The correction of extinction of observed measurements has been done using G.A.S. (using the two grey cuves of Fig. 2). The solid grey line shows our predicted UV luminosity function from G.A.S., which is in good agreement with the observed data points coming from Bouwens et al. (2015) and Finkelstein et al. (2015). The UV luminosity function predicted from our [CII] luminosity function assuming a fixed [CII] to UV luminosity ratio is shown in blue.

In the text
thumbnail Fig. 15

Cumulative [CII] luminosity functions predicted at z ≃ 4 and z ≃ 6 (thick continuous blue and orange lines, respectively, with the thin lines representing ), compared to the current observational constraints (data points with error bars). The dot-dashed lines are derived by relaxing the criterium req<rg which concerns only <0.8% of the galaxies. The dashed lines show the [CII] luminosity functions derived from the carbon mass function of our model.

In the text
thumbnail Fig. C.1

Same as Fig. 6 but without taking into account CMB effects (both heating and attenuation).

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.