THE LOCAL [C ii] 158 μm EMISSION LINE LUMINOSITY FUNCTION

Ideally, constraining the [C ii] line luminosity function down to faint luminosities (∼10^(7–8) L_⊙), would require far-IR spectroscopy of a complete sample of low IR luminosity galaxies. However, such a data set currently does not exist. The alternative best approach is to predict the [C ii] luminosities for RBGS galaxies fainter than LIRGS in the GOALS sample (${L}_{\mathrm{IR}}\lt {10}^{11}\,{L}_{\odot }$). The [C ii] luminosities are calculated using the established correlation between f([C ii])/f(FIR) versus far-IR color, e.g., dust temperature (T_dust) (see Díaz-Santos et al. (2013) and references therein). This relation, as shown on the left panel of Figure 1, is expressed as:

$$\begin{eqnarray}&&\displaystyle \frac{[{\rm{C}}\,{\rm{II}}]}{\mathrm{FIR}}=0.016(\pm 0.001)\times \exp \left(\displaystyle \frac{\tfrac{{S}_{63\mu {\rm{m}}}}{{S}_{158\mu {\rm{m}}}}}{0.60(\pm 0.038)}\right)\end{eqnarray} \tag{ 1 }$$

with a dispersion of 0.0017 dex. We used a modified blackbody function with an emissivity index of β = 1.8 and reference wavelength of 100 μm, that reproduces the observed S_{63 μm}/S_{158 μm} color with the dust temperature shown in the upper axis of Figure 1. The far-IR fluxes covering the 40–500 μm are calculated as: $\mathrm{FIR}=1.26\times {10}^{-14}(2.58{S}_{60\mu {\rm{m}}}+{S}_{100\mu {\rm{m}}})[{\mathrm{Wm}}^{-2}]$ (Helou et al. 2000) where S_ν are in [Jy]. To measure the far-IR color on the right-hand side of Equation (1), we use the correlation between the PACS-based far-IR color ${S}_{63\mu {\rm{m}}}/{S}_{158\mu {\rm{m}}}$ and the more commonly used IRAS-based color ${S}_{60\mu {\rm{m}}}/{S}_{100\mu {\rm{m}}}$.

$$\begin{eqnarray}\mathrm{Log}\left(\displaystyle \frac{{S}_{60\mu {\rm{m}}}}{{S}_{100\mu {\rm{m}}}}\right) & = & -0.161(\pm 0.004)+0.539(\pm 0.018)\\ & & \times \,\mathrm{Log}\left(\displaystyle \frac{{S}_{63\mu {\rm{m}}}}{{S}_{158\mu {\rm{m}}}}\right)\end{eqnarray} \tag{ 2 }$$

As expected, and shown on the right panel of Figure 1, the two colors correlate well, with a dispersion of 0.052 dex.

Zoom In Zoom Out Reset image size

The 1σ uncertainties of [C ii] luminosities for GOALS galaxies are taken from Díaz-Santos et al. (2014). For the non-GOALS RBGS galaxies, however, all the above-mentioned assumptions need to be accounted for. We calculate these errors using a bootstrapping technique. We perturb S_{60 μm}, S_{100 μm} flux densities far-IR fluxes as well as the fitting parameters (from Equations (1) and (2)) by randomly drawing values from normal distributions with their corresponding standard errors as widths of the distributions to measure the [C ii] line luminosity standard error. These 1σ measured uncertainties are largely dominated by the uncertainties in the fitted parameters in Equations (1) and (2) and are an upper limit to the true error as the uncertainty in the fitting parameters already account for uncertainties in fluxes and flux densities. Our method allows us to quantify the uncertainties in the derived [C ii] luminosities for non-GOALs galaxies. Fundamentally, these uncertainties are driven by the intrinsic variations in [C ii]-to-FIR flux ratios. [C ii]-to-FIR ratios are affected by several physical conditions, such as far-IR temperature, neutral gas density, or equivalently surface brightness (e.g., Lutz et al. 2016). [C ii] emission is collisionally excited and can be suppressed in very high-density regions. Another possible physical process leading to [C ii] deficit is the reduction in the photoelectric heating efficiency due to the charging of the dust grains. Therefore, the scattering in the observed [C ii]/FIR versus ${S}_{63\mu {\rm{m}}}/{S}_{158\mu {\rm{m}}}$ relation reflects the physical diversity of ISM in different galaxies.

After excluding the very nearby galaxies (luminosity distances less than 1 Mpc) as well as those with predicted [C ii] line luminosity less than the lowest PACS observed [C ii] line luminosity (10^6.73 L_⊙) from this sample, we are left with 200 GOALS and 395 RBGS non-GOALS galaxies spanning the redshift range of 0.00023–0.076 and [C ii] line luminosities in the range 10^6.73–9.33 L_⊙

The sample we have selected here is complete at flux density ${S}_{60\mu {\rm{m}}}$ > 5.24 Jy. However, to estimate the [C ii] line luminosity function, we need to know how incomplete the sample is at each [C ii] line flux. We therefore calculate a completeness function at each [C ii] line flux (C([C ii], D_L)). This is critical for determining the faint end of the line luminosity function. To estimate the amount of completeness at each [C ii] flux, we make a grid of [C ii] line luminosity and distance. At each cell of this grid we randomly draw a hundred far-IR colors (${S}_{60\mu {\rm{m}}}/{S}_{100\mu {\rm{m}}}$) and use these ratios along with Equations (1) and (2) and their dispersion as well as the definition of far-IR to calculate the corresponding S_{60 μm}. Then the completeness fraction at each cell of the grid is simply the ratio of galaxies with ${S}_{60\mu {\rm{m}}}$ > 5.24 Jy (the selection criteria for the RBGS) to all 100 galaxies in that grid. The completeness function is shown in the left panel of Figure 2. As we are going to use this completeness function in the luminosity function, to speed up the calculation we fit an analytical form to it. The best-fit function as shown on the right panel of Figure 2 is expressed with a sigmoid function as:

$$\begin{eqnarray}&&C([{\rm{C}}\,{\rm{II}}],{D}_{L})\\ &&\quad =\displaystyle \frac{100}{1+\exp (-2.4(\mathrm{Log}([{\rm{C}}\,{\rm{II}}])-4.86{({D}_{L}-4.256)}^{0.2}))}\end{eqnarray} \tag{ 3 }$$

where the unit of [C ii] luminosities is L_⊙ and distances are in Mpc. Completeness for different surveys is often just measured by the turnover in the source counts as a function of brightness. Here, we also show the [C ii] line flux at which the source count starts to drop with a white dashed line on Figure 2 which agrees well with where our completeness function starts to drop.

Zoom In Zoom Out Reset image size

In Figure 2 we also show where our galaxies sit with black crosses (GOALS objects) and plus symbols (the rest of the RBGS) on top of the rainbow-colored completeness values. The galaxies are not originally selected based on their [C ii] luminosity and the majority of the sample sits above 80% completeness. We note that the handful of galaxies with completeness values below 50% are all from the GOALS sample with observed [C ii] luminosities. The completeness values will be used as weights for probabilities of individual sources in the sample. We explain the details of implementing the completeness in the next section.

In addition to the source detection incompleteness, the sky coverage percentage is also accounted for in our calculation. This galaxy sample spans the entire sky except for a thin strip within a galactic latitude of $| b| \lt 5^\circ$. The effective sky coverage is 37,657 square degrees, 91.3% of the full sky (Sanders et al. 2003). We include this multiplicative factor in the luminosity function estimation.

We start with the 1/V_max method to estimate the [C ii] line luminosity function. This method was first discussed by Schmidt (1968) and revised later by Felten (1976). With this method, when we calculate the volume number density of galaxies per Δ(L) for a flux-limited sample, the relevant quantity is the maximum volume (V_max) within which a galaxy could lie and still be detected by the survey. The underlying concept is that a brighter galaxy can be seen further away than an intrinsically fainter one, thus probing a larger volume. This maximum volume is usually constrained by both the maximum and minimum redshift an object could have and still be included in the survey sample.

The minimum redshift (z_min) for all the galaxies in the sample is set by the cut on the luminosity distance (D_L > 1 Mpc) as mentioned in the previous section. The maximum redshift (z_max) a galaxy can have and still be included in the sample is the maximum between the galaxy's actual redshift and that measured by comparing the [C ii] line flux (F_{[C ii]}) to the limiting line flux (Flim_{[C ii]}) of the sample:

$$\begin{eqnarray}&&{z}_{\max }=\max \left(z,0.5\left(\sqrt{1+4z(1+z)\sqrt{\displaystyle \frac{{F}_{[{\rm{C}}{\rm{II}}]}}{{\mathrm{Flim}}_{[{\rm{C}}{\rm{II}}]}}}}-1\right)\right).\end{eqnarray} \tag{ 4 }$$

All cosmological calculations, such as co-moving distance and volume, are simplified for galaxies with small redshifts. Here we adopt the following equations: luminosity distance D_L = $z(1+z)\tfrac{c}{{H}_{0}};$ co-moving volume V_c = $\tfrac{4\pi }{3}{D}_{M}^{3}$, with D_M being proper distance ($\tfrac{{cz}}{{H}_{0}}$) (Hogg 1999). The co-moving maximum volume for each galaxy in the sample is then calculated as:

$$\begin{eqnarray}&&{V}_{\max }=\displaystyle \frac{4\pi }{3}({D}_{{\rm{M}}}^{3}({z}_{\max })-{D}_{{\rm{M}}}^{3}({z}_{\min })).\end{eqnarray} \tag{ 5 }$$

Finally, the luminosity function is the sum of (1/V${}_{\max ,i}$) over all galaxies, divided by the luminosity interval of $\mathrm{Log}({L}_{[{\rm{C}}{\rm{II}}]})=0.2\,{L}_{\odot }$ with which the luminosity function is binned (Table 2).

$$\begin{eqnarray}&&\phi (L)=\displaystyle \frac{1}{{\rm{\Delta }}L}\displaystyle \sum _{i}\displaystyle \frac{1}{{V}_{\max ,i}}.\end{eqnarray} \tag{ 6 }$$

The maximum likelihood estimator is a powerful tool in statistics which estimates the parameters of a model, given the data. In observational cosmology, it was first used by Sandage et al. (1979) (hereafter STY) for deriving luminosity functions of different types of galaxies. As mentioned earlier, this method assumes a functional form for the luminosity function which elliminates the need for binning the data. The most common model used in UV and optical studies is the Schechter function. This is expected as the model was originally derived from a stellar mass function (Press & Schechter 1974; Schechter 1976). In the IR, however, more galaxies have been found in the bright end making a double power-law a better fit to the luminosity function (e.g., Soifer et al. 1987, Patel et al. 2013). For our analysis, we tested both functional forms and found that the double power-law is a much better fit to the shape of the luminosity function, defined as:

$$\begin{eqnarray}&&\phi {(L)={\phi }^{* }{((L/{L}^{* })}^{\alpha }+{(L/{L}^{* })}^{\beta })}^{-1}\end{eqnarray} \tag{ 7 }$$

where ϕ* is the normalization factor, L* is the characteristic luminosity, and α and β are the slopes of the faint and bright end of the luminosity function. To find the optimized parameters of the double-power-law luminosity function given our [C ii] line luminosities and uncertainties, we have to maximize the product (${ \mathcal L }$) of probabilities of finding each galaxy in the sample:

$$\begin{eqnarray}&&{ \mathcal L }=\displaystyle \prod _{i}P({L}_{i}| \alpha ,\beta ,{L}^{* })\end{eqnarray} \tag{ 8 }$$

where the probability of each galaxy in the sample is defined as:

$$\begin{eqnarray}&&P({L}_{i}| \alpha ,\beta ,{L}^{* })={\left[\displaystyle \frac{{\int }_{{L}_{\min }}^{\infty }\phi ({L}^{\prime })\times F({L}^{\prime };{L}_{i},{\sigma }_{i}){{dL}}^{\prime }}{{\int }_{{L}_{\min }}^{\infty }\phi (L){dL}}\right]}^{{w}_{i}}.\end{eqnarray} \tag{ 9 }$$

In measuring the probability for each galaxy the normalization factor ϕ* cancels out and at each luminosity (L_i) this probability can be measured by the three parameters α, β, and L*. Equation (9) is a modified version of what was used in the original STY. Here we account for both incompleteness in the sample and line uncertainties of each galaxy. We account for incompleteness in the sample by including weights (w_i) into Equation (9). These weights are the inverse of the completeness at each luminosity with the simple idea that galaxies with the same luminosity at the same distance should have the same probability of being found. This is similar to incompleteness corrections in spectroscopic samples (see for example Zucca et al. 1994). The uncertainties in [C ii] line luminosity estimates are taken into account by adding the error function ($F({L}^{\prime };{L}_{i},{\sigma }_{i})$) to the probabilities. This is done by assuming a normal distribution for luminosities with 1σ errors as the width of the distribution and summing over all luminosities. While this is negligible for Herschel-detected GOALS sources which have small error bars, it is essential for the rest of the sample. This factor is in essence similar to what was introduced in Chen et al. (2003) for accounting for photometric redshift errors.

In practice, we need to vary parameters L*, α, and β to find the maximum ${ \mathcal L }$ as well as the posterior distribution for each of the parameters. However, this is very computationally expensive for a large enough grid with multiple integration, three free parameters, and over 600 galaxies. Therefore, rather than measuring the probabilities over the whole grid, we perform a random walk Markov Chain Monte Carlo (MCMC) to derive the desired parameters and their uncertainties. To further speed up the calculations, we go to the logarithm space and use the summation of probabilities rather than multiplication to measure ${ \mathcal L }$.

We run our MCMC program with 100,000 steps randomly chosen from the three-dimensional grid (L*, α, and β). The step size in each parameter is not fixed and is drawn from normal distributions. We start by measuring the probabilities as described by Equation (9) with an initial guess for the parameters. These initial guesses do not need to be precise as the first 5%–10% of steps are thrown away (burn-in process) and as long as the order of magnitude is correct the process will converge to the optimized values. At each step, with the jump to the new parameters we calculate the new ${ \mathcal L }$ and compare it to the previous one. In addition, a random acceptance rate scheme is adopted, where the new parameters are accepted and added to the chain if the new likelihood is either larger than the previous one or if their ratio is larger than a random number drawn from a normal distribution. We choose the jump size (width of the normal distributions) to get an acceptance rate of 23%, which is shown theoretically to be the optimal value for an N-dimentional distribution (Roberts et al. 1997).

Figure 3 represents our derived [C ii] line luminosity function from both 1/V_max (blue circles) and maximum likelihood estimator (solid cyan line) methods as well as the posterior distributions of α, β, and L*. We also show on Figure 3 the estimated [C ii] LF based on the GOALS measurements only without any completeness correction (purple squares) using the 1/V_max method. It can be seen that while there is good agreement at the very bright end, the LF starts to be very incomplete and drops at luminosities below $\mathrm{log}10({L}_{[{\rm{C}}{\rm{II}}]})\sim 8.7$. This clarifies the importance of adding the estimated [C ii] fluxes from the rest of the RBGS sample. We note here that the faintest and brightest [C ii] emitters in the sample are from the GOALS galaxies with Herschel PACS observations and by including the rest of the RBGS sample we did not extrapolate to fainter or brighter [C ii] luminosities.

Zoom In Zoom Out Reset image size

The maximum likelihood curve (solid cyan line) is calculated from skewed Gaussian function fits to the posterior distributions with best-fitted parameters shown in Table 1. As can be seen from the figure, there is a good agreement within uncertainties between the maximum likelihood estimator and the 1/V_max methods. In the maximum likelihood estimator method we accounted for sources that might be missing using our derived completeness function, in which we assigned random far-IR color to hundreds of galaxies at each [C ii] luminosity and distance to determine whether they will be detected by our survey. We note that we have drawn the random far-IR colors from a uniform distribution in the color range of our sample due to lack of prior knowledge of the true distribution of sources that we might be missing. The agreement between the two methods of measuring the luminosity function demonstrates the validity of this assumption.

Table 1.  [C ii] Luminosity Function Parameters

Method	α	β	L*(L⊙) × 108	$\phi * ({\mathrm{Mpc}}^{-3}{\mathrm{log}}_{10}{({L}_{[{\rm{C}}{\rm{II}}]})}^{-1})$
Maximum Likelihood	2.36 ± 0.25	0.42 ± 0.09	2.173 ± 0.743	0.003 ± 0.002

Download table as:  ASCIITypeset image

Our completeness simulation is designed to account for potential [C ii] emitters that might be missed due to having S_{60 μm} less than the flux limit of the RBGS sample. However, we have assumed similar far-IR color and dust temperature properties for the galaxies in this simulation to those in the RBGS sample. An important question is whether there exist galaxies with vastly different properties that could change our results. Many recent studies have shown the important role of dwarf galaxies in understanding how galaxies form and evolve in general (e.g., Dekel & Silk 1986; Ferguson & Binggeli 1994, Walter et al. 2007; Tolstoy et al. 2009). More specifically, studies of optical and UV luminosity functions found that dwarf galaxies can contribute significantly to the faint end of the luminosity function (e.g., Liu et al. 2008; Alavi et al. 2016). Madden et al. (2013), using Herschel, provides a rich far-IR and submm photometric and spectroscopic database for local dwarf galaxies covering a large range in metallicity (Dwarf Galaxy Survey—DGS). Focusing on the Herschel PACS spectroscopic data of the DGS, Cormier et al. (2015) found an increasing trend of [C ii] luminosity with increased metallicity among the dwarfs, where they cover a large range of [C ii] luminosites (10^4–9 L_⊙) and metallicities (1/50–1 Z_⊙). The low-metallicity dwarfs which have low [C ii] luminosities have different dust temperature and far-IR colors compared to normal local galaxies (see Rémy-Ruyer et al. 2015). As we are only constraining the [C ii] luminosity function down to $\sim {10}^{7}\,{L}_{\odot }$, the dwarf population in that [C ii] luminosity regime will not have low metallicities ($12+\mathrm{log}({\rm{O}}/{\rm{H}})\gt \sim 8.0$). These galaxies have similar far-IR colors and dust temperatures as the galaxies in our sample and therefore are taken care of in the completeness measurement. However, low-metallicity dwarf galaxies need to be taken into account if the [C ii] luminosity function extends to very faint [C ii] luminosities and they can play an important role in constraining the very faint end of the [C ii] luminosity function.

Understanding the precise evolution of the [C ii] luminosity function would require a larger sample of high-redshift galaxies than already exists. ALMA with full capability is ideal for acquiring such a statistical sample. There however exists limited observations and limited high-redshift [C ii] detections (e.g., Swinbank et al. 2012; Capak et al. 2015; Matsuda et al. 2015). We show on the right panel of Figure 4 where these measurements sit as well as a very rough estimate of the [C ii] luminosity function at high redshifts based on the UV and IR observations.

The cyan circle on the right panel of Figure 4 represents a lower limit from the ALMA detection of [C ii] in two z ∼ 4 galaxies by Swinbank et al. (2012). In that study, they suggested a dramatic increase from z = 0 to z = 4 in the number density at the bright end of the luminosity function in contrast to what we see here. This is solely due to their lower number density estimate at the bright end of [C ii] luminosity function at z = 0. To estimate the z = 0 luminosity function Swinbank et al. (2012) used the Sanders et al. (2003) far-IR luminosity function and the [C ii]/FIR with far-IR luminosity correlation of Brauher et al. (2008). To test the reliability of their method they also used [C ii] luminosities of 227 galaxies compiled by Brauher et al. (2008). As the data come from a complex mix of observations, completeness measurement becomes a big issue.

Using ALMA cycle 1 archival data (in band 7), Matsuda et al. (2015) looked for [C ii] emission in z ∼ 4 galaxies and found no significant emission. They presented upper limits to the z = 4 [C ii] luminosity function which is at least two orders of magnitude larger than the [C ii] luminosity function expected from the UV luminosity function. Capak et al. (2015) observed nine z ∼ 5 normal (∼1–4L_*) star-forming galaxies using ALMA and detected [C ii] in all of them. They reported enhancement in the [C ii] emission relative to the far-IR continuum and therefore a strong evolution in the ISM properties in the very early universe. Blue circles on the right panel of Figure 4 represent a very rough estimate of where these measurements sit compared to the local luminosity function. To do this, we measure the volume for each observation using the area and the redshift width of each ALMA pointing and, as this was a targeted observation of Lyman break galaxies (LBGs), we correct the volume using the number density of these galaxies (Bouwens et al. 2015). Aside from these factors and the low number statistics which make these estimates very sensitive on choice of bins and therefore uncertain, it should be noted that there might exist classes of galaxies that are faint in the UV and optical and therefore not selected as LBGs at high redshifts that are bright in the far-IR and can contribute to the luminosity function of the [C ii] line. Recently, Aravena et al. (2016) identified 14 [C ii] line-emitting candidates using ALMA observations of optical dropout galaxies in the Hubble Ultra-Deep Field in the range 6 < z < 8. Their data points are overplotted on the right panel of Figure 4.

Also shown on the right panel of Figure 4 are estimates of [C ii] luminosity function from UV observations at three different redshifts (z = 0, 2, and 5). As a crude estimate, we start with the UV luminosity function, and we adopt the values and uncertainties from Wyder et al. (2005) for local galaxies, from Alavi et al. (2014) at z = 2 and from Bouwens et al. (2015) at z = 5. To convert this to an IR estimate we use the IRX–β relation (Meurer et al. 1999), which states that the ratio of dust emission in the IR to UV emission (IRX) correlates with the UV spectral slope β. We chose the ratio from the literature (Takeuchi et al. 2012; Alavi et al. 2014; Capak et al. 2015) at each redshift based on the reported average β value and the best corresponding IRX–β curve (i.e., Calzetti-like dust (Calzetti et al. 2000) for z = 0, 2 and SMC-like dust (Gordon et al. 2003) for z = 5). These correspond to log(L_IR/L₁₆₀₀) of 1.0, 0.7, and −0.2 for z = 0, 2, and 5 respectively. The final factor is an assumption for the [C ii]/IR ratio which yields an estimate of the [C ii] luminosity function. Again this average ratio and its range are taken from the literature to be 0.01 [0.001–0.03] at z = 2 (Stacey et al. 2010) and 0.01 [0.003–0.03] at z = 5 (Capak et al. 2015). At z = 0 (orange dashed line) we only show the average curve to compare with our derived local [C ii] luminosity function. Overall the agreement between the UV estimate and the actual derived [C ii] luminosity function is good despite all the assumptions that went into the estimate from the UV. At z = 2 (z = 5), the cyan (purple) dashed line shows the estimated [C ii] luminosity function with the median assumptions and the cyan (purple) shaded region corresponds to the whole possible range based on errors of the UV luminosity function, the IR/UV, and the [C ii]/IR. As can be seen from the width of the shaded regions, these are very uncertain estimates and ALMA observations are needed to constrain the picture. However, assuming the median values (dashed lines) we see a similar evolutionary trend as predicted by Popping et al. (2016) simulations, where the high redshift and local estimates are similar and the rise in the number densities is seen at intermediate redshifts (z ∼ 2).

IR luminosity functions can exclude the uncertainy from the IRX–β but, unfortunately, at high redshifts (z > 3) they only overlap with the brightest part of our [C ii] luminosity function. For example, Gruppioni et al. (2013), using Herschel PACS-selected galaxies, estimated the IR luminosity function out to z = 4, but their highest redshift bins only cover [C ii] luminosities outside the range of our study. However, we use their IR luminosity function at z = 0 and z = 2 and again assume a [C ii]/IR ratio (as in our UV test above) and find perfect agreement for the local measurement and agreement within the errors at z = 2. These estimates of the [C ii] luminosity function are shown on Figure 4 with gray squares on the left panel at z = 0 and at z = 2 with cyan squares on the right panel.

In conclusion, we find tentative evidence which suggests that redshift evolution of [C ii] LF may not be simply linear. The volume density of [C ii] emitters may increase significantly from z = 0 to z = 2, but at z = 5–6, [C ii] LF seems to return back to the similar level as z = 0. This redshift evolution behavior is similar to that of the cosmic star formation rate density which rises from early epochs to its peak value between z ∼ 3 and 1 and drops toward the present time (e.g., Hopkins & Beacom 2006; Khostovan et al. 2015). This evolution of the cosmic star formation rate density is partly explained by the ISM masses and the accretion/consumption of the molecular/total gas (e.g., Walter et al. 2014; Genzel et al. 2015; Scoville et al. 2016).