The IRAM/GISMO two-millimeter survey in the COSMOS field

TL;DR

This paper presents the deepest 2mm continuum survey in the COSMOS field, detecting high-redshift, star-forming galaxies and constraining galaxy evolution models, highlighting the effectiveness of 2mm surveys in selecting massive, high-z galaxies.

Contribution

It provides the widest deep 2mm survey to date, identifies high-redshift galaxies, and compares observations with galaxy evolution models, emphasizing the role of 2mm surveys in high-z galaxy selection.

Findings

Detected four high-significance sources with low false detection probability.

Most detected galaxies are at redshifts above 3, with high star formation rates.

Survey results align with most galaxy evolution models, except those predicting many dusty galaxies at z>7.

Abstract

We present deep continuum observations at a wavelength of 2mm centered on the COSMOS field using the Goddard IRAM Superconducting Millimeter Observer (GISMO) at the IRAM 30m-telescope. These data constitute the widest deep 2mm survey to-date, reaching a uniform $\sigma\sim0.23$mJy beam$^{-1}$ sensitivity over $\sim250$arcmin$^2$ at $\sim\,24''$ resolution. We detect four sources at high significance (S/N$\,\geq\,$4.4) with an expected number of false detection of 0.09 sources, and five sources at $4.4>\,$S/N$\,\geq\,3.7$ with an expected number of false detection of 1.65 sources. Combined with deep GISMO observations in GOODS-N, we constrain the 2mm number counts over one decade in flux density. These measurements agree with most galaxy evolution models tested here, except those with large population of dusty star-forming galaxies at $z>7$. Five GISMO sources have counterparts in (sub-)millimeter catalogs available in COSMOS. Their redshifts suggest that all but one lie above $z\sim3$. These four high-redshift ($z>3$) galaxies have $\tilde{z}=3.9$, SFRs $\sim\,$400 - 1200M$_{\odot}\,$yr$^{-1}$ and $M_{\rm dust}\sim10^{9.5}\,$M$_{\odot}$. They provide a relatively complete selection ($\sim66\%$) of the most luminous ($L_{\rm IR}>10^{12.6}\,$L$_{\odot}$) and highest redshift ($z>3$) galaxies detected within our survey area by AzTEC at 1.1mm. We thus conclude that 2mm surveys favor the selection of massive, vigorously star-forming, high-redshift galaxies. This is corroborated by GISMO-C4, a source with a low false detection probability ($\sim\,6.2\%$), for which the absence of a (sub-)millimeter counterpart supports a high redshift origin ($z\gtrsim3$).

Tables

Table 1. Table 1: GISMO-Detected Source List

ID	R.A.	Decl.	S/N	$S_{2mm}$	$S_{2mm}$			$P_f$b
				raw	deboosteda
					Zavala+19	Béthermin+17	Staguhn+14
	(J2000)	(J2000)		(mJy)	(mJy)	(mJy)	(mJy)	(%)
GISMO-C1	09:59:59.2	+2:34:41.89	5.8	1.29 $\pm$ 0.22	${1.11}_{-0.31}^{+0.20}$	${1.07}_{-0.31}^{+0.20}$	${1.15}_{-0.25}^{+0.25}$	0.0
GISMO-C2	10:00:08.0	+2:26:11.90	5.0	1.09 $\pm$ 0.22	${0.90}_{-0.35}^{+0.18}$	${0.80}_{-0.33}^{+0.19}$	${0.90}_{-0.27}^{+0.27}$	1.4
GISMO-C3	09:59:57.0	+2:27:26.89	4.6	1.02 $\pm$ 0.22	${0.81}_{-0.35}^{+0.18}$	${0.69}_{-0.32}^{+0.20}$	${0.80}_{-0.28}^{+0.28}$	4.0
GISMO-C4	9:59:48.2	+2:26:29.88	4.4	1.01 $\pm$ 0.23	${0.76}_{-0.35}^{+0.21}$	${0.65}_{-0.34}^{+0.20}$	${0.75}_{-0.30}^{+0.30}$	6.2
GISMO-C5	10:00:25.0	+2:23:38.90	3.9	0.86 $\pm$ 0.22	${0.55}_{-0.35}^{+0.22}$	${0.48}_{-0.33}^{+0.13}$	${0.55}_{-0.30}^{+0.30}$	23.3
GISMO-C6	10:00:05.2	+2:25:23.90	3.9	0.84 $\pm$ 0.22	${0.53}_{-0.35}^{+0.21}$	${0.47}_{-0.32}^{+0.13}$	${0.52}_{-0.30}^{+0.30}$	24.1
GISMO-C7	10:00:19.8	+2:32:02.90	3.8	0.84 $\pm$ 0.22	${0.53}_{-0.35}^{+0.21}$	${0.47}_{-0.32}^{+0.13}$	${0.52}_{-0.30}^{+0.30}$	24.1
GISMO-C8	10:00:23.4	+2:29:05.90	3.7	0.82 $\pm$ 0.22	${0.50}_{-0.35}^{+0.19}$	${0.43}_{-0.31}^{+0.13}$	${0.49}_{-0.29}^{+0.29}$	31.7
GISMO-C9	09:59:49.8	+2:32:20.88	3.7	0.82 $\pm$ 0.22	${0.50}_{-0.35}^{+0.19}$	${0.43}_{-0.31}^{+0.13}$	${0.49}_{-0.29}^{+0.29}$	31.7

Table 2. Table 2: Raw and corrected cumulative number counts at 2 mm

Raw number counts			Corrected number counts
$S_{2\mathrm{mm}}^{\mathrm{obs}}$	$N(>$$S)$		$S_{2\mathrm{mm}}^{\mathrm{deboosted}}$	$N(>$$S){}^a$	$N(>$$S){}^b$
(mJy)	(deg-2)		(mJy)	(deg-2)	(deg-2)
This Survey
0.8	96$\pm$19		0.4	353$\pm$147	324$\pm$147
0.9	76$\pm$21		0.5	208$\pm$77	169$\pm$72
1.0	55$\pm$19		0.6	129$\pm$45	96$\pm$43
1.1	38$\pm$17		0.7	87$\pm$32	58$\pm$27
1.2	25$\pm$15		0.8	59$\pm$24	41$\pm$21
1.3	14$\pm$13		0.9	37$\pm$19	26$\pm$16
1.4	18		1.0	23$\pm$15	15$\pm$12
GDF
			0.25	1823$\pm$415
			0.35	997$\pm$291
			0.45	548$\pm$185
			0.55	306$\pm$128
			0.65	158$\pm$105
			0.75	73$\pm$81

Table 3. Table 3: GISMO (sub)mm counterparts

	GISMO-C1	GISMO-C2	GISMO-C3	GISMO-C6	GISMO-C7	GISMO-C4o
	AzTEC8	AzTEC2	AzTEC9	m450.173/850.104	AzTEC5
(Sub-)mm name	AzTEC/C2	AzTEC/C3	AzTEC/C14		AzTEC/C42	…
		450.03/850.00	850.01		450.04/850.03
ALMA namea	C2a	C3a	C14	$\dots$	C42	…
Redshift	3.179i	3.89${}_{-0.67}{}^{}{}_{}{}^{+3.11}$ j	4.58${}_{-0.68}{}^{}{}_{}{}^{+0.25}$ l	1.003m	3.63${}_{-0.56}{}^{}{}_{}{}^{+0.37}$ l	…
$S_{24\mu m}^a$	0.054	0.181	0.054	0.26$\pm$0.02	0.189$\pm$0.013	0.054
$S_{100\mu m}^a$	7.7	5.0	5.0	6.9	5.0	8.5
$S_{160\mu m}^a$	33.5	15.6	10.2	36.1	10.2	17
$S_{250\mu m}^a$	62.7	33.3	14.4	47.6	54.1	14
$S_{350\mu m}^a$	62.0	47.1	22.9	41.4	57.5	19
$S_{450\mu m}^b$	$\dots$	25.3	17.9	8.11$\pm$5.23	25.35$\pm$6.06	25
$S_{500\mu m}^a$	69.8	39.2	30.1	20.7	43.4	25
$S_{850\mu m}^{b,c}$	6.30$\pm$1.43c	11.73$\pm$1.08b,k	11.49$\pm$1.1b	2.12$\pm$0.98b	11.42$\pm$1.38b	5
$S_{870\mu m}^d$	12.3$\pm$3.6	11.47$\pm$2.39k	16.4$\pm$3.3	$\dots$	7.2$\pm$0.2	…
$S_{890\mu m}^e$	21.6$\pm$2.3	9.87$\pm$0.79	7.4$\pm$3.0	$\dots$	9.3$\pm$1.3	…
$S_{1.1\mathrm{mm}}^{f,g}$	5.5$\pm$1.3f	6.6$\pm$1.0f,k	5.8$\pm$1.3f	5.5g	6.5$\pm$1.2f	5.5g
$S_{1.3\mathrm{mm}}^a$	4.07$\pm$0.15	4.50$\pm$0.15	5.01$\pm$0.10	$\dots$	2.39$\pm$0.1	…
$S_{2\mathrm{mm}}^h$	1.09$\pm$0.25	0.68$\pm$0.21	0.75$\pm$0.26	0.50$\pm$0.32	0.50$\pm$0.32	0.70$\pm$0.27
log($L_{\mathrm{IR}}/L_{\odot }$)	$12.6\pm 0.3$	$13.0\pm 0.3$	$13.0\pm 0.2$	$11.6\pm 0.3$	$13.1\pm 0.1$	…
log($M_{\mathrm{dust}}/M_{\odot }$)	$9.9\pm 0.1$	$9.6\pm 0.1$	$9.6\pm 0.1$	$9.2\pm 0.3$	$9.3\pm 0.1$	…
log($M_{\ast }/M_{\odot }$)${}^{a,m}$	${10.97}_{-0.01}^{+0.01}$ ${}^a$	$\dots$	${10.82}_{-0.10}^{+0.01}$ ${}^a$	${10.34}_{-0.08}^{+0.06}$ ${}^m$	${11.46}_{-0.00}^{+0.00}$ ${}^a$	…
$t_{\mathrm{depletion}}$${}^n$ [Gyr]	${2.5}_{-1.3}^{+2.2}$	${0.6}_{-0.3}^{+0.6}$	${0.6}_{-0.2}^{+0.4}$	${5.7}_{-3.8}^{+7.3}$	${0.3}_{-0.1}^{+0.1}$	…

Equations

$$\Delta\alpha=9\farcs 4\times\bigg{(}\frac{{\rm SNR}}{5}\bigg{)}^{-1.4}$$

$$N(>S)=\sum ^{S_{\mathrm{in}}^i>S}\frac{1-P_f(S_{\mathrm{out}}^i/N_{\mathrm{out}}^i)}{A_{\mathrm{eff}}\times C(S_{\mathrm{in}}^i)},$$

$$M_{\mathrm{dust}}/N(M_0<M<M_1),$$

$$N(M_0<M<M_1)=M_{\ast }\frac{{\int }_{M_0}^{M_1}\mathrm{IMF}(M)dM}{{\int }_{M_{\min }}^{M_{\max }}\mathrm{IMF}(M)MdM},$$

Full text

galaxies: evolution – galaxies: formation – galaxies: high redshift – galaxies: starburst – galaxies: photometry – galaxies: luminosity function, mass function

††facilities: IRAM: 30m

1 Introduction

One of the most pressing questions in extragalactic astronomy concerns the production of dust in the very early universe. In particular, the amount of dust that could have reasonably been formed in primordial galaxies within the first few hundred million years after the Big Bang is still heavily debated theoretically and also observationally (e.g. Michalowski, 2015; Mancini et al., 2015; Wang et al., 2017a). Indeed, dust formation in asymptotic giant branch (AGB) star atmospheres is expected to take too long ($\gtrsim$ 500 Myr) to explain the large dust reservoir observed in some z$$\,\gtrsim\,$$5 galaxies; while a significant fraction of the dust formed on short timescales ($<$ 500 Myr) by supernovae (SNe) is expected to be destroyed by the associated shocks. Unfortunately, to-date only a few objects could be studied thoroughly in their cold dust properties for a meaningful estimate of their individual dust masses, including bolometrically bright quasar host galaxies, $\gamma$-ray burst host galaxies, or gravitationally lensed galaxies at redshifts z$$\,\gtrsim\,$$5 (e.g. Michałowski et al., 2010b; Riechers et al., 2013; Hjorth et al., 2013; Michalowski, 2015; Strandet et al., 2017; Zafar et al., 2018; Marrone et al., 2018; Decarli et al., 2018; Venemans et al., 2018). Such observational constraints do not, however, qualify for robust extrapolation to the total dust mass in the early universe, as these galaxies are not necessarily representative of the broader population.

Despite some limitations, work over the past two decades has led to the discovery of considerable populations of dust-rich star-forming galaxies (DSFGs) at high redshifts (though largely found below z$$\,=\,$$4), through deep field survey campaigns at (sub-)millimeter (hereafter (sub)mm) wavelengths (e.g., see reviews of Blain, Smail, Ivison, Kneib, & Frayer, 2002; Casey, Narayanan, & Cooray, 2014). Surveys of such DSFGs benefit from a strongly negative $K$-correction at (sub)mm wavelengths, which renders galaxies of equal total infrared luminosity and cold dust temperature equally likely to be detected in a flux-limited survey, no matter if they reside at $z\sim 1$ or $z\sim 10$. This effect is due to the steep slope of the Rayleigh-Jeans continuum emitted by interstellar dust grains heated by the radiation field of young stars in star-forming galaxies. Unfortunately, distinguishing in (sub)mm surveys the dustiest galaxies at the highest redshifts ($z>4$) from lower-redshift galaxies is exceedingly challenging due to the difficulty in making spectroscopic identifications at those epochs; or is observationnally expensive if one used multiple (sub)mm color criteria such as the Herschel-Red sources (e.g., Cox et al., 2011; Riechers et al., 2013; Ivison et al., 2016). In addition, current (sub)mm surveys being mostly performed at 850$$\,\mum\,<\,$$\lambda_{\rm obs}<1.3\,mm, at the highest redshifts ($z>4$) they do not probe sufficiently long rest-frame wavelengths to yield accurate dust mass estimates (i.e., $\lambda_{\rm rest}\gtrsim 250\,\mu$m; Scoville et al., 2016).

In principle, the 2 mm window provides an ideal setup for deep field surveys aiming at statistical probes of dust-rich objects at these early epochs. It probes sufficiently long rest-frame wavelengths for accurate dust mass estimates, while “filtering out” lower redshift dusty sources at $z\sim 2$, which “contaminate” the 850 $\mu$m and 1.3 mm surveys (Casey et al., 2018, 2018b). Unfortunately, there have not been to-date the observational datasets to test this latter hypothesis. Indeed, until recently the 2 mm window has mainly been explored by large-scale surveys, such as those performed by the South Pole Telescope (SPT; Vieira et al., 2010), with detection limits of few milli-Janskys and thus restricted to the study of gravitationally-lensed galaxies.

Here we present a pioneering survey in this largely unexplored 2 mm window using the Goddard IRAM Superconducting Millimeter Observer (GISMO) as guest instrument at the IRAM 30m-telescope (Staguhn et al., 2014, hereafter S14). In order to complement the deep pencil-beam observations with GISMO towards the northern Great Observatories Origins Survey (GOODS-N; described in S14) field, we targeted a \sim\,$$0.1\,deg2 area towards the central portion of the equatorial 2 deg2 COSMOS field (Scoville et al., 2007). Our GISMO 2 mm survey covers a part of the COSMOS field with the deepest multi-wavelength coverage, including Hubble/WFC3 data from the CANDELS survey (Koekemoer et al., 2011). It overlaps with all (sub)mm surveys undertaken in COSMOS, including deep AzTEC/JCMT (Scott et al., 2008) and AzTEC/ASTE (Aretxaga et al., 2011) 1.1 mm surveys and Scuba-2/JCMT 450 and 850 $\mu$m surveys (Casey et al., 2013; Geach et al., 2017; Wang et al., 2017b). In addition, sub-millimeter galaxies (SMGs) selected from these surveys benefited from numerous interferometric follow-ups (SMA/CARMA/PdBI/ALMA; Younger et al., 2007, 2009; Aravena et al., 2010; Smolčić et al., 2012; Brisbin et al., 2017), offering precise locations and redshift estimates for a large fraction of them (Brisbin et al., 2017; Miettinen et al., 2017a), including some distant (z$$\,>\,4) starbursts, a sub-class of SMGs identified only in recent years (Smolčić et al., 2012; Riechers et al., 2014; Smolčić et al., 2015; Miettinen et al., 2015; Brisbin et al., 2017). This wealth of deep pan-chromatic ancillary data available within our GISMO survey is vital to test the hypothesis that the 2 mm window unveils the dustiest and most distant starbursts of the Universe.

The focus of this paper is two-fold. While we will analyze the individual detections with respect to their redshift and multi-wavelength properties, we will thoroughly discuss statistical implications arising from a novel constraint of the bright-end of the 2 mm source counts. The latter have been previously predicted by several authors using different semi-empirical techniques (Béthermin et al., 2017; Casey et al., 2018; Zavala et al., 2018) and are supposed to be a useful tool to distinguish between different model universes – from the very dust-rich to the very dust-poor (Casey et al., 2018). Providing an observational constraint on the 2 mm source counts is thus paramount to elaborate on the best future strategies to statistically explore the interstellar medium at the earliest cosmic epochs. Throughout this paper we assume a Planck cosmology, adopting $H_{0}=67.8$ (km/s)/Mpc, $\Omega_{M}=0.308$ and $\Omega_{\Lambda}=0.692$ (Planck Collaboration et al., 2016). A Chabrier (2003) initial mass function (IMF) is used for all stellar mass and SFR measurements in this article.

2 The COSMOS-GISMO 2 mm survey

2.1 Observations

The GISMO observations of the COSMOS field were obtained in pooled campaigns over the course of three years (April 2012, April 2013, April 2014, October 2014 and February 2015)222The corresponding IRAM project IDs are 247-11, 227-12, 242-13, 117-14 and 232-15. at the IRAM 30-m telescope as a combined open/guaranteed time program for a total of 113.6 hrs (incl. $\sim\,$37% calibration/instrumental overheads, with a total of $\sim\,$71 hrs on target).

GISMO consists of $8\times 16$ pixels (Staguhn et al., 2008) with super conducting transition edge sensors (TES). The TES are read out by time domain SQUID multiplexers built at the National Institute for Standards (NIST), in Boulder, Colorado (Irwin et al., 2002). The GISMO bandpass has a full width at half maximum (FWHM) of $\sim\,$25 GHz around its peak at 150 GHz (i.e., 2 mm). Pixels are spaced by $13\farcs 75$ and the instantaneous field-of-view is 1\farcm 8$$\,\times\,$$3\farcm 7. More details about the instrument can be found in S14.

The beamsize of GISMO observations at the 30-m telescope is $16\farcs 6$ FWHM. The observations were generally carried out under stable atmospheric conditions and $<\,$4 mm perceptible water vapor, in other words, with a zenith opacity of \tau_{\rm{2mm}}$$\,<\,$$0.11. On average we found 100 pixels to be working during our observations. Focus in $z$-direction was regularly monitored (four times a day) and pointing was frequently checked (once per observing hour) using the nearby bright quasars J1055+018, J0851+202 and J0823+033. Fluxes were calibrated to $<\,$10% accuracy by monitoring Mars, Uranus, and Neptune and employing the atmospheric transmission model of the Caltech Submillimeter Observatory333http://www.submm.caltech.edu/cso/weather/atplot.shtml and the 30-m telescope 225 GHz radiometer readings. We employed an on-the-fly raster scan pattern comprised of 41 $\times$ 10–11s (depending on target elevation) subscans and additional 3s turnaround overhead to cover a total square area of 20\arcmin$$\,\times\,$$20\arcmin per full scan. This observing mode works in total power, hence the signal is not modulated by switching the secondary mirror. Per subscan, additional signal was obtained over an area corresponding to the array extent. The short subscan duration resulted in fast scanning speeds between 110–120″ per second and was chosen to minimize the impact of $1/f$ noise on our large map. We changed the scanning direction by 45 deg about the current azimuth between every two scans in order to assure a homogeneous coverage of the target area and reduce systematic effects in the data. We typically observed for contiguous 7–8 hr blocks and hence profited from earth rotation to further enhance a homogeneous target coverage. Consequently, the resulting multi-season map is circular and reaches a uniform \sigma$$\,\sim\,$$0.23 ($0.3$) mJy/beam sensitivity over an area with $\sim\,$17$\farcm$8 (21$\farcm$4) diameter, i.e. $\sim\,$0.07 (0.1) deg2. At the radius of the full map (i.e., $12\arcmin$) around the central coordinates ($\alpha=10^{\rm h}\,00^{\rm m}\,19\fs 75,\ \delta=+02\arcdeg\,32\arcmin\,04\farcs 40$; J2000) a \sigma$$\,\sim\,$$0.4\,mJy/beam sensitivity is reached.

2.2 Data Reduction

The data have been reduced using version 2.30-4 of the CRUSH software (Kovács, 2008) in deep mode used for the accurate recovery of point sources (for details see S14). In this process the map is spatially filtered above 60″ FWHM to remove spatially variant atmospheric residuals and smoothed with a Gaussian FWHM of 17$\farcs$5 (i.e., $\sim$matched-filtered) to optimize the detection of point sources. The resulting effective image resolution is 24″ FWHM.

To estimate the noise, we produced “jackknifed” maps with CRUSH by randomly multiplying each scan by $+1$ or $-1$, eliminating any stationary noise (including sources and foregrounds) but retaining random noise (including that from the atmosphere). From these jackknifed realizations, we generated the noise map associated with the regular smoothed and filtered map. Fig. 1 shows the histogram of the signal-to-noise ratio (S/N) of one of these jackknifed, filtered and smoothed realization compared to the S/N histogram of the regular filtered and smoothed map. The jackknifed distribution is well fit by a Gaussian with $\sigma=1$, demonstrating that noise in the GISMO map is indeed Gaussian. The histogram for the regular map is also well described by this Gaussian, except at high S/N where resolved sources create an asymmetric excess. This agreement demonstrates that the regular map does not contain a significant stationary noise contribution from unresolved sources (i.e., confusion noise) that would appear as a symmetric widening of the distribution with respect to the jackknifed expectation. This is in clear contrast to the GISMO Deep Field (GDF; S14) in which a significant fraction of the total noise comes from confusion. Indeed, in the GDF, the integration time per beam is long enough for the random noise ($\sigma_{\rm rand.}\propto t^{-1/2}$) and the confusion noise ($\sigma_{\rm confusion}$) to be comparable and thus for a significant fraction of the total noise ($\sigma_{\rm tot}$) to come from confusion ($\sigma_{\rm tot}^{2}=\sigma_{\rm rand.}^{2}+\sigma_{\rm confusion}^{2}$). In our COSMOS survey, $\sigma_{\rm rand.}$ dominates and $\sigma_{\rm confusion}$ can be neglected. Note that the subtraction of sources truncates the S/N histogram for the residual map at our 3.7$$\,\sigma detection threshold (see Sect. 2.3).

2.3 Source Extraction

Source extraction was performed on the smoothed and filtered signal and noise maps produced by CRUSH. The value of these beam-smoothed maps provides the amplitude and uncertainty of a fitted point spread function (PSF) at each position, which in case of point sources corresponds directly to their fluxes and associated uncertainties (see S14). To extract these fluxes, we employed a standard top-down peak-finding approach: (i) we search for the highest peak in the signal-to-noise ratio map; (ii) the flux, noise and position of this source is cataloged; (iii) it is removed from the smoothed and filtered map using a scaled version of our effective PSF (see below) and a new signal-to-noise ratio map is generated; (iv) we repeat this procedure until the highest peak in the signal-to-noise ratio map is below 3$\sigma$, i.e., a level under which sources are deemed to be unreliable. The effective PSF used in this process was generated by convolving the instrumental 16$\farcs$6 FWHM Gaussian beam with a 17$\farcs$5 FWHM Gaussian ($\sim$matched-filtering) and a negative 60″ FWHM Gaussian bowl (resulting from our spatial-filtering above 60″ FWHM; Sect. 2.2). The resulting PSF has a peak of one, zero integral (no DC sensitivity due to filtering; S14) and an 24″ FWHM.

In Sect. 2.3.3, we further refined the S/N threshold above which sources are considered as reliable. We cut our final catalog at 3.7$\sigma$, where the overall false-detection rate is 20%, but particularly highlight sources above 4.4$\sigma$ where the overall false-detection rate is only 2.5%. A total of nine sources were found with S/N$\geq 3.7$, among which four have S/N$\geq 4.4$, within the uniform \sigma$$\,\sim\,$$0.23 mJy/beam sensitivity area of our map (Fig. 2; Tab. 1). In the rest of the paper, we strictly limit our analysis to this homogeneously deep part of the map, which covers $\sim\,$250 arcmin2.

2.3.1 Flux boosting

Low angular resolution maps in the millimeter have the potential to suffer from a set of boosting effects which result in measured fluxes for sources higher than their intrinsic values. Flux boosting can be caused by Eddington bias (in the event that the intrinsic number counts at the observed wavelength are steep) and confusion noise (the high density of sources per beam boosts the measured flux density of any single measured source). To evaluate this flux boosting as a function of the observed signal-to-noise ratio (where the signal is the measured “raw” flux density), we performed extensive Monte-Carlo simulations. We inserted fake sources following a realistic number count distribution into our jackknifed maps and recovered their positions, fluxes and flux uncertainties applying the same source extraction method as that used for our real map.

Here, we adopted two different number count distribution models, that of Zavala et al. (2018) and Béthermin et al. (2017). The Béthermin et al. (2017) SIDES model uses an halo-occupation distribution modeling and infrared spectral energy distributions (SEDs) drawn from a two star-formation modes galaxy evolution model (Béthermin et al., 2012); while the Zavala et al. (2018) model infers the number counts at 2 mm using the first 3 mm number count measurements combined with the infrared luminosity function models discussed in Casey et al. (2018) and Casey et al. (2018b). These models have slightly different number count distributions as a result of different luminosity function assumptions at $z>4$. This give us the ability to explore how those assumptions affect the flux boosting statistic. It is worth noting, however, that both models agree within the uncertainties with the number counts from our measurements (see Sect. 3.1). Injected sources of both models have fluxes as low as $0.05\,$mJy (i.e., well below our detection threshold) and are randomly positioned into our jackknifed maps. We verified that the signal-to-noise histogram of these simulated maps peaks at zero and has a dispersion of one, similarly as our real map (see Fig. 1). Sources recovered in these simulated maps above 3$\sigma$ are associated with the brightest input source within a radius of 18″. This radius corresponds to the GISMO astrometric accuracy as inferred from Eq. 1 using S/N$\,=\,$3 (Sect. 2.3.4; see also Eq. 9 of S14). For each number count model, we generated 10,000 mock maps, recovering about half a million fake sources above 3$\sigma$.

The output-to-input flux ratio (i.e., $S_{\rm out}/S_{\rm in}$) as a function of the observed signal-to-noise ratio (i.e., $S_{\rm out}$/N) inferred from our simulations is shown in Fig. 3. While at high S/N flux boosting converges towards zero and thereby becomes insignificant, it steadily increases at lower S/N, reaching up to $\sim$0.3 dex at 3.7$\,\sigma$. In the S/N range probed by real detections (i.e., 5.8$\,\geq\,$S/N$\,\geq\,$3.7), the model of Zavala et al. (2018) predicts lower flux boosting statistics ($\sim\,$0.05 dex), as its number count distribution has a shallower slope in the corresponding flux density range than that of Béthermin et al. (2017) (see Sect. 3.1). These differences are, however, well within the uncertainties of each model and their effects on the measured number counts are further discussed in Sect. 3.1.

In Tab. 1, we tabulate the “deboosted” fluxes and uncertainties of our real detections, as inferred from both models. Deboosted fluxes correspond to the mode of the distributions at a given S/N, while the upper and lower uncertainties correspond to their 16th and 86th percentile, respectively. We also tabulated in Tab. 1 deboosted fluxes and uncertainties obtained using the analytical approach described in S14 (see their Sect. 4.2) and assuming the number count distribution of Béthermin et al. (2011). These deboosted fluxes are fully consistent within their uncertainties with those inferred here.

2.3.2 Completeness

We also used our simulations to evaluate the source detection completeness of our final catalog. The completeness was defined as the ratio of recovered sources above a given S/N to the total number of input sources at a given flux density. Completeness was evaluated in bins of input flux density and for a detection significance $\geq\,$3.7$\,\sigma$, which was the final cut of our catalog described in the next subsection. The completeness was evaluated for both simulations but only that inferred from the model of Zavala et al. (2018) is shown in Fig. 4. Indeed, because the completeness does not depend on the input number count distribution but simply on the map noise properties, both models yield exactly the same source detection completeness function.

The completeness is close to zero at low input fluxes (i.e., $\sim\,$0.3 mJy) but increases rapidly, reaching 40% and 80% at $\sim\,$0.7 mJy and $\sim\,$1.0 mJy, respectively. The completeness shown in Figure 4 is used in Sect. 3.1 to measure the number counts.

2.3.3 False-detection rate

The number of false detections in our map corresponds to the number of pure noise fluctuations reaching a given significance (i.e., S/N) and thus are mistakenly identified as real sources. In the case of Gaussian noise fluctuations, this number depends only on the considered significance and number of independent Gaussian variables in the map. However, because the number of independent Gaussian variables cannot be trivially calculated for our smoothed and filtered map, we evaluated the expected number of false detection using an empirical approach. We run our source extraction algorithm on 10,000 pure jackknifed maps – which by definition contain only noise – and measured the average number of false detection recovered per map in bins of signal-to-noise ratio (see inset table of false sources identified per map in Fig. 5). At S/N$\,\geq\,$4.4, where our real catalog contains four sources, the expected number of false detection is relatively low and equal to 0.09, corresponding to an overall false-detection rate of $2.5$%. However, at lower significance, the number of false detection rises rapidly, and, for example, no less than 13.74 false detections are expected with S/N$\,\geq\,$3.1. In the real map, we recovered only 27 sources with S/N$\,\geq\,$3.1, leading to an unacceptably high overall false-detection rate of about 51%. We thus decided to cut our final catalog at 3.7$\,\sigma$, where the overall false-detection rate is at an acceptable value of 20%, i.e., 1.74 false detections are expected with S/N$\,\geq\,$3.7 while we detected nine sources at this significance in the real map. Note that the number of false detections predicted by Eq. 2 of S14 is in very good agreement with our empirical estimates.

While measuring number counts, a false detection rate correction needs to be applied to each source according to its detection significance (see Sect. 3.1). Unfortunately, the low number of sources detected in the map does not allow us to infer a false-detection rate at any given significance (i.e., S/N$\,\pm\,\Delta$S/N) without being considerably affected by low number statistics. Therefore, we evaluated these false-detection rates from our simulations, dividing the number of false detections recovered at a given significance (i.e., S/N$\,\pm\,\Delta$S/N) in the 10,000 pure jackknifed maps by the number of sources recovered at this significance in the 10,000 simulated maps of a given model (Fig. 5). While the numbers of false detections remain the same for either models, the numbers of sources recovered in the simulated maps depend on the input number count distribution. Therefore, the two models yield slightly different false-detection rates. At 3.7$\,\sigma$, the false-detection rate measured from the models of Zavala et al. (2018) and Béthermin et al. (2017) are 28% and 38%, respectively, in good agreement with that obtained by dividing the number of false detection with the number of sources recovered in the real map at 4.4$\,>\,$S/N$\,\geq\,$3.7, i.e., 1.65 false detections for five sources in the real map (i.e., 33%). Note that in Sect. 3.1, the false-detection rate of each model is used when appropriate, but in Tab. 1 we tabulated, for each source, the mean false-detection rate from these two models.

2.3.4 Positional uncertainties

Finally, we evaluated the positional uncertainties of our catalog as the difference between the input and recovered positions of injected sources in our simulations. Figure 6 shows the mean, 1$\,\sigma$ and 2$\,\sigma$ positional uncertainties of sources recovered from simulations following the model of Zavala et al. (2018) in bins of observed signal-to-noise ratio. Note that the last signal-to-noise ratio bin is affected by our maximum matching radius of 18″ (see Sect. 2.3.1), and thus likely underestimated. As estimates from the model of Béthermin et al. (2017) are very similar, we did not show them in Fig. 6.

As noted in, e.g., Ivison et al. (2007), S14 and Geach et al. (2017), the positional uncertainties of point-like sources vary with their detection significance (i.e., S/N). S14 defined in their Eq. 9 the maximum allowable separation between a GISMO source and that from other catalogs, taking into account the GISMO positional uncertainties and the 1$\,\sigma_{\rm p}$ catalog position errors. In Fig. 6, we plotted their predictions, setting \sigma_{\rm p}$$\,=\,[math] as there are naturally no intrinsic position errors in our simulated catalog. These predictions should be compared to our 2$\,\sigma$ positional uncertainties, as they correspond to the radius where the counterpart must fall with $\sim$98% confidence. While in broad agreement, these predictions do not perfectly capture the trend with S/N observed in our simulations. We thus fitted these positional uncertainties with a simple power law, parametrized with the signal-to-noise ratio (SNR):

[TABLE]

When matching GISMO-COSMOS sources with source catalogs from the literature, these positional uncertainties should be added in quadrature with the 1$\,\sigma_{\rm p}$ catalog position errors.

2.4 Associations with known COSMOS (sub)mm sources

We searched for the counterparts of our GISMO sources in all relevant (sub)mm catalogs publicly available in the COSMOS field, i.e., the 1.1mm AzTEC/JCMT (Scott et al., 2008) and AzTEC/ASTE (Aretxaga et al., 2011) surveys, and the deep Scuba-2/JCMT 450 $\mu$m and 850 $\mu$m surveys (Casey et al., 2013; Geach et al., 2017). Counterparts were deemed robust if their separations with our sources were lower than the quadratic combination of the GISMO positional uncertainties and the 1-$\sigma$ position errors of the relevant catalog444The source density in the AzTEC and Scuba-2 catalogs are low enough that at this radius the probability that a counterpart is a random association is lower than 5%.. We found counterparts for five GISMO sources: GISMO-C1/AzTEC8, GISMO-C2/AzTEC2, GISMO-C3/AzTEC9, GISMO-C7/AzTEC5 and GISMO-C6/ SCUBA-2 m450.173/850.104. Amongst those, four benefit from intermediate resolution ($\sim 1\arcsec-2\arcsec$) (sub)mm interferometric follow-up of AzTEC sources with the SMA at 890 $\mu$m (Younger et al., 2007, 2009) and ALMA at 1.3mm (Brisbin et al., 2017). This follow-up is key to obtaining robust multi-wavelength and subsequently redshift identification of these galaxies, otherwise impossible with the coarse resolution of single-dish (sub)mm observations (Fig. 7). However, despite this effort, reliable spectroscopic redshift estimates are not yet available for all our sources, mostly due to their faint UV/optical/near-infrared counterparts (Casey et al., 2017). Nevertheless, based on optical/NIR/FIR/mm/radio photometric information, there is an emerging consensus in the literature about their high-redshift nature (Koprowski et al., 2014; Miettinen et al., 2015; Brisbin et al., 2017). In the following, we summarize the current knowledge on the five GISMO galaxies with (sub)mm counterparts and discuss the possible nature of the remaining four sources with no counterparts.

2.4.1 GISMO-C1 / AzTEC8

The brightest galaxy in our GISMO 2 mm survey corresponds to the eighth and second brightest galaxies in the 1.1mm AzTEC/JCMT (AzTEC8 – \theta_{\rm offset}$$\,=\,$$8\farcs 4; Scott et al., 2008) and AzTEC/ASTE (AzTEC-C2 Aretxaga et al., 2011) surveys, respectively. It thus benefited from extensive (sub)mm interferometric follow-up (e.g., Younger et al., 2007, 2009; Brisbin et al., 2017), the latest being performed with ALMA at 1.3mm (Brisbin et al., 2017; Miettinen et al., 2017b). It revealed two millimeter counterparts, with southern and northern component separated by $\sim\,$18″ (Fig. 7). Smolčić et al. (2012) and Koprowski et al. (2014) reported a photometric redshift of $3.17^{+0.29}_{-0.22}$ and $3.15^{+0.05}_{-0.15}$ for the southern component, which in combination with a single CO line detection with CARMA leads to a best available solution of $z=3.179$ (Smolčić et al., 2012; Brisbin et al., 2017, Riechers et al. in prep). The northern component is fainter at 1.3mm than the southern component, and has only a very uncertain photometric redshift estimate ($z_{\rm phot}\lesssim 3$; Brisbin et al., 2017). This component fell below the AzTEC/JCMT detection threshold (Scott et al., 2008) but contributed significantly to the flux density reported by Aretxaga et al. (2011) using the coarse resolution of the AzTEC/ASTE survey (34″ vs. 17″). Our GISMO-C1 detection also exhibits a slight extention towards the northern component. Applying a PSF-fitting analysis to our GISMO map at the position of these two components, we found a very similar 2 mm flux density for the southern component as that reported in Tab. 1 (1.33$\pm$0.22 mJy vs. 1.29$\pm$0.22 mJy), while the north component falls below our detection threshold, i.e., 0.65$\pm$0.23 mJy. Fitting only the southern component did not change its 2 mm flux density, while fitting only the north component yields a $<\,$3.7$\,\sigma$ detection with bad residuals. We thus conclude that the northern component did not significantly contribute to the 2 mm flux density of GISMO-C1. GISMO-C1 and AzTEC8 are most likely colocated near the southern component at a redshift of 3.179.

As described in Younger et al. (2007, 2009), about $\sim\,$2″ east from GISMO-C1/AzTEC8 lies a bright and compact radio galaxy (Fig. 7). This galaxy is not associated with any millimeter emission but to a 24$\,\mu$m-bright low-redshift counterpart. Spitzer and Herschel observations ($\sim\,$5–36″) are likely dominated by emission from this low-redshift interloper, and it is thus impossible to measure reliable 24-to-500$\,\mu$m flux densities for GISMO-C1/AzTEC8. Using Spitzer and Herschel measurements as upper limits, we derived the dust mass and infrared luminosity of GISMO-C1/AzTEC8 via FIR-to-mm SED fitting using the dust model of Draine & Li (2007) (Tab. 3; Fig. 12; see Sect. 3.4 for details). We found $\log(M_{\rm dust}/{\rm M_{\odot}})=9.9\pm 0.1$ and $\log(L_{\rm IR}/{\rm L_{\odot}})=12.6\pm 0.3$, yielding a ${\rm SFR=400^{+390}_{-200}\,{\rm M_{\odot}\,yr^{-1}}}$. Finally, high-resolution ($<0\farcs 05$) observations at $870\mu{\rm m}$ with ALMA have revealed asymmetric structures with multiple clumps in the central kilo-parsec of AzTEC8 (Iono et al., 2016).

2.4.2 GISMO-C2 / AzTEC2/ SCUBA-2 450.03/850.00

GISMO-C2 is associated with the second and third brightest galaxies in the 1.1mm AzTEC/JCMT (AzTEC2 – \theta_{\rm offset}$$\,=\,$$2\farcs 3; Scott et al., 2008; Younger et al., 2007) and AzTEC/ASTE (AzTEC-C3; Aretxaga et al., 2011) surveys, respectively. The source is also coincident with 450.03/850.00 in the SCUBA-2 450$\mu$m and 850$\mu$m maps of Casey et al. (2013). Follow-up with ALMA at 1.3mm revealed two components within the contours of our GISMO detection (Fig. 7; Brisbin et al., 2017). These two components are separated by 3″, the eastern component being \sim\,$$\times\,$$4 brighter at 1.3mm than the western component. Both have counterparts in the 0$\farcs$75 JVLA 3 GHz COSMOS survey (Smolčić et al., 2017). Unfortunately, there are no reliable spectroscopic redshift estimates of either component. A preliminary optical/near-infrared redshift solution of $z=1.123$ has been used in the literature (e.g., Smolčić et al., 2012, 2017; Miettinen et al., 2015, 2017a, 2017c) for the eastern component; however, this solution corresponds to an optical counterpart 1″ offset to the south of the eastern component. That galaxy, visible in the optical behind what appears to be a much lower redshift ($z=0.3$) galaxy whose centroid is another 1″ to the south, could have an associated tentative CO line emission detected with CARMA (see discussion of the ambiguity of this source in Casey et al., 2017, E.F. Jiménez-Andrade priv. comm.). However, further analysis of this source’s obscured SED argues against a $z=1.123$ redshift solution, given the unusually cold dust temperature that such a redshift would imply. Without direct optical/near-infrared counterparts, the only redshift constraint we can place on GISMO-C2 are based on the ALMA-1.3 mm and JVLA-3 GHz radio photometry. For the eastern component, Brisbin et al. (2017) found $z=3.89^{+3.11}_{-0.67}$, while for the western component they reported $z=2.03^{+1.19}_{-0.31}$. Assuming that flux density ratios of these components are the same at 1.3mm and 2mm (i.e., 4.5/1.15), we concluded that GISMO-C2 is dominated by emission from the eastern component, i.e., a galaxy potentially at $z>3$.

It is clear that Spitzer and Herschel observations at the position of GISMO-C2/AzTEC2 are significantly contaminated by emission from the foreground galaxies, potentially both the $z=1.12$ system and the $z\sim 0.3$ galaxy further to the south (Fig. 7). Using Spitzer and Herschel measurements as upper limits, scaling all single-dish (sub)mm flux densities by the flux density ratio observed at 1.3mm, and using $z=3.89$, we found $\log(M_{\rm dust}/{\rm M_{\odot}})=9.6\pm 0.1$ and $\log(L_{\rm IR}/{\rm L_{\odot}})=13.0\pm 0.3$, that yields a ${\rm SFR=1000^{+1000}_{-500}\,{\rm M_{\odot}\,yr^{-1}}}$ (Tab. 3; Fig. 12) for GISMO-C2/AzTEC2.

2.4.3 GISMO-C3 / AzTEC9/ SCUBA-2 850.01

GISMO-C3 is associated with AzTEC9 (\theta_{\rm offset}$$\,=\,$$7\farcs 1; Scott et al., 2008; Younger et al., 2009), AzTEC/C14 (Aretxaga et al., 2011), and SCUBA-2 850.01 (Casey et al., 2013). ALMA follow-up at 1.3mm revealed a single component, well within our GISMO contours (Fig. 7; Brisbin et al., 2017). A tentative spectroscopic redshift of 1.357 was obtained for AzTEC9 from a relatively weak spectrum (Salvato et al., in prep) with DEIMOS at the Keck Telescope. Although this estimate is consistent with the photometric redshift derived by Smolčić et al. (2012) of $1.07^{+0.11}_{-0.10}$, and the photometric redshift of $z=1.45$ from Laigle et al. (2016), it has been suggested that this spectroscopic redshift corresponds to a nearby source unassociated with the submm emission (Koprowski et al., 2014). It has also not been verified in deep near-infrared observations with Keck/MOSFIRE where one might have expected to cleanly detect H$\alpha$ emission (Casey et al., 2017). In contrast, Brisbin et al. (2017) have derived optical/NIR and FIR based photometric redshifts of $4.58^{+0.25}_{-0.68}$ and $4.39\pm 1.39$, respectively, which agree with those derived by Koprowski et al. (2014) of $4.85^{+0.50}_{-0.15}$ and $4.60^{+0.50}_{-0.31}$, respectively. Based on the inconsistency of the low-redshift solution, we suspect the high-redshift solution for GISMO-C3/AzTEC9 is more likely.

GISMO-C3/AzTEC9 does not have any 24$\,\mu$m counterpart. It is associated with emission in the Herschel images but its non-detection in the Scuba-2/JCMT 450 $\mu$m image suggests that a significant fraction of these Herschel flux densities comes from nearby galaxies. Assuming $z=4.58$ and treating Herschel flux densities as upper limits, we found $\log(M_{\rm dust}/{\rm M_{\odot}})=9.6\pm 0.1$ and $\log(L_{\rm IR}/{\rm L_{\odot}})=13.0\pm 0.2$, that yields a ${\rm SFR=1000^{+580}_{-370}\,{\rm M_{\odot}\,yr^{-1}}}$ (Tab. 3; Fig. 12) for GISMO-C3/AzTEC9.

2.4.4 GISMO-C6 / SCUBA-2 m450.173/850.104

GISMO-C6 could be associated with SCUBA-2 m450.173/850.104 (\theta_{\rm offset}$$\,=\,$$7\farcs 5), one of the marginal 3\,<\,$$\sigma$$\,<\,3.6 450$\,\mu$m identified sources with $>\,$3$\sigma$ 850$\,\mu$m counterparts reported by Casey et al. (2013). Though there are only marginal detections at all of these wavelengths, the source is unlikely to be spurious for having been detected in multiple independent datasets. The source is ambiguous because GISMO-C6 is situated $\sim\,$7$\farcs$5 away from the SCUBA-2 position, though this is within the positional uncertainties of both sources. Casey et al. (2013) associated SCUBA-2 m450.173/850.104 to an optical counterpart with 24$\,\mu$m counterpart and $z_{\rm phot}=1.01$ (Fig. 7). Near-infrared spectroscopic follow-up of this source with Keck MOSFIRE in Casey et al. (2017) yields a spectroscopic redshift of 1.003 and subsequently a infrared luminosity of $\log(L_{\rm IR}/{\rm L_{\odot}})=11.81\pm 0.25$. However, deep radio 3 GHz imaging of the COSMOS field (Smolčić et al., 2017), did not show any radio counterpart to this source, leading to a far-infrared-to-radio luminosity ratio, $q$, of $2.8$ at odds with the expected $2.47\pm 0.26$ value at this redshift (Magnelli et al., 2015; Delhaize et al., 2017). The 3 GHz COSMOS image shows, however, two other possible counterparts (Fig. 7) with a photometric redshift of $0.78$ and $1.00$, $\sim\,$6″ northwest and $\sim\,$6″ northeast from the source identified by Casey et al. (2017), respectively. Spitzer and Herschel observations exhibit bright but confused emission most likely associated with these low-redshift galaxies. Similarly, the (sub)mm emission of GISMO-C6/SCUBA-2 m450.173/850.104 could well originate from the combined emission of these low-redshift galaxies, which might be part of a single system. However, it could also be emitted by a yet unknown optically-faint high-redshift galaxy. Deep interferometric (sub)mm observations are needed to further investigate the nature of GISMO-C6 / SCUBA-2 m450.173/850.104.

Assuming a low redshift origin of the (sub)mm emission, i.e., $z=1.003$ as in Casey et al. (2017), and using Herschel measurements as upper limits, we found $\log(M_{\rm dust}/{\rm M_{\odot}})=9.2\pm 0.3$ and $\log(L_{\rm IR}/{\rm L_{\odot}})=11.6\pm 0.3$, that yields a ${\rm SFR=40^{+40}_{-20}\,{\rm M_{\odot}\,yr^{-1}}}$ (Tab. 3; Fig. 12). These estimates are fully consistent with that inferred in Casey et al. (2017).

2.4.5 GISMO-C7 / AzTEC5/ SCUBA-2 450.04/850.03

GISMO-C7 is associated with AzTEC5 (\theta_{\rm offset}$$\,=\,$$1\farcs 6; Scott et al., 2008), AzTEC/C42 (Aretxaga et al., 2011), and SCUBA-2 450.04/850.03 (Casey et al., 2013). The source has interferometric follow-up from both the SMA at 890$\,\mu$m (Younger et al., 2007, 2009) and from ALMA at 1.3mm (Fig. 7; Brisbin et al., 2017) revealing a single counterpart. Although there still is no spectroscopic redshift reported in the literature for AzTEC5, Casey et al. (2013) quote a optical/near-infrared photometric redshift, derived in Ilbert et al. (2013), of $z_{\rm phot}=3.82^{+0.44}_{-0.69}$, which is consistent with other far-infrared/radio determinations of the photometric redshift (Smolčić et al., 2012; Koprowski et al., 2014; Brisbin et al., 2017). Recently, thanks to high-resolution observations with HST towards AzTEC5, Gómez-Guijarro et al. (2018) reported that it is has three primary rest-frame UV components, with stellar masses of $\log(M_{\star}/M_{\odot})=9.92^{+0.10}_{-0.10},\,9.78^{+0.08}_{-0.10}$ and $9.59^{+0.08}_{-0.06}$ and 3D-HST-based redshifts of $z=3.63^{+0.14}_{-0.15}$, $z=4.02^{+0.08}_{-0.08}$ and $z=3.66^{+0.40}_{-0.43}$ (Brammer et al., 2012; Skelton et al., 2014; Momcheva et al., 2016). Since these redshifts are consistent within the uncertainties, it is likely that these three components belong to the same system (Gómez-Guijarro et al., 2018), favoring a scenario in which GISMO-C7/AzTEC5 is a merger-driven star-forming galaxy at $z\sim 3.6$.

Assuming $z=3.63$ and treating Herschel measurements as upper limits, we found $\log(M_{\rm dust}/{\rm M_{\odot}})=9.3\pm 0.1$ and $\log(L_{\rm IR}/{\rm L_{\odot}})=13.1\pm 0.1$, that yields a ${\rm SFR=1250^{+320}_{-260}\,{\rm M_{\odot}\,yr^{-1}}}$ (Tab. 3; Fig. 12) for GISMO-C7/AzTEC5.

2.4.6 GISMO sources with no (sub)mm counterparts

Four out of nine GISMO \geq\,$$3.7$$\sigma detections are not associated with any known SMGs, even though they are within deep and available (sub)mm coverage of the COSMOS field (Scott et al., 2008; Aretxaga et al., 2011; Casey et al., 2013; Geach et al., 2017). Such (sub)mm dropouts could be false 2 mm detections. However, the probability of all four being false detections is relatively low, as only 9% of our pure jackknifed maps have $\geq$4 false detections with S/N$\,\geq 3.7$. These (sub)mm dropouts could instead be unidentified high-redshift galaxies (e.g., S14).

GISMO-C4 is the fourth brightest source in our survey. It is detected with high significance and is thus very unlikely to be a false detection, i.e., P_{\rm f}$$\,=\,$$6.2% (Tab. 1). It is within deep Scuba-2/JCMT and AzTEC/ASTE coverages (i.e., $\sim 0.8$ and $\sim 1.26\,$mJy/beam rms at 850$\,\mu$m and 1.1mm, respectively; Casey et al., 2013; Aretxaga et al., 2011). Taking the 80% completeness detection limits of these surveys555these 80% completeness limits are read off Fig. 6 of Casey et al. (2013) – taking the mean value from their 3$\,\sigma$ and 4$\,\sigma$ curves as the final 850 $\mu$m catalog is cut at 3.6$\,\sigma$ – and Fig. 5 of Aretxaga et al. (2011) (i.e., 5 and 5.5 mJy at 850$\,\mu$m and 1.1 mm, respectively), yields 850$\,\mu$m-to-2 mm and 1.1 mm-to-2 mm flux densities ratios $<\,$7.1 and $<\,$7.8, respectively. Such low flux densities ratios indicate that the 850$\,\mu$m and 1.1 mm broadbands do not probe the Rayleigh-Jeans emission of GISMO-C4, because otherwise the 850$\,\mu$m-to-2 mm and 1.1 mm-to-2 mm flux densities ratios would be in the range 20–30 and 8 – 11, respectively. Instead, the 850$\,\mu$m and 1.1 mm broadbands likely probe rest-frame wavelengths shortward of 300$\,\mu$m, supporting a high-redshift solution for this galaxy, i.e., $z>3-4$. This high-redshift candidate will require dedicated follow-up efforts with, e.g., ALMA or NOEMA.

GISMO-C5, GISMO-C8 and GISMO-C9 are all detected at 4.0$\,\geq\,$S/N$\,\geq\,$3.7 at 2 mm, a detection significance range in which we expect 1.22 false detections (Sect. 2.3.3). Given that GISMO-C6 and GISMO-C7 (the other two sources detected in this significance range) have known (sub)mm counterparts, it is likely that at least one of these three (sub)mm dropouts is a false detection. GISMO-C5, GISMO-C8 and GISMO-C9 are within deep AzTEC/ASTE coverage (Aretxaga et al., 2011). However, only GISMO-C5 and GISMO-C8 sit in the deep central 850$\mu$m map of Casey et al. (2013), while GISMO-C9 has shallower 850$\mu$m upper limit from Geach et al. (2017). The 80% completeness detection limits from these surveys666the 80% completeness from Geach et al. (2017) is read off their Fig. 8, i.e., 6.2 mJy, lead to 850$\,\mu$m-to-2 mm and 1.1 mm-to-2 mm flux densities ratios in the range 9.8–13.3 and 10.7–11.8, respectively. Such low flux densities (yet not as low as observed in GISMO-C4) suggests that these galaxies reside at high redshift (i.e., $z\gtrsim 3$). Further follow-up efforts are needed, keeping in mind that at least one of these sources is likely a false detection.

3 Results

3.1 2 mm source counts

In the absence of robust redshift determinations for all our sources, number counts are the most powerful tool to constrain models of galaxy evolution from our catalog. “Corrected” cumulative number counts, i.e., N(>$$\,S), are given by the number of galaxies with a flux density higher than $S$:

[TABLE]

where $S_{\rm in}^{i}$ and $S^{i}_{\rm out}$ are respectively, the deboosted and observed flux densities of the $i$th source, while $N^{i}_{\rm out}$ is its observed flux densities error (Tab. 1); $P_{f}(S^{i}_{\rm out}/N^{i}_{\rm out})$ is its probability to be a false detection (Tab. 1); $C(S^{i}_{\rm in})$ is the completeness of our source extraction at this deboosted flux density (Fig. 4); and, finally, $A_{\rm eff}$ is the area for the source extraction (i.e., 250 arcmin2 where the rms is better than 0.23 mJy/beam). Deboosted flux densities and $P_{f}(S^{i}_{\rm out}/N^{i}_{\rm out})$ being model-dependent, “corrected” cumulative number counts must be evaluated for each model. To avoid these model-dependencies, we will also report here “raw” number counts, by setting $P_{f}(S^{i}_{\rm out}/N^{i}_{\rm out})=0.0$, $C(S^{i}_{\rm in})=1$ and $S_{\rm in}^{i}=S^{i}_{\rm out}$. “Raw” number counts are not easily comparable with past and future literature measurements as these latter naturally suffer from different observational biases (i.e., flux boosting, contamination and completeness). However, “raw” number counts can be compared to model predictions using the mock maps produced in Sect. 2.3 and which reproduce the observational biases affecting our survey.

To properly account for the large, asymmetric and non-Gaussian deboosted flux uncertainties of our sources, we measured $N(>S)$ using a Monte-Carlo approach. We created 1,000 realizations of our catalog, drawing the deboosted flux densities of each source following the flux boosting distribution measured at their respective S/N and shown in Fig. 3. The “corrected” cumulative number counts and associated uncertainties are then given by the average and dispersion of the $N(>S)$ distribution measured over these 1,000 realizations. To validate this approach, we applied this methodology to our simulations. In the upper-left and lower-right panels of Fig. 8, we show the mean and dispersion of $N(>S)$ as measured by applying this methodology to 200 map realizations. On average, this methodology retrieved perfectly the input number count distribution and for 68% of our map realizations it provides measurements within $\sim\,$0.15 dex. However, these tests correspond to an ideal case in which the corrections (i.e., flux boosting statistics and $P_{f}(S^{i}_{\rm out}/N^{i}_{\rm out}$)) were measured on the same model. In reality, we obviously do not know the intrinsic sky model and are thus limited by the model-dependencies of our corrections. To test this effect, we applied corrections measured on one model to the other model. Results are shown in the upper-right and lower-left panels of Fig. 8. In these more realistic cases, we naturally do not retrieve perfectly the input number count distributions. The corrections from Zavala et al. (2018), with their lower flux boosting statistics, slightly over-estimate the number count distribution when applied to catalogs extracted from Béthermin et al. (2017) simulated maps. On the other hand, corrections from Béthermin et al. (2017) slightly under-estimate the number counts when applied to the catalogs extracted from Zavala et al. (2018) simulated maps. These over/under-estimations are, however, comparable to the uncertainties, and certainly, do not lead to number count distributions equal to that from the correction model. Note that the number count distributions followed by these models bracket that inferred from our real catalog, suggesting that they provide a realistic representation of the range of possible corrections.

Finally, using the same simulations, we tested the methodology advocated in S14, i.e., plotting the number of sources at each deboosted flux density, divided by the effective area for the detection of sources. Results of this test for one of our simulations are shown by grey circles in Fig. 8. Irrespective of the corrections used in this methodology, it systematically overestimates the number counts at low flux densities, where the uncertainties on the deboosted fluxes are large. This is understandable when considering the case of two sources with the same deboosted flux densities but large uncertainties: at their common flux density, the number of sources is equal to two but the probability that both have fluxes greater than this value is equal to 0.5, leading to an $\times$2 over-estimation of the cumulative number counts. Note that this test does not imply that the deboosted flux densities quoted in S14 are incorrect; simply that the number counts inferred in S14 from these deboosted flux densities are overestimated at faint flux densities.

Having validated our Monte-Carlo methodology on simulations, we now measure the “raw” and “corrected” GISMO COSMOS 2 mm cumulative number counts using this approach. For our “raw” measurements, this simply implies to create 1,000 realizations of our catalog, drawing the observed flux density of each source following Gaussian distributions characterized by their observed flux uncertainties.

In the left panel in Fig. 9 we show our “raw” 2 mm cumulative number counts, while on the right panel we show our “corrected” cumulative number counts. Table 2 provides the same data in tabular form. In the left panel of Fig. 9, we compare these “raw” measurements with the “raw” predictions from the models of Zavala et al. (2018) and Béthermin et al. (2017). These predictions correspond to the 16th and 86th percentiles of the “raw” cumulative number counts distribution measured using 1,000 mock source catalogs retrieved from 1,000 simulated maps generated in Sect. 2.3. On the right-hand panel, we compare our “corrected” 2 mm cumulative number counts with the intrinsic model predictions. Finally, on the right-hand panel we also compare our “corrected” measurements with those from S14, applying our Monte-Carlo methodology to their Béthermin et al. (2011)-based deboosted flux densities, assuming that their associated uncertainties follow a Gaussian distribution.

Our “raw” cumulative number counts are nicely bracketed by the “raw” predictions from the models of Zavala et al. (2018) and Béthermin et al. (2017). This suggests that these models provide us with reasonable range of possible corrections, and thus robust “corrected” cumulative number count measurements. This also suggests that models predicting significantly more or significantly fewer 2 mm sources in the $\sim$mJy regime are inconsistent with our observations.

Our two different sets of corrections yield very consistent “corrected” cumulative number count measurements, well within their uncertainties. These measurements agree also within their uncertainties with those from S14 in the flux density range where they overlap. This comparison clearly illustrates the advantage of the ”wedding cake” observing strategy followed by the GISMO team. On the one hand, the COSMOS map with its large sky coverage provides critical constraints on the 2 mm number counts at high flux densities, inaccessible to the pencil-beam survey of the GDF. On the other hand, the GDF provides critical constraints at low flux densities, well below the detection threshold of our COSMOS map. Combining these two surveys, we obtain robust measurements of the 2 mm number counts over almost one decade in flux density.

Both the models of Zavala et al. (2018) and Béthermin et al. (2017) are consistent within the uncertainties with our “corrected” measurements. We also include two additional models drawn from Casey et al. (2018), which are identical to the Zavala et al. (2018) model but assume a different evolution in $\Phi_{\star}$ of the obscured luminosity function beyond $z>2$: Model A represents a “dust-poor” early Universe, while Model B represents a “dust-rich” early Universe. Among all models plotted in Fig. 9, the Casey et al. (2018) Model A provides the best description of our estimates777Performing our Monte-Carlo simulations using the model A of Casey et al. (2018), we ended up with deboosted flux densities and “corrected” number count measurements intermediate to those inferred from the models of Zavala et al. (2018) and Béthermin et al. (2017), leaving our results unchanged. Although with our number counts we cannot fully discriminate between the model of Zavala et al. (2018), Béthermin et al. (2017) and Model A of Casey et al. (2018), our data are inconsistent with the most extreme extrapolation of the Casey et al. (2018) Model B, the “dust-rich” early Universe. It is important to point out, however, that such “dust-rich” models are very dependent on the model cutoff redshift (discussed more extensively in Zavala et al., 2018). This cutoff redshift represents a model instantaneous redshift above which no more DSFGs can be found, and letting it vary between $6<z_{\rm cutoff}<12$ changes the expected 2 mm number counts of a “dust-rich” Universe substantially. By tuning down this parameter, starting from its original value of 12 in Casey et al. (2018) and until we find a broad agreement with our number count measurements888The results of this fine tuning is, however, not shown in Fig. 9., we can rule out Model B with $z_{\rm cutoff}\geq 7$.

3.2 The observed (sub)mm-to-2 mm colors

The observed 850$\,\mu$m-to-2 mm and 1.1 mm-to-2 mm flux density ratios provide constraints on the nature and dust properties of galaxies. On the one hand, for galaxies at $z\lesssim 2$, these broadbands probe their Rayleigh-Jeans dust emission, providing a measure of their dust emissivity spectral index, $\beta$. While $\beta$ varies on Galactic scales, extragalactic measurements converge to 1.5\,\leq\,$$\beta$$\,\leq\,2.0 (e.g., Dunne & Eales, 2001; Magnelli et al., 2012), corresponding to 850$\,\mu$m-to-2 mm and 1.1 mm-to-2 mm flux density ratios in the range 8–11 and 20–30, respectively. On the other hand, for galaxies at $z\gtrsim 2$, the 850$\,\mu$m and 1.1 mm broadbands probe closer to the peak of the dust emission, yielding much lower 850$\,\mu$m-to-2 mm and 1.1 mm-to-2 mm flux density ratios, which in absence of robust redshift measurements can be used to support their high-redshift nature.

The observed 850$\,\mu$m-to-2 mm and 1.1 mm-to-2 mm flux density ratios of our nine GISMO galaxies as a function of their 1.1 mm flux densities and, when available, redshifts (see Tab. 3) are displayed in Fig. 10. For sources with no counterparts in the AzTEC 1.1mm or Scuba-2 850$\,\mu$m catalogs, we used as upper limits the 80% completeness limits of these respective surveys (see Sect. 2.4).

The four GISMO galaxies with millimeter counterparts (i.e., GISMO-C1/-C2/-C3/-C7) have a flux density ratio distribution consistent with that observed in the GDF (S14) and in the gravitationally-lensed galaxy sample from the SPT (i.e., South Pole Telescope; Strandet et al., 2016). This agreement demonstrates the quality of the GISMO calibration but also suggests that our galaxies have intrinsic properties similar to galaxies in those samples, i.e, high-redshift highly star-forming galaxies. In addition, our galaxies as well as those from the GDF and SPT, have at their redshifts higher 850$\,\mu$m-to-2 mm and 1.1 mm-to-2 mm flux density ratios than predicted from local SED templates. This implies that high-redshift star-forming galaxies have in average hotter dust temperatures than in the local Universe (see also, e.g., Magdis et al., 2012; Magnelli et al., 2014; Béthermin et al., 2015; Faisst et al., 2017).

GISMO-C6 has a 1.1 mm-to-2 mm upper limit consistent with the rest of the distribution but its very low 850$\,\mu$m-to-2 mm flux density ratio is at odds with the rest of the distribution if this galaxy is at $z=1.003$. This suggests that the redshift association of this source is likely incorrect, as already discussed in Sect. 2.4.

All GISMO galaxies without (sub)mm counterparts have low 850$\,\mu$m-to-2 mm upper limits, consistent with a high-redshift nature. Among those, GISMO-C4 exhibits the lowest 850$\,\mu$m-to-2 mm upper limits. Combined with its low false detection probability ($\sim\,$6%), this makes GISMO-C4 a robust high-redshift ($z>4$) candidate.

Finally, as a sanity check, we evaluated lower limits on the 1.1 mm-to-2 mm flux density ratios of all AzTEC-1.1 mm sources within our deep coverage (\sigma$$\,\sim\,0.23 mJy/beam) but undetected by GISMO at S/N$\,\geq\,$3.7. Here, we used as 2 mm upper limits, $5\,\times$ the value of the GISMO noise at the position of the AzTEC-1.1 mm sources. None of these lower limits are at odds with the rest of the distribution. This implies that our GISMO survey did not “miss” any plausible 2 mm emitters.

3.3 Redshift Distribution

Recent observations as well as phenomenological galaxy evolution models suggested that flux density and wavelength selection are crucial determining factors in the redshift distribution of (sub)mm-selected galaxy samples (e.g, Younger et al., 2007, 2009; Béthermin et al., 2015; Strandet et al., 2016; Brisbin et al., 2017; Casey et al., 2018). In particular, these studies found that selecting bright galaxies at long wavelengths ($\lambda_{\rm obs}>1.1$mm) provides the most favorable criterion for picking out high-redshift star-forming galaxies. With our GISMO-2 mm wide survey, we can further explore these findings.

The three brightest galaxies (i.e., S/N$\,\geq\,$4.6) in our catalog have all a (tentative) redshift measurement: GISMO-C1 has a spectroscopic redshift of 3.179; GISMO-C2 a (sub)mm-to-radio-based photometric redshift of 3.89${}^{+3.11}_{-0.67}$; and GISMO-C3 an optical/NIR-based photometric redshift of $4.58^{+0.25}_{-0.68}$ (Brisbin et al., 2017, see Sect. 2.4 and Tab. 3). In Fig. 11, we show the cumulative redshift distribution of this redshift-complete 2 mm-bright (i.e., S/N$\,\geq\,$4.6) galaxy sample. To account for the large redshift uncertainties for GISMO-C2 and GISMO-C3, we used 1,000 realizations of our catalog, each time drawing the redshift of our sources randomly and uniformly within their redshift uncertainties as given in Tab. 3 (see also Sect. 2.4). In Fig. 11, we plotted the 16th, median and 86th percentiles of these 1,000 cumulative redshift distributions. The median redshift of our 2 mm-bright galaxy sample is $\tilde{z}=4.1$, significantly higher than that of the COSMOS AzTEC/ASTE 1.1 mm sample analyzed by Brisbin et al. (2017), i.e., $\tilde{z}=2.45$. However, restricting their sample to the five galaxies with $S_{\rm 1.3\,mm}\geq 4.07\,$mJy, (the faintest 1.3 mm flux density found within our 2 mm-bright galaxy sample; Tab. 3), their median redshift increases to $\tilde{z}=4.3$, following their conclusion that brighter millimeter sources are preferentially found at higher redshifts. While based on very small numbers, this suggests that bright sources at 1.1 mm and 2 mm surveys yield very similar redshift distribution. However, the possible advantage of 2 mm surveys in picking out higher redshift star-forming galaxies than that at 1.1 mm might only be revealed by probing even larger co-moving volumes. In addition, one has to bear in mind that GISMO-C4 – our next brightest detection – has no AzTEC-1.1 mm counterpart and could thus potentially lie at very high redshifts.

In Fig. 11, we also compare our findings to predictions from the model A of Casey et al. (2018) that best describes our number count measurements (see Sect. 3.1). To this end, we used mock maps generated from this model, similarly as in Sect. 2.3. In Fig. 11, we plotted the 16th and 86th percentiles of the cumulative redshift distributions of sources retrieved with S/N$\,\geq\,$4.6 in 1,000 mock maps. The redshift distribution as well as the median redshift predicted by this model (i.e., $\tilde{z}_{\rm model}=3.8$) is consistent with our observations.

Unfortunately we could not explore the effect of different flux density cuts on the median redshift of our 2 mm-selected samples. Indeed, cutting our catalog at lower S/N would include galaxies for which no redshift informations are yet available (i.e., GISMO-C4, GISMO-C5, GISMO-C8 and GISMO-C9). Such analysis is postponed until further follow-up observations on these galaxies are made.

3.4 FIR-to-mm Spectral Energy Distribution

We derived the infrared luminosities, SFRs and dust masses ($M_{\rm dust}$) of all GISMO galaxies with tentative redshift measurements (i.e., GISMO-C1, GISMO-C2, GISMO-C3, GISMO-C6 and GISMO-C7) via mid-infrared-to-(sub)mm SED fitting using the Draine & Li (2007) dust model999The low number of models compatible with GISMO-C7 is due to its particular photometry. It has a high-significance MIPS-24$\,\mu$m detection and very constraining PACS-100$\,\mu$m and -160$\,\mu$m upper limits with respect to its (sub)mm flux densities. (Fig. 12; Tab. 3). The mid-infrared-to-(sub)mm photometry and redshift used in these fits are summarized in Tab. 3 and were discussed in Sect. 2.4.

As suggested by the observed (sub)mm-to-2 mm colors (Sect. 3.2), the 2 mm flux densities measured by GISMO are consistent with the overall SEDs of these galaxies while putting additional constraints on their Rayleigh-Jeans dust emission. Only GISMO-C6 exhibits an unusual dust SED, which peaks at very long wavelength (\lambda_{\rm rest}^{\rm peak}$$\,\sim\,150$\,\mu$m), corresponding to a low luminosity-weighted dust temperature of $\sim\,$20 K. This low dust temperature could be explained by the well-known $L_{\rm IR}-T_{\rm dust}$ selection bias affecting (sub)mm surveys (e.g., Magnelli et al., 2012). However, this could also suggest that the redshift association for this source is incorrect, as already discussed in Sect. 2.4.

All our galaxies except GISMO-C6, are very luminous (L_{\rm IR}$$\,>\,$$10^{12.6}\,L⊙), corresponding to SFRs in the range 400 – 1200 M${}_{\odot}\,$yr*-1*; assuming a Chabrier IMF and the relation SFR [M${}_{\odot}\,$yr*-1*]$\,=\,10^{-10}\times L_{\rm IR}\,$[L⊙] (Kennicutt, 1998). Even with their relatively large stellar masses ($10^{10.8-11.5}\,$M*⊙*; Tab. 3), these high SFRs locate these galaxies on the upper-part or above the z$$\,\sim\,$$4 main-sequence of star-forming galaxies (e.g., Schreiber et al., 2015). These are thus “starburst” galaxies with respect to the bulk of the star-forming galaxy population at these redshifts.

When compared to the AzTEC-1.1 mm-selected galaxies within our map (Fig. 13; Miettinen et al., 2017b), it becomes clear that our GISMO 2 mm survey picked out the brightest and highest redshift galaxies among them. Among the six AzTEC-1.1 mm sources within our map with L_{\rm IR}$$>\,10${}^{12.6}\,$L*⊙* and $z>3$, four are detected by our GISMO survey (i.e., 66%). The only AzTEC-1.1 mm-selected galaxy not detected by our survey (S/N$\,\geq\,$3.7) and with $L_{\rm IR}>10^{12.6}\,$L⊙ and $z>3$ are AzTEC1 (a.k.a. AzTEC/C5) and AzTEC3 (a.k.a. AzTEC/C138). However, the GISMO flux density upper limits for these galaxies yield 1.1 mm-to-2 mm lower limits consistent with the rest of the distribution (see the black upward pointing arrows with $S_{1.1{\rm mm}}=9.3$ mJy and 5.9 mJy for AzTEC1 and AzTEC3, respectively, in the bottom-left panel of Fig. 10). These non-detections could thus be simply explained by the inherent incompleteness of our catalog at faint flux densities (see Fig. 4) or by particularly hot dust emission, especially in the case of AzTEC3 (Riechers et al., 2014). In the GISMO map, at the position of AzTEC1, we find a S/N$\,\sim\,$3.4 detection with $S_{\rm 2\,mm}=0.78\pm 0.23\,$mJy, while at the position of AzTEC3 we find a S/N$\,\sim\,$2.4 detection with $S_{\rm 2\,mm}=0.53\pm 0.22\,$mJy.

We measured massive dust content in all GISMO-detected galaxies ($10^{9.3-9.9}\,$M*⊙*; Tab. 3). From the dust and stellar masses of these galaxies, we can calculate the required dust yields per AGB star and per SNe, following Michałowski et al. (2010a, b), i.e.,

[TABLE]

where $N(M_{0}<z<M_{1})$ is the number of stars with masses between $M_{0}$ and $M_{1}$ in the stellar population with a total mass of $M_{\ast}$, i.e.,

[TABLE]

and where IMF$(M)$ is the initial mass function from Chabrier (2003) with $M_{\rm min}=0.15\,$M⊙ and $M_{\rm max}=120\,$M⊙. As in Michałowski et al. (2010a), for AGB stars we assumed $M_{0}=2.5\,$M⊙ and $M_{1}=8\,$M⊙, whereas for SNe we assumed $M_{0}=8\,$M⊙ and $M_{1}=120\,$M⊙. We measure dust yields per AGB star of 1.9$\pm$0.7, 1.3$\pm$0.5 and 0.15$\pm$0.05 M⊙ for GISMO-C1, GISMO-C3 and GISMO-C7, respectively, while theoretical works predict dust yields \lesssim\,$$4\times 10^{-2}\,M⊙ (Michałowski et al., 2010a, b, and reference thererin). We calcutate dust yields per SNe of 6.9$\pm$2.5, 4.9$\pm$1.8 and 0.56$\pm$0.21 M⊙ for GISMO-C1, GISMO-C3 and GISMO-C7, respectively. As for AGB stars, these yields are in tension with theoretical expectations which are \lesssim\,$$1.32\,M⊙ without dust destruction and \lesssim\,$$0.1\,M⊙ with dust destruction (Michałowski et al., 2010a, b, and reference thererin). As pointed out in Michałowski et al. (2010a, b), such unrealistically high dust yields suggest efficient dust production in the interstellar medium of these galaxies. Note that GISMO-C2 is excluded from this analysis because it has no stellar mass estimate. Indeed, its optical/NIR photometry remains very uncertain as it appears very obscured and unfortunately situated in the vicinity of a bright optical/NIR foreground galaxy (see Sect. 2.4 and Fig. 7). GISMO-C6 is also excluded from this analysis because of its most probably wrong redshift identification (see Sect. 2.4).

Assuming a standard gas-to-dust ratio of 100 appropriate for massive systems (Leroy et al., 2011), these dust masses also translate into large gas reservoirs several time more massive than the stellar component of these galaxies. Yet, these gas-rich galaxies with their extreme star-formation activities deplete these reservoirs in $<$1–2 Gyr, in agreement with lower redshift observations (e.g., Tacconi et al., 2018). GISMO-C6 is again the only galaxy with an usually large depletion time scale of $5.7^{+7.3}_{-3.8}$ Gyr.

Note that deriving the gas masses of these galaxies using the methodology advocated in Scoville et al. (2016) (i.e., from our 2 mm deboosted flux densities and their Eq. 16), yields measurements in perfect agreement (within $\sim\,$0.1 dex) with those inferred from the dust model of Draine & Li (2007) and a gas-to-dust ratio of 100. The only exception is again GISMO-C6 for which the methodology of Scoville et al. (2016) leads to $\sim\,$0.4 dex lower gas mass.

3.5 Cosmic infrared luminosity and dust mass densities of galaxies

Using our high-redshift sample, we constrained the bright-end of the infrared luminosity function and cosmic infrared luminosity density at z$$\,\sim\,$$4; and in a pioneer effort, the massive-end of the dust mass function and cosmic dust mass density in galaxies at z$$\,\sim\,$$4 (Fig. 14). To this end, we summed up the contributions of GISMO-C1, -C2, -C3 and -C7, and considered as comoving volume that probed by our survey (250 arcmin2) between $z=3.1$ and $z=4.6$.

Our constraint on the z$$\,\sim\,$$4 infrared luminosity function (i.e., \phi$$\,>\,$$2.5$$\,\times\,$$10^{-6}\Mpc*-3* dex*-1* at $L_{\rm IR}=10^{12.85}\,$L*⊙*) is consistent with that obtained by Gruppioni et al. (2013, upper-left panel of Fig. 14) using far-infrared (100–500$\,\mu$m) observations from the Herschel Space Observatory. This implies that even though our Rayleigh-Jeans 2 mm-selection does not in principle provide a luminosity-limited sample, the population of luminous galaxies with hot dust emission missed by this selection (e.g., AzTEC3) does not dominate at z$$\,\sim\,$$4.

Our study adds also an interesting lower limit to the cosmic infrared luminosity density at z$$\,\sim\,$$4 (lower-left panel of Fig. 14; \rho_{\rm IR}$$\,>\,$$3.3$$\,\times\,$$10^{7} L*⊙* Mpc*-3* at $z=3.9$), which to date remains uncertain owing to the sensitivity limits of current far-infrared surveys (Madau & Dickinson, 2014). Assuming that about half of the cosmic star formation rate density inferred at these redshifts from UV-selected surveys should be seen in form of infrared emission (e.g. Cucciati et al., 2012; Madau & Dickinson, 2014), we estimate that our sample contribute for $\sim\,$20% of the cosmic infrared luminosity density expected at z$$\,\sim\,$$4.

Our pioneer measurement of the z$$\,\sim\,$$4 dust mass function of galaxies (upper-right panel of Fig. 14; \phi$$\,>\,$$5.8^{+4.9}_{-3.3}$$\,\times\,$$10^{-6}\Mpc*-3* dex*-1* at $M_{\rm dust}=10^{9.6}\,$M*⊙*) is consistent with that inferred at z$$\,\sim\,$$2.5 in Dunne et al. (2003) by estimating the dust masses of submillimeter galaxies assuming a flat redshift distribution ranging from z$$\,=\,$$1 to z$$\,=\,$$5. Both studies suggest a mild evolution of the massive-end of the dust mass function of galaxies from z$$\,\sim\,$$0.4 to z$$\,\sim\,$$4.

Finally, we provide a lower limit on the cosmic dust mass density in galaxies at z$$\,\sim\,$$4 (lower-right panel of Fig. 14; \rho_{\rm dust}^{\rm gal.}$$\,>\,$$1.6$$\,\times\,$$10^{4} M*⊙* Mpc*-3* at $z=3.9$). This estimate is consistent with recent theoretical expectations, even though theory predicts a drastic decrease of the cosmic dust mass density in galaxies from z$$\,=\,$$2 to z$$\,=\,$$5 (Gioannini et al., 2017).

4 Conclusions

With the GISMO array at the IRAM 30 m-telescope, we performed the widest deep 2 mm survey to date, reaching a uniform \sigma$$\,\sim\,0.23 mJy beam*-1* sensitivity over an area of \sim\,$$250 arcmin2 in the COSMOS field. Within this map, we detected four sources with a high detection significance (i.e., S/N$\,\geq\,$4.4), corresponding to a low false detection rate per map of only 0.09 sources. Five sources detected with 4.4$\,>\,$S/N$\,\geq\,$3.7 are also added to our catalog, among which 1.65 are supposed to be false detections. With this catalog in hand, we found that:

• •

Combined with the GISMO deep field (\sigma$$\,\sim\,$$0.135 mJy beam*-1* over 13 arcmin2; S14), it provides robust and consistent measurements of the 2 mm number counts over one decade in flux densities. These give critical constraints for current and upcoming galaxy evolution models. For example, while model A of Casey et al. (2018), which represents a “dust-poor” early Universe, best describes our 2 mm number counts, we can begin to rule out subsets of their model B “dust-rich” Universe, whereby only cutoff redshifts $z_{\rm cutoff}<7$ are consistent with our data.

• •

Five sources in our map have counterparts in other deep (sub)mm catalogs available for the COSMOS field (Scott et al., 2008; Aretxaga et al., 2011; Casey et al., 2013; Geach et al., 2017). The redshifts of these sources found in the literature suggest that all but one lie above z$$\,\sim\,3. For these four high-redshift galaxies, their GISMO-2 mm flux densities are consistent with their overall mid-to-far-infrared SEDs, while providing additional constraints on their Rayleigh-Jeans dust emission. These high-redshift galaxies are found to be ultra-luminous infrared galaxies with SFRs in the range 400 – 1200 M${}_{\odot}\,$yr*-1*. They are associated with large dust/gas reservoirs but their extreme SFRs yield gas depletion time scales of $<$1–2 Gyr, in agreement with lower redshift observations (e.g., Tacconi et al., 2018). GISMO-C6 is the only galaxy at odds with this picture, as it is a relatively low redshift galaxy ($z=1$) with a moderate infrared luminosity (L_{\rm IR}$$\,=\,10^{11.6}\,L⊙) and a very cold luminosity-weighted dust temperature ($\sim\,$20 K). The detection of this galaxy by GISMO could be explained by the well-known $L_{\rm IR}-T_{\rm dust}$ selection bias affecting (sub)mm surveys. However, this could also suggest that the redshift association for this source is incorrect.

• •

Comparing the redshift-infrared luminosity distribution of our galaxies to that of the AzTEC-1.1 mm-selected galaxies within our map (Miettinen et al., 2017b; Riechers et al., 2014), we found that our GISMO-2 mm survey is picking out a relatively complete sample ($\sim\,66\%$) of the most luminous (L_{\rm IR}$$>\,10${}^{12.6}\,$L*⊙*) and highest redshift ($z>3$) galaxies among them. This suggests that the selection of bright galaxies at long wavelengths provides the most favorable criterion for finding massive vigorously star-forming high-redshift galaxies, in agreement with recent observations and galaxy evolution models (e.g., Younger et al., 2007, 2009; Strandet et al., 2016; Brisbin et al., 2017; Casey et al., 2018; Zavala et al., 2018). Unfortunately, due to small number statistics, it is not yet possible to fully quantify the observational advantage of 2 mm over 1.1 mm selection for high-redshift studies.

• •

GISMO-C4 is the fourth brightest source in our catalog and is thus very unlikely to be a false detection (P_{\rm f}$$\,=\,$$6.2\%). Yet, it has no (sub)mm counterpart. Such a (sub)mm dropout could be an unidentified very high-redshift galaxies ($z>4$), as suggested by its unusually low 850$\,\mu$m-to-2 mm flux density ratio. This very high-redshift candidate will require future dedicated follow-up with ALMA or NOEMA. Three other sources in our catalog have potentially no (sub)mm counterparts and low 850$\,\mu$m-to-2 mm flux density ratios, however, 1.22 sources are supposed to be false detections.

Our wide GISMO 2 mm survey, combined with the pencil-beam confusion-limited GISMO deep field (S14), has unambiguously demonstrated the advantage of long wavelength surveys for studying the rare massive high-redshift highly star-forming galaxies. Such surveys provide valuable constraints on the yet very uncertain bright-end of the infrared luminosity function and massive-end of the dust mass function at z$$\,\sim\,$$4. However, our 2 mm surveys are still limited by their relatively small sky coverage. The ALMA 2 mm continuum survey in COSMOS (PI Casey, Cycle 6) combined with current and future 2 mm instruments on single dish facilities – like NIKA-II on the IRAM 30m and GISMO-2 and TolTEC on the LMT 50m – will certainly demonstrate further the utility of long-wavelength selection in solving for the relative abundance of dusty star-forming galaxies in the $z>3$ Universe by mapping large areas of sky to sub-mJy depths. They will provide thereby invaluable samples of efficiently-selected high-redshift sources for interferometric spectral scan follow-up with ALMA and NOEMA, enabling the study of dust production within the first $\sim$2 Gyr of cosmic time.

We would like to thank the referee for their comments which have helped improving the paper. BM, AK, EJA and FB acknowledge support from the Collaborative Research Centre 956, sub-project A1, funded by the Deutsche Forschungsgemeinschaft (DFG). We would like to thank Carsten Kramer, Santiago Navarro, David John, Albrecht Sievers, and the entire IRAM Granada staff for their support during the instrument installation and observations. CMC thanks the University of Texas at Austin College of Natural Sciences and NSF grants AST-1714528 and AST-1814034 for support. D.R. acknowledges support from the National Science Foundation under grant number AST-1614213. Based on observations carried out under project number 247-11, 227-12, 242-13, 117-14 and 232-15 with the IRAM 30-meter telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). This work was also supported through NSF ATI grants 1020981 and 1106284.