The VLA-COSMOS 3 GHz Large Project: Continuum data and source catalog release

TL;DR

This paper presents the deepest and largest 3 GHz radio continuum survey of the COSMOS field, providing a detailed catalog of over 10,800 sources with high positional accuracy and spectral information, advancing radio astronomical studies.

Contribution

The study delivers the first large-scale, high-resolution 3 GHz radio source catalog for the COSMOS field, combining extensive observations with advanced imaging and analysis techniques.

Findings

Catalog of 10,830 sources down to 5 sigma detection threshold

Median spectral index of -0.7 derived from combined data

Astrometry accurate to 0.01 arcseconds for bright sources

Abstract

We present the VLA-COSMOS 3 GHz Large Project based on 384 hours of observations with the Karl G. Jansky Very Large Array (VLA) at 3 GHz (10 cm) toward the two square degree Cosmic Evolution Survey (COSMOS) field. The final mosaic reaches a median rms of 2.3 uJy/beam over the two square degrees at an angular resolution of 0.75". To fully account for the spectral shape and resolution variations across the broad (2 GHz) band, we image all data with a multiscale, multifrequency synthesis algorithm. We present a catalog of 10,830 radio sources down to 5 sigma, out of which 67 are combined from multiple components. Comparing the positions of our 3 GHz sources with those from the Very Long Baseline Array (VLBA)-COSMOS survey, we estimate that the astrometry is accurate to 0.01" at the bright end (signal-to-noise ratio, S/N_3GHz > 20). Survival analysis on our data combined with the VLA-COSMOS 1.4~GHz Joint Project catalog yields an expected median radio spectral index of alpha=-0.7. We compute completeness corrections via Monte Carlo simulations to derive the corrected 3 GHz source counts. Our counts are in agreement with previously derived 3 GHz counts based on single-pointing (0.087 square degrees) VLA data. In summary, the VLA-COSMOS 3 GHz Large Project simultaneously provides the largest and deepest radio continuum survey at high (0.75") angular resolution to date, bridging the gap between last-generation and next-generation surveys.

Figures

Tables

Table 1. Table 1: Amplitude of the phase calibrator (J1024-0052) in each SPW for all observing blocks

SPW	Frequency	Mean flux	Median flux	Standard
	(GHz)	density (Jy)	density (Jy)	deviation (Jy)
1	2.060	0.739	0.735	0.029
4	2.444	0.707	0.704	0.025
5	2.572	0.700	0.696	0.023
6	2.700	0.684	0.680	0.023
7	2.828	0.668	0.665	0.024
8	2.956	0.652	0.648	0.026
9	3.084	0.645	0.642	0.020
10	3.212	0.635	0.632	0.020
11	3.340	0.625	0.622	0.021
12	3.468	0.615	0.611	0.020
13	3.596	0.603	0.600	0.021
14	3.724	0.539	0.579	0.129
15	3.852	0.525	0.569	0.153
16	3.980	0.535	0.566	0.124

Table 2. Table 2: Completeness and bias correction factors for the VLA-COSMOS 3 GHz catalog as a function of flux density

Flux density	Completeness and bias	Error
($\mu$Jy)	correction factor (Ccompl)
10.4	0	-
11.0	0.08	0.01
11.6	0.40	0.02
12.3	0.55	0.03
13.0	0.58	0.03
13.8	0.56	0.04
14.6	0.57	0.03
15.5	0.57	0.04
16.4	0.57	0.04
17.3	0.57	0.04
18.4	0.59	0.04
19.4	0.56	0.04
21.1	0.68	0.05
23.3	0.73	0.04
25.8	0.82	0.05
28.6	0.85	0.05
31.7	0.89	0.06
35.1	0.91	0.07
38.8	0.94	0.06
43.0	0.92	0.06
47.6	0.95	0.07
53.9	0.94	0.05
62.4	0.95	0.07
72.2	0.93	0.06
83.5	0.92	0.08
96.7	0.97	0.08
$>{100}^a$	1.00a	0.05a

Table 3. Table 3: False detection probability as a function of S/N and flux density in the COSMOS two square degree field

S/N	Fraction
5.05	0.24
5.15	0.15
5.25	0.11
5.35	0.09
5.45	0.06
5.55	0.03
5.65	0.02
5.75	0.03
5.85	0.01
5.95	0.02
6.05	0.01
6.15	0.01
6.25	0.02
6.35	0.01
6.45	0.00
6.55	0.00
6.65	0.01
6.75	0.01
6.85	0.00
6.95	0.01

Table 4. Table 4: Radio source counts at 3 GHz within the COSMOS two square degree field, normalized to Euclidean geometry

Flux density	Countsa	Errorb	N	Correction
(mJy)	(Jy1.5sr-1)	(Jy1.5sr-1)		factora
0.011	0.866	0.068	631	3.27
0.013	0.952	0.056	1109	1.64
0.015	1.10	0.078	991	1.70
0.018	1.16	0.094	849	1.67
0.020	1.40	0.11	888	1.54
0.024	1.38	0.086	811	1.33
0.028	1.47	0.10	780	1.18
0.032	1.57	0.12	702	1.12
0.037	1.59	0.13	587	1.08
0.045	1.58	0.13	753	1.07
0.059	1.53	0.12	505	1.05
0.076	1.67	0.14	366	1.08
0.098	1.60	0.17	250	1.03
0.13	1.65	0.15	181	1.00
0.16	1.62	0.17	121	1.00
0.21	1.82	0.21	93	1.00
0.27	1.93	0.25	67	1.00
0.35	2.03	0.31	48	1.00
0.78	3.95	0.37	159	1.00
3.0	10.4	1.5	56	1.00
11	37.4	7.4	27	1.00
44	72.3	28	7	1.00

Table 5. Table 5: Catalog sample page

ID	NAME	RA	RA_ERR	DEC	DEC_ERR	FLUX	FLUX_ERR	RMS	SNR	NPIX	RES	MULTI
		(J2000, deg)	(arcsec)	(J2000, deg)	(arcsec)	($\mu$Jy)	($\mu$Jy)	($\mu$Jy beam-1)
78	COSMOSVLA3 J095709.33+020940.7	149.288886	0.01	2.161331	0.01	13400.0	670.0	28.7	385.0	104	1	0
1110	COSMOSVLA3 J095709.83+015457.4	149.290996	0.021	1.915946	0.021	1190.0	64.0	22.1	29.0	102	1	0
5144	COSMOSVLA3 J095710.49+013644.7	149.29372	0.065	1.612418	0.065	145.0	19.0	17.4	8.32	20	0	0
3749	COSMOSVLA3 J095710.53+025132.2	149.29389	0.053	2.858966	0.053	345.0	32.0	26.8	10.4	31	1	0
10099	COSMOSVLA3 J095710.57+022657.0	149.294072	0.1	2.449174	0.1	92.0	18.0	17.4	5.3	11	0	0
4979	COSMOSVLA3 J095710.88+020929.8	149.295339	0.076	2.158297	0.076	2820.0	140.0	21.4	7.13	332	1	0
3366	COSMOSVLA3 J095711.09+023031.3	149.296231	0.048	2.508722	0.048	318.0	27.0	21.5	11.4	35	1	0
8753	COSMOSVLA3 J095711.15+021104.2	149.296459	0.095	2.184514	0.095	121.0	22.0	21.2	5.68	18	0	0
4046	COSMOSVLA3 J095711.64+021236.8	149.298529	0.055	2.210248	0.055	223.0	26.0	22.6	9.83	22	0	0
6546	COSMOSVLA3 J095711.67+021401.3	149.298655	0.077	2.233706	0.077	135.0	21.0	19.3	7.02	18	0	0
…
10942	COSMOSVLA3 J100034.76+014635.7	150.144846	-99.0	1.776607	-99.0	374.0	-99.0	2.4	-99.0	828	1	1
9229	COSMOSVLA3 J100034.78+025027.4	150.144947	0.1	2.840947	0.1	25.6	2.7	2.4	5.18	25	1	0
8777	COSMOSVLA3 J100034.80+021421.1	150.145009	0.095	2.239195	0.095	13.0	2.4	2.27	5.7	14	0	0
489	COSMOSVLA3 J100034.81+025515.6	150.145042	0.013	2.921023	0.013	229.0	13.0	3.66	62.5	49	0	0
439	COSMOSVLA3 J100034.83+014247.1	150.145157	0.013	1.713091	0.013	162.0	9.0	2.27	71.3	51	0	0
2850	COSMOSVLA3 J100034.94+020234.9	150.145599	0.042	2.043054	0.042	29.9	2.8	2.29	13.1	29	0	0
2439	COSMOSVLA3 J100034.99+024524.5	150.145815	0.038	2.756821	0.038	32.7	2.8	2.22	14.7	29	0	0
2837	COSMOSVLA3 J100035.00+024614.3	150.145836	0.041	2.770643	0.041	39.6	2.9	2.18	13.3	41	1	0
4092	COSMOSVLA3 J100035.05+024154.6	150.146065	0.055	2.69852	0.055	28.3	2.6	2.23	9.93	31	1	0
4148	COSMOSVLA3 J100035.07+020350.5	150.14614	0.057	2.064044	0.057	23.4	2.7	2.43	9.64	23	0	0
…
5351	COSMOSVLA3 J100345.81+015420.5	150.940886	0.067	1.905719	0.067	117.0	16.0	14.6	8.05	18	0	0
5225	COSMOSVLA3 J100346.48+023458.3	150.943703	0.067	2.582882	0.067	127.0	17.0	15.5	8.14	21	0	0
8447	COSMOSVLA3 J100346.52+022031.2	150.943867	0.11	2.342001	0.11	290.0	20.0	14.4	5.04	48	1	0
10677	COSMOSVLA3 J100346.56+015500.1	150.944027	0.11	1.916707	0.11	75.7	15.0	14.9	5.1	11	0	0
6533	COSMOSVLA3 J100346.63+022415.8	150.944306	0.076	2.404402	0.076	143.0	16.0	14.2	7.09	28	1	0
2195	COSMOSVLA3 J100347.12+022510.6	150.946345	0.036	2.419631	0.036	548.0	31.0	15.3	15.6	83	1	0
10084	COSMOSVLA3 J100347.27+020117.6	150.94699	0.1	2.021575	0.1	73.1	15.0	14.1	5.18	16	0	0
9695	COSMOSVLA3 J100347.51+023539.7	150.947979	0.1	2.594384	0.1	93.2	18.0	17.6	5.3	10	0	0
10652	COSMOSVLA3 J100347.80+024951.9	150.949175	0.11	2.831096	0.11	80.7	17.0	16.1	5.01	14	0	0
10484	COSMOSVLA3 J100348.53+021102.3	150.952242	0.11	2.183987	0.11	82.3	17.0	16.3	5.05	13	0	0

Equations

$$\frac{S_t}{S_p}=\frac{\sqrt{{\theta }_M^2+{\theta }_b^2}\sqrt{{\theta }_m^2+{\theta }_b^2}}{{\theta }_b^2},$$

$$\frac{S_t}{S_p}=\sqrt{\frac{{\theta }_M^2+{\theta }_b^2}{{\theta }_b^2}}.$$

$$\frac{S_t}{S_p}=\frac{{\theta }_M^2+{\theta }_b^2}{{\theta }_b^2}.$$

$${\theta }_{min}=k_1{e}^{-(S_t/k_2{)}^2},$$

Full text

1 Introduction

One of the main quests in modern cosmology is understanding the formation of galaxies and their evolution through cosmic time. In the past decade it has been demonstrated that a panchromatic, X-ray to radio, observational approach is key to develop a consensus on galaxy formation and evolution (e.g., Dickinson et al. 2003; Scoville et al. 2007; Driver et al. 2009, 2011; Koekemoer et al. 2011; Grogin et al. 2011). In this context, the radio regime offers an indispensable window toward star formation and supermassive black hole properties of galaxies as radio continuum emission i) provides a dust-unbiased star formation tracer at high angular resolution (e.g., Condon 1992; Haarsma et al. 2000; Seymour et al. 2008; Smolčić et al. 2009b; Karim et al. 2011), and ii) directly probes those active galactic nuclei (AGN) that are hosted by the most massive quiescent galaxies and deemed crucial for massive galaxy formation (e.g., Croton et al. 2006; Bower et al. 2006; Best et al. 2006; Evans et al. 2006; Hardcastle et al. 2007; Smolčić et al. 2009a; Smolčić 2009; Smolčić & Riechers 2011; Smolčić et al. 2015a).

In recent decades, radio interferometers, such as the Karl G. Jansky Very Large Array (VLA), Australia Telescope Compact Array (ATCA), and Giant Meterwave Radio Telescope (GMRT), have surveyed fields of different sizes (ranging from tens of square arcminutes to thousands of square degrees), depths (microjansky to Jansky), and multiwavelength coverage (e.g., Becker et al. 1995; Condon et al. 1998; Ciliegi et al. 1999; Georgakakis et al. 1999; Bock et al. 1999; Prandoni et al. 2001; Condon et al. 2003; Hopkins et al. 2003; Schinnerer et al. 2004; Bondi et al. 2003, 2007; Norris et al. 2005; Schinnerer et al. 2007, 2010; Afonso et al. 2005; Tasse et al. 2007; Smolčić et al. 2008a; Owen & Morrison 2008; Miller et al. 2008, 2013; Owen et al. 2009; Condon et al. 2012; Smolčić et al. 2014; Hales et al. 2014). These past surveys have shown that deep observations at high angular resolution ($\lesssim 1\arcsec$) with exquisite panchromatic coverage are critical to comprehensively study the radio properties of the main galaxy populations, avoiding cosmic variance with large area coverage (e.g., Padovani et al. 2009; Padovani 2011; Smolčić et al. 2008a, 2009b, 2009a; Smolčić 2009; Smolčić & Riechers 2011; Seymour et al. 2008; Bonzini et al. 2012, 2013). In this context, large area surveys down to unprecedented depths are planned with new and upgraded facilities (e.g., VLA, Westerbork, ASKAP, MeerKAT, and SKA; e.g., Jarvis 2012; Norris et al. 2011, 2013, 2015; Prandoni & Seymour 2015). Figure 1 shows the $1\sigma$ sensitivity of each survey as a function of the area covered for past, current, and future radio continuum surveys. The VLA-Cosmic Evolution Survey (COSMOS) 3 GHz Large Project bridges the gap between past and future radio continuum surveys by covering an area as large as two square degrees down to a sensitivity reached to date only for single pointing observations. This allows for individual detections of $>10,000$ radio sources, further building on the already extensive radio coverage of the COSMOS field at 1.4 GHz VLA (VLA-COSMOS Large, Deep and Joint projects; Schinnerer et al. 2004, 2007, 2010), 320 MHz VLA (Smolčić et al., 2014), 325 MHz and 616 MHz GMRT data (Karim et al., in prep.; Brady et al., in prep.), 6 GHz VLA (Myers et al., in prep.), and the deep multiwavelength X-ray to mm photometry (Scoville et al. 2007; Koekemoer et al. 2007; Hasinger et al. 2007; Capak et al. 2007; Sanders et al. 2007; Bertoldi et al. 2007; Elvis et al. 2009; Ilbert et al. 2013; McCracken et al. 2012; Scott et al. 2008; Aretxaga et al. 2011; Smolčić et al. 2012; Miettinen et al. 2015; Civano et al. 2016; Laigle et al. 2016, Capak et al., in prep.) and more than 97,000 optical spectroscopic redshifts (Salvato et al., in prep.; zCOSMOS, Lilly et al. 2007, 2009; Trump et al. 2007; Prescott et al. 2006; Le Fèvre et al. 2015; Aihara et al. 2011; Nagao et al., priv. comm.). This further makes the survey part of one of the richest multiwavelength data sets currently available.

Radio continuum surveys at 3 GHz with the upgraded VLA are still sparse in the literature. Condon et al. (2012) performed single-pointing observations targeting the Lockman hole for 50-hours on-source with the VLA in C-array configuration. The observations resulted in a confusion-limited map with an $rms$ of $1\leavevmode\nobreak\ \mu$Jy beam*-1*. Based on this they constrained the counts of discrete sources in the $1-10\leavevmode\nobreak\ \mu$Jy range via a $P(D)$ analysis. A more complex $P(D)$ analysis using the same data was applied by Vernstrom et al. (2014) who probed the counts down to $0.1\leavevmode\nobreak\ \mu$Jy. Both results are qualitatively in agreement with the already well-known flattening of the radio source counts (normalized to the $N(S)\propto S^{-3/2}$ of a static Euclidian space) below flux densities of $S_{\mathrm{1.4GHz}}\approx 1$ mJy, and a further decrease of the counts with decreasing flux density below $S_{\mathrm{1.4GHz}}\approx 60\leavevmode\nobreak\ \mu$Jy. Such a shape of radio source counts is expected owing to the cosmic evolution of galaxy populations (e.g., Hopkins et al. 2000; Wilman et al. 2008; Béthermin et al. 2012), but this shape is contrary to that obtained based on i) the previous Lockman hole observations at 1.4 GHz (Owen & Morrison, 2008), and ii) a comparison of the sky brightness temperature measured by the ARCADE 2 experiment (Fixsen et al., 2009) with that derived from the integral of the observed radio source counts (Vernstrom et al., 2011). The latter results instead point to a rise of the counts with decreasing flux density at these levels. To investigate this further, we here derive the radio source counts using our VLA-COSMOS 3GHz Large project data, yielding the deepest radio counts derived to date based on direct source detections.

In Sect. 2 we describe the VLA 3 GHz observations, calibration, and imaging. We present the catalog extraction in Sect. 3, an analysis of the radio spectral indices in Sect. 4, the radio source count corrections in Sect. 5 and the radio source counts in Sect. 6. We summarize our products and results in Sect. 7. We define the radio spectral index $\alpha$ as $S_{\nu}\propto\nu^{\alpha}$, where $S_{\nu}$ is flux density at frequency $\nu$.

2 Observations and data reduction

2.1 Observations

A total of 384 hours of observations toward the COSMOS field were taken in S band using the S 3s full width setup covering a bandwidth of 2048 MHz centered at 3 GHz, and separated into 16 $128$ MHz-wide spectral windows (SPWs hereafter), with full polarization, and a 3s signal-averaging time. The observations were taken from November 2012 to January 2013, June to August 2013, and February to May 2013 in A-array (324 hours) and C-array configurations (60 hours; Legacy ID AS1163). Sixty-four pointings, separated by $10\arcmin$ in right ascension (RA) and declination (DEC), corresponding to two-thirds of the half-power beam width (HPBW) at the central frequency of 3 GHz, were chosen to cover the full two square degree COSMOS field. Three sets of 64 pointings in such a grid were used to achieve a uniform $rms$ over the field; this resulted in a total of 192 pointings (shown in Fig. 2). The first set of pointings is nominal, the second is shifted by $5\arcmin$ in RA and DEC, while the third set is shifted by $-5\arcmin$ in RA. Observing runs of 5 and 3 hours in length were conducted. In each observing run J1331+3030 was observed for flux and bandpass calibration for about 3-5 minutes on-source (J0521+166 was used only for the first day of observations) at the end of every run, J1024-0052 was observed every 30 minutes for 1m 40s on-source for gain and phase calibration, while the source J0713+4349 was observed for 5 minutes on-source at the beginning of each run for polarization leakage calibration. During the five hour observing runs each pointing was visited twice, while the order of the pointing coverage blocks during the fixed 5-hour observing blocks was changed between the different observing runs to optimize the $uv$ coverage. During the 3-hour observing blocks each pointing was visited once, and a good $uv$ coverage was assured via dynamic scheduling. Typically, 26 antennas were used during each observing run. The A-array configuration observations were mostly conducted under good to excellent weather conditions. The C-array configuration observations were partially affected by poor weather conditions (Summer thunderstorms), yielding on some days up to 30% higher $rms$ than expected based on the VLA exposure calculator.

2.2 Calibration

Calibration of the data was performed via the Astronomical Image Processing System (AIPS) based data reduction pipeline AIPSLite (Bourke et al., 2014) developed for the Caltech-National Radio Astronomy Observatory (NRAO) Stripe 82 Survey (Mooley et al., 2016). This pipeline was adapted for the VLA-COSMOS 3 GHz Large Project (as described below) and it follows, in general, the procedures outlined in Chapter E of the AIPS Cookbook111http://www.aips.nrao.edu/cook.html.

In brief, the data are first loaded with the Obit222http://www.cv.nrao.edu/~bcotton/Obit.html task BDFIn. Band edges, and to a larger extent IF edges, were then flagged with the task UVFLG. SPWs 2 and 3, found to be irreparably corrupted by radio frequency interference (RFI) in all observations (see Fig. 3), were entirely flagged using the task UVFLG. After flagging, FRING, BPASS, SETJY, CALIB, GETJY, and CLCAL tasks were used to derive the delay, bandpass, and complex gain solutions. Polarization calibration was performed using the tasks RLDLY, PCAL, and RLDIF as detailed in Sect. 7 of Chapter E in the AIPS Cookbook. The task RFLAG was used to flag all target pointings and the flags were applied using the UVCOP task. The derived calibration was applied and the calibrated data set was produced with the SPLAT task. Finally, the calibrated $uv$ data were saved to disk using the task FITTP. During the pipeline process several diagnostic plots were generated to assess the quality of the calibration: bandpass solutions, antenna gains as a function of time, calibrated spectrum of the gain calibrator, and calibrated amplitude versus phase plots of the gain calibrator per pointing. In Table 1 we list the statistics for the amplitude of the phase calibrator in each SPW for all observing blocks. The average amplitude scatter around the mean is typically $\leavevmode\nobreak\ 2-3$%, with the exception of the highest frequency SPWs, for which it is higher than 10%. The combined typical scatter around the mean is $\sim 5$%. This assures a good flux calibration.

The highest frequency SPWs marked 14, 15, and 16 have low amplitude RFI and the phases are significantly affected for some observations. The C-array configuration data at the upper end of the S band are mostly unusable due to this RFI. These data have been manually flagged, and we additionally ran the RFLAG task on the rest of the C-array configuration data to further remove bad data and extend flags in frequency and time. The A-array configuration data for these SPWs are generally good. Our imaging tests show that the data from these SPWs generally improves the sensitivity, but limits the dynamic range for certain pointings. Looking at the overall imaging performance, we decided to retain these SPWs. Despite the data drop-outs, the median flux density values of the phase calibrator (J1024-0052; Table 1) are consistent with the spectral parameters inferred from the other SPWs. Through our tests we find that, in the majority of observations, RFI adversely affects the system temperature measurements, and hence we have left out the correction for the system temperature from the calibration process.

At this point the pipeline diverges in two directions to: i) image the target fields and ii) produce and export a calibrated data set in preparation for mosaicking. To image the target fields, they were split out with calibration applied (using the task SPLIT). The fields were then further auto-flagged (using the task RFLAG), imaged (using the task IMAGR), and exported (using the task FITTP) in parallel. The calibrated data set was generated by applying RFLAG and imaging the target fields, including applying flags (using the task UVCOP), calibration (using the task SPLAT), and exporting the $uv$ data and maps (using the task FITTP).

The pipeline performance and output were tested by i) manually reducing separate blocks of VLA-COSMOS observations and comparing the results with the pipeline output, and ii) comparing the output to the Common Astronomy Software Applications (CASA333CASA is developed by an international consortium of scientists based at the NRAO, the European Southern Observatory (ESO), the National Astronomical Observatory of Japan (NAOJ), the CSIRO Australia Telescope National Facility (CSIRO/ATNF), and the Netherlands Institute for Radio Astronomy (ASTRON) under the guidance of NRAO. See http://casa.nrao.edu; (McMullin et al., 2007))-based NRAO reduction pipeline for randomly selected data taken in the A- and C-array configurations. No obvious differences were found. As the pipeline used here was tailored specifically to the COSMOS field (e.g., it includes polarization calibration), after this verification it was further applied to the remaining VLA-COSMOS data sets.

The calibrated $uv$ data sets output by the pipeline for each observing block were first run through the AIPS task UVFIX to assure accurately computed positions. We note that applying UVFIX at the end of calibration has the same effect as applying it at the beginning of calibration. They were then further processed in CASA by clipping each calibrated $uv$ data set in amplitude (above 0.4 Jy) using the task FLAGDATA444In total, about 30-35% of the data were flagged (using the tasks RFLAG and FLAGDATA)., splitting the individual pointings using the task SPLIT, and concatenating all existing observations of the same pointing using the task CONCAT. The concatenated $(u,v)$ data for each pointing were then imaged prior to being combined into the final mosaic, as described in detail in the next section.

2.3 Self-calibration, imaging, and mosaicking

To image our data we used the multiscale multifrequency synthesis (MSMF) algorithm developed by Rau & Cornwell (2011) and implemented in CASA. This method uses the entire 2 GHz bandwidth at once to calculate the monochromatic flux density at 3 GHz and a spectral index between 2 and 4 GHz. After extensive testing of various imaging methods (see Novak et al. 2015) we settled for the MSMF method as it allows for a combination of the best possible resolution, $rms$, and image quality. Because of the large data volume, joint deconvolution was not practical and we imaged each pointing individually and then combined them into a mosaic in the image plane.

We found sources that were bright enough (peak surface brightness higher than 5 mJy beam*-1*) to allow for self-calibration in 44 out of 192 pointings. To prevent artifacts affecting the model used for self-calibration small clean masks were centered around bright sources. An integration time of 3 min, which roughly corresponds to one scan length, was used to obtain phase-gain solutions for these pointings (i.e., only the phase part of the complex gain was solved for and applied). It was typically not possible to find a solution using self-calibration for 10% of the data with the fraction increasing to 20% for a few pointings, which was the maximum value we allowed. We applied gain solutions to the $uv$ data but did not apply the flags calculated in the self-calibration process as that usually increased the noise in the map. For the remaining pointings we applied phase gains obtained by self-calibrating the phase-calibrator J1024-0052, as it further reduced artifacts and sidelobes around brighter sources as illustrated in Fig. 4.

We used the CLEAN task with Briggs weighting scheme for gridding of visibilities with a robust parameter of 0.5 to obtain the best compromise between the resolution and the noise. Two Taylor terms (nterms=2: TT0 and TT1) were used for multifrequency synthesis, which allows the reconstruction of the total intensities and spectral slopes (Rau & Cornwell, 2011). Each pointing was tapered with its own Gaussian to achieve a circular beam, where the difference between the major and the minor axis is 3% at maximum (see Fig. 5). Prior to this step the beam was slightly elliptical, but the position angle changed considerably between different pointings. A cyclefactor of 3 was applied for a more robust deconvolution and to prevent artifacts in the map possibly caused by sidelobe intersections. Widefield imaging was necessary to produce correct astrometry far from the pointing center and we used 128 projection planes. We cleaned on three additional spatial scales corresponding to $2\times$, $5\times,$ and $10\times$ the synthesized beam size to better handle extended sources such as radio jets and lobes. A gain parameter of 0.3 was used to speed up this multiscale algorithm. Each pointing map was set to 8,000 pixels on-the-side with a pixel size of $0.2\times 0.2$ arcsec2. Cleaning was performed down to $5\sigma$ in the entire map and further down to 1.5$\sigma$ using tight masks around sources. These masks were defined manually across the entire observed field by visually inspecting the mosaic555A preliminary mosaic was generated with pointings cleaned down to $5\sigma$ and then used to define cleaning masks. Masks were usually circles with 0.7$\arcsec$ radius, but they were modified where necessary to accommodate larger (resolved) extended sources. It was not necessary to set clean boxes around known strong sources outside of the imaged area. . Synthesized beam size variations between different pointings were about $0.03\arcsec$, which was small enough to allow restoration of every cleaned pointing to an average circular beam of 0.75$\arcsec$. Finally, each pointing was corrected for the frequency-dependent primary beam response down to a value of 20% (corresponding to a radius of 10.5$\arcmin$) using the WIDEBANDPBCOR task. The noise level in the phase center of an individual pointing was usually around 4-5 $\mu$Jy beam*-1*.

To construct the mosaic of all pointings, we used our custom IDL procedure combined with the AIPS task FLATN to carry out noise weighted addition of the individually imaged pointings. Every pixel in the sum was weighted by the inverse square of the local $rms$, which was determined in the pointing itself via the AIPS task RMSD (see below). We mosaicked both Taylor terms individually using the noise weights calculated from the total intensity maps. The 3 GHz continuum mosaic is shown in Fig. 6, where we overplot Gaussian fits to the pixel surface brightness distributions across the mosaic. Cutouts of several extended sources and a mosaic zoom-in are presented in Fig 7. The visibility function showing the covered area at a given $rms$ is presented in Fig. 8. In summary, the final mosaic has a resolution of $0.75\arcsec$, with a median $rms$ of 2.3 $\mu$Jy beam*-1* over the COSMOS 2 square degrees.

3 Cataloging

3.1 Extracting source components

To extract source components from the VLA-COSMOS MSMF mosaic and catalog their properties we employed blobcat developed by Hales et al. (2012). Extractor blobcat uses the flood fill algorithm to find islands of pixels (blobs) above a certain signal-to-noise ratio (S/N) threshold. The local noise map used to evaluate the S/N at each pixel position was created from the total intensity mosaic with the AIPS task RMSD with a circular mesh size of 100 pixels. Once blobcat locates islands, it measures the peak surface brightness ($S_{p}$) by fitting a two-dimensional (2D) parabola around the brightest pixel, while the total flux density ($S_{t}$) is obtained by summing up the pixel values inside the island and dividing the sum by the beam size in pixels. In the next step blobcat takes into account a small positive peak surface brightness bias created by the presence of noise peaks in the map and also corrects for a negative integrated surface brightness bias caused by the finite island size used for integration. We used default parameters when running blobcat (as Hales et al. 2012 ran extensive simulations to optimize them; see also Hales et al. 2014), where the required size of a blob is at least 3 pixels in RA and 3 pixels in DEC. This was necessary to detect low S/N sources, which would have otherwise been missed owing to our relatively coarse pixel grid. With this setup we recovered 10,899 radio source components with local S/N greater or equal to 5 across the entire observed area. As detailed in Sect. 3.3, 67 components have been merged into unique, multicomponent sources resulting in a total of 10,830 radio sources.

3.2 Resolved versus unresolved sources

In order to determine whether our identified source components are resolved (i.e., extended, larger than the synthesized beam) we make use of the ratio between total flux density ($S_{t}$) and peak surface brightness ($S_{p}$) as this is a direct measure of the extension of a radio source. The flux densitites were computed by blobcat as described in the previous section. For a perfect Gaussian unresolved source, the peak surface brightness in Jy beam*-1* equals the integrated flux density in Jy or $S_{t}/S_{p}=1$. The extension of a radio source increases its total flux density when compared to its peak surface brightness, however, background noise can lower the total flux density (see Bondi et al. 2003). Therefore, in Fig. 9 we plot the ratio between the total flux density and the peak surface brightness as a function of the S/N (=$S_{p}/rms$) for all 10,899 components in the catalog. To select the resolved components, we determined the lower envelope of the points in Fig. 9, which contains 95% of the components with $S_{t}<S_{p}$ and mirrored it above the $S_{t}/S_{p}=1$ line (upper envelope in Fig. 9). The shape of the envelope was chosen following Bondi et al. (2008) and the fit to our data is given as $S_{t}/S_{p}=1+6\times(\mathrm{S/N})^{-1.44}$. We consider the 3,975 components above the upper envelope as resolved. These resolved components were flagged in the catalog. For the unresolved components the total flux density was set equal to the peak surface brightness in the catalog.

3.3 Multicomponent sources

Large sources with diffuse structures, such as radio galaxies (see Fig. 7) or resolved star-forming disks, can be listed in a component catalog as multiple entries. This can happen for example if there is no significant radio emission between the two radio lobes, or if the local $rms$ noise is overestimated because of large-scale faint radio emission, which affects the ability of blobcat to properly detect a contiguous blob. We identified 10,899  components in our mosaic, as described above. In order to generate a source catalog rather than a source component catalog, we aimed to identify such sources and convert the multiple entries into one entry that described the entire source, i.e., listing its proper total flux density and position. For this purpose we visually inspected over 2,500 components. The inspected sample was a combination of the i) brightest 2,500 components, ii) all known multicomponent sources that were identified and listed in the 1.4 GHz joint catalog (126 components), and iii) sources with $R_{\mathrm{EST}}>1+30\times(\mathrm{S/N})^{-1}$ (351 components). The $R_{\mathrm{EST}}$ parameter is a size estimate reported by blobcat, which can be used to find sources with non-Gaussian morphology; see Hales et al. (2012, 2014) for more details. Following the procedure already applied to the VLA-COSMOS 1.4 GHz survey sources (Schinnerer et al., 2007) these components were visually inspected with respect to the near-infrared (NIR) images, i.e., the $z^{++}YJHK$ stacked maps (Laigle et al., 2016). In total, we identifed 67 multicomponent sources. As for the previous VLA-COSMOS survey catalogs, we computed their total flux densities using the AIPS task TVSTAT in the area encompassed by $2\sigma$ contours, where $\sigma$ is the local $rms$ measured as the average $rms$ from a 100-300 pixel wide area around the source, ensuring that the $rms$ is not biased by the influence of the strong sources. The source position is then taken to be the radio core or optical counterpart position (if identifiable) or the luminosity weighted mean. In our catalog we then excluded all the components combined into the multicomponent sources, and listed instead the multicomponent source with the above-defined position and total flux density, setting all other cataloged values to -99. A further column multi was added designating the multicomponent sources (multi=1 for a multicomponent source, and multi=0 for a single-component source). We note that the number of multicomponent sources is smaller than that identified in the shallower VLA-COSMOS 1.4 GHz survey. This is due to the higher frequency of breaking-up large sources into multiple components within the latter as it used the AIPS Search and Destroy source finding algorithm, when compared to the performance of the blobcat algorithm used here. A full assessment of large sources in the survey will be presented by Vardoulaki et al. (in prep.).

3.4 Astrometric accuracy

To assess our astrometric accuracy at the bright end we have compared the positions of 443 sources at 3 GHz with S/N $>$ 20, also detected in the Very Long Baseline Array (VLBA)-COSMOS 1.4 GHz survey (PI: Middelberg; N. Herrera Ruiz et al., in prep.). The results, shown in Fig. 10, yield an excellent agreement with a mean offset of $0.01\arcsec$ in $\Delta$ RA and $0.00\arcsec$ in $\Delta$ DEC and a standard deviation of 0.01$\arcsec$ for both. We note that we did not correct the catalog entries for the $0.01\arcsec$ offset in $\Delta$ RA. We took the standard deviation value (0.01$\arcsec$) as the calibration error in RA and DEC to compute the positional uncertainties for our sources using the equations reported in Hales et al. (2012). We note that these are estimated to be accurate for point-sources, but likely underestimated for resolved sources (see Hales et al. 2012 and references therein for details).

3.5 Bandwidth smearing

Owing to the finite bandwidth of the antenna receiver, bandwidth smearing (BWS) occurs and radially smears peak surface brightness while conserving the integrated flux density. The effect is a function of distance from the phase center in a given pointing while it reaches a constant smearing value in the combined mosaic (see, e.g., Bondi et al. 2008). Although the bandwidth of the antenna receiver is large ($\sim 4$ GHz), the relevant bandwidth for the smearing effect is only the $2$ MHz channel width used to image the data.

To empirically test BWS in our data, we selected 106 point-like ($0.9\leq S_{\mathrm{t}}/S_{\mathrm{p}}\leq 1.1$) radio sources with S/N$>$200. Since each source can be observed in up to 11 neighboring pointings, we can compare the peak surface brightnesses obtained in various pointings ($S_{\mathrm{P}}$) relative to the peak surface brightness retrieved from the mosaic ($S_{\mathrm{M}}$) as a function of distance from the pointing center. If our data were affected by BWS, $S_{\mathrm{P}}/S_{\mathrm{M}}$ would exhibit a declining trend with increasing distance from the pointing center. This surface brightness ratio, obtained by fitting an inverted parabola at the 106 bright source positions in the individual pointings, and the mosaic is shown in the top panel of Fig. 11. The median ratio stays constant ($S_{P}/S_{M}\approx 1$) across all distance ranges, with increasing scatter toward higher distances where the noise is amplified by the primary beam correction. This demonstrates that there are no empirical bandwidth-smearing issues. This is also in accordance with theoretical expectations. A theoretical prediction for BWS can be made using the Condon et al. (1998) equation [12] for the reduction of peak response $I/I_{0}\approx 1/\sqrt{1+0.46\beta^{2}}$, where $\beta=(\Delta\nu/\nu_{0})\times(\theta_{0}/\theta_{HPBW})$ equals fractional bandwidth times offset in synthesized beam-widths. Using the VLA channel width $\Delta\nu=2$ MHz, central frequency $\nu_{0}=3$ GHz, distance from the phase center $\theta_{0}=300\arcsec$, and beam size of $\theta_{HPBW}=0.75\arcsec$ the estimated peak reduction amounts to about 2%. The distance was chosen as a minimal distance between two different pointing centers. This is illustrated in Fig. 12 where we show the peak over total flux density for $\mathrm{S/N}>200$ sources in different pointings. An offset of $\sim 2.5\%$ is present in this diagram, however, it is not distance dependent, and thus unlikely to be related to bandwidth smearing. Thus, for the reasons outlined above, we do not apply any corrections for the BWS effect.

3.6 The 3 GHz VLA-COSMOS Large Project catalog

A sample page of the catalog is shown in Table 5. For each source, we report its ID, 3 GHz name, RA and DEC position (weighted centroid) and error on the position, total flux density with relative error666 The flux errors reported do not depend on the number of pixels used for integration, but scale with the source brightness (see Hales et al. 2012, 2014). , 3 GHz $rms$ calculated at the position of the source, S/N, number of pixels used in flux density integration, flag for resolved sources, and flag for multicomponent sources. The peak surface brightness of resolved sources can be obtained by multiplying the S/N with the $rms$ value. The catalog is available in electronic format in the COSMOS IRSA archive777http://irsa.ipac.caltech.edu/frontpage/.

4 Radio spectral indices

Given the wide bandwidth of our VLA-COSMOS 3 GHz survey and the existence of previous COSMOS radio surveys, we approached radio spectral index calculations in two ways. The first method uses the MSMF algorithm to construct spectral indices directly from our observed data by fitting a two-term Taylor polynomial to amplitudes between 2 and 4 GHz (MSMF-based spectral index or $\alpha_{\scriptsize\mbox{MSMF}}$ hereafter). The second method uses the cataloged monochromatic flux densities at 3 GHz in combination with the values taken from the 1.4 GHz joint catalog (Schinnerer et al., 2010) to calculate spectral indices between these two frequencies (1.4–3 GHz spectral index or $\alpha_{\scriptsize 1.4-3\leavevmode\nobreak\ \mbox{GHz}}$ hereafter). In Sect. 4.1 we investigate systematics in the MSMF spectral index maps, and compare the differently derived spectral indices. In Sect. 4.2 we derive the 1.4–3 GHz spectral index distribution for the full sample of the 3 GHz sources.

4.1 MSMF-based versus 1.4–3 GHz spectral indices

We can calculate the MSMF-based spectral indices defined for each source using the wide bandwidth of our observations if the source has a sufficient S/N between 2 and 4 GHz. These spectral indices should be viable for point sources that have S/N$>$10, and for diffuse emission that has S/N$>$100. To do so, a mosaic of spectral indices ($\alpha$-map) was generated by dividing the Taylor term 1 (TT1) mosaic by the Taylor term 0 (TT0) mosaic (see Rau & Cornwell, 2011). For each source, its spectral index was extracted from the pixel in the $\alpha$ map that corresponds to the pixel containing the peak surface brightness in the total intensity mosaic.

To investigate possible systematics in the $\alpha$ map due to wideband primary beam corrections we utilized the 106 bright, point-like sources introduced in Sect. 3.5. We derived MSMF-based spectral indices both in the mosaic and individual pointings for these sources. In the bottom panel of Fig. 11 we show the difference between such derived spectral indices as a function of distance from the pointing center. The MSMF spectral indices show a systematic steepening with increasing distance, which likely arises due to an imperfect primary beam correction of TT1.888The MSMF algorithm is still in active development and the upcoming software versions should correct for this. To correct for this effect a posteriori (as necessary here), we performed a linear fit to the trend. We then applied this distance-dependent correction to each $\alpha$-map pointing prior to mosaicking to generate an $\alpha$ mosaic corrected for this effect.

In Fig. 13 we compare the corrected spectral indices from MSMF with those derived from the cataloged flux densities at 3 GHz and 1.4 GHz (joint catalog, Schinnerer et al. 2010). The catalogs were cross-matched using a search radius of 1.3$\arcsec$, which is half of the beam size of the lower resolution (1.4 GHz) survey. The sample was further limited to single-component sources with S/N$>5$ in the 1.4 GHz catalog yielding a total of 2,191 sources. Although there are no systematic offsets within the error margins, there is a rather large scatter between the spectral indices obtained with these two methods. A non-negligible portion of this spread is due to the large uncertainty on the in-band (i.e., MSMF derived) spectral indices; a point-like source with S/N $\sim 50$ and $\alpha=0.7$ has an uncertainty of $\sim 0.1$ in its in-band spectral index (see Condon, 2015). If the MSMF spectral indices had not been corrected, there would have been a systematic offset of -0.2 across the entire S/N range.

In summary, the MSMF-based spectral indices require further corrections after PB corrections are applied to the data. As a result of this, and the large scatter observed between the MSMF-based and 1.4–3 GHz spectral indices, we do not include the MSMF-based spectral indices in the final catalog. New CASA software versions should intrinsically correct for this. For the further analysis of spectral indices presented here we have, therefore, used only the values based on flux density measurements at 3 and 1.4 GHz.

4.2 1.4–3 GHz spectral indices

A high percentage of 3 GHz sources do not have a counterpart in the 1.4 GHz survey because of the better sensitivity of our 3 GHz survey. We employed the survival analysis to properly constrain the distribution of spectral indices for our 3 GHz selected sample without introducing any bias due to neglecting sources not detected in one of the surveys. This is a statistical method that takes into account both direct detections as well as upper (or lower) limits (see Feigelson & Nelson 1985 and Schmitt 1985 for details).

We first cross-correlated and combined our 3 GHz catalog with the 1.4 GHz joint catalog (Schinnerer et al., 2010) with a maximum separation of 1.3$\arcsec$, but also including sources without counterparts in one survey or the other. We then removed all sources that fell outside of the area observed at 1.4 GHz as the area observed at 3 GHz is larger (2.6 square degrees.). This was performed to ensure the same area for both surveys. We also removed 1.4 GHz multicomponent sources (80) and their 3 GHz counterparts. The final sample contains 10,523 entries out of which 23% were detected in both surveys, 74% were detected only at 3 GHz, and 3% were detected only at 1.4 GHz, as illustrated in Fig. 14 (top panel). If a source was not cataloged in one of the surveys we used five times the local $rms$ value at the coordinates of the source as an upper limit on the flux density. Each nondetection at 1.4 GHz yielded one lower limit on spectral index, and similarly, each nondetection at 3 GHz yielded one upper limit.

A Gaussian fit to the distribution of spectral indices detected in both surveys (green line in Fig. 14, top panel) results in the peak at $\alpha=-0.84$ and a standard deviation of $\sigma=0.35$. As this result is valid only for the subsample of 3 GHz sources also detected at 1.4 GHz, we employed the survival analysis to account for the full 3 GHz detected sample. We therefore ran the survival analysis on a single-censored data set that only included detections in both surveys and lower limits. The method assumes that limits follow the same distribution as direct detections and generates a cumulative distribution for all sources in the sample. This is shown in the middle panel of Fig. 14. There was enough overlap between direct detections and lower limits enabling the survival analysis to properly constrain the median of the total distribution to $\alpha=-0.68\pm 0.02$, even though there were three times more limits than detections. This method however cannot constrain all lower limits and the cumulative function does not converge to 0, yielding a total of 6% unconstrained sources. To constrain these (as needed to derive the probability density function, PDF, for spectral indices; see below) we employed a physical argument that a radio source exhibiting standard synchrotron self-absorption cannot have a spectral index higher than $\alpha_{\mbox{\scriptsize max}}=2.5$ (Rees 1967; unless it is extremely rare and exotic; for example see Krishna et al. 2014). Our data can also constrain the distribution of spectral indices only up to value of $\alpha=0.8$, since this interval contains 99.5% of sources detected at both 1.4 and 3 GHz. With these limits we can at best assume a flat probability that unconstrained sources have $0.8<\alpha<2.5$, and we can formally extrapolate the cumulative distribution function to 0 (red dashed line in Fig. 14, middle panel). Having constrained this, we then derived the PDF for spectral indices of our 3 GHz sources by calculating the first derivative of the cumulative distribution function extrapolated to 0. The PDF is shown in the bottom panel of Fig. 14. The best-fit Gaussian to the PDF yields a mean of $\alpha=-0.73$ and a standard deviation of $\sigma=0.35$. Both the median of the distribution and the mean of the Gaussian fit agree very well with previous work carried out on spectral indices (e.g., Condon 1984; Lisenfeld & Völk 2000; Kimball & Ivezić 2008) and we conclude that our catalog flux densities do not show any significant systematics.

5 Radio source counts corrections

A well-established approach to estimate the combined effects of noise bias, source extraction and flux determination systematics, inhomogeneuos noise distribution over the imaged field, and resolution bias on the measured source counts (completeness and bias corrections hereafter) is to rely on mock samples of radio sources, as described in Sect. 5.1. As these corrections do not take into account the fraction of spurious sources (as the mock sources are always inserted into the same mosaic) in Sect. 5.2,  we separately derive the false detection rate. The combination of the two corrections then yields the net radio source count corrections.

5.1 Completeness and bias corrections

We here describe the Monte Carlo simulations used to derive the completeness and bias corrections. Mock sources were injected over the imaged field and then recovered using the same technique adopted for the real radio sources, as detailed in Sect. 5.1.1. The flux density and size distributions assumed are described in Sect. 5.1.2 and Sect. 5.1.3, respectively. The results of the final simulations yielding the adopted completeness and bias corrections are detailed in Sect. 5.1.4, and a summary of the effects taken into account by the completeness and bias corrections is given in Sect. 5.1.5.

5.1.1 Retrieval of mock radio sources

The procedure adopted to insert and retrieve mock sources in and from the mosaic is as follows. Since blobcat does not produce a residual map, we inserted mock sources (Gaussian in shape) directly into the continuum map avoiding already cataloged components. This procedure was limited to the central two square degrees of the mosaic. For each mock source, a square shape with a width of $6\sigma+21$ pixel on the side was required to be free of any cataloged emission (real or mock), where $\sigma$ is the standard deviation along the Gaussian major axis. The positions were randomly chosen until this was satisfied. At a resolution of 0.75$\arcsec$ the continuum map is mostly empty of sources and confusion is negligible. We did not observe any systematic clustering of mock sources toward the less populated parts of the mosaic (more noisy parts closer to the edge for example) by requiring no overlap between the components. After all mock sources were inserted, we ran blobcat with the exact parameters as performed for the real sources. Since the extraction was carried out on a map containing both real and mock emission, all the 10,899  real components were always recovered and then removed from the extracted catalog, prior to further processing. To generate realistic mock catalogs of radio sources, however, we needed to assume i) a flux density distribution in (and below) the range tested by the observations and ii) an angular size distribution of the radio sources. This is described in detail in the following sections.

5.1.2 Flux density distribution

We simulated the flux density distribution using both a simple power-law model (PL model) and a multinode power law (MPL model) that better reproduces the observed source counts below $500\leavevmode\nobreak\ \mu$Jy. In the former case we used the 1.4 GHz source counts from previous surveys scaled to 3 GHz (see Bondi et al. 2008). The multinode power-law model is that derived by Vernstrom et al. (2014) (see their Table 4, Zone 1). For both models the mock catalogs were generated down to a total flux density of 5 $\mu$Jy and contained more than 40,000 (65,000) objects in the PL (MPL) model. This also allowed us to count sources with flux densities below the S/N threshold as positive noise fluctuations might lead these to have a measured peak flux density above our source detection threshold. As shown below, the results of our simulations do not yield differences between the two models, and we adopted the MPL model for our final simulations.

5.1.3 Angular size distribution

We needed to assign an intrinsic angular size to each mock source. Unfortunately, a satisfying description of the intrinsic source angular size distribution at sub-mJy flux density is still missing and we needed to rely on extrapolations from higher flux densities. Bondi et al. (2008) used a simple power-law parametrization distribution of the angular sizes of the sources as a function of their total flux density. We followed the same method with some adaptations, as described below.

The angular size ($\theta$) distribution was simulated, assuming a power-law relation between angular size and flux density ($\theta\propto S^{n}$). This relation was normalized using the cumulative angular size distribution derived at $\sim 1$ mJy from the VVDS 1.4 GHz observations with a resolution of $6\arcsec$ (Bondi et al., 2003). The relatively low resolution of the VVDS survey allowed us to avoid bias against sources with angular sizes of up to $15\arcsec$ in our simulations (Bondi et al., 2003). We explored the range of $n$ values between 0.3 and 1.0 in steps of 0.1. To infer the best $n$ value, the angular size distribution of the sources from the catalog in a specific total flux density range was compared with the corresponding distribution derived from the mock samples with different $n$ values. The value of $n$ that gave the best match between the angular size distribution of observed and mock sources was then chosen as the best approximation for the intrinsic source size versus total flux density relation.

Since the observed source angular sizes are not provided by blobcat, these were estimated using the relation between the ratio of the total flux density and peak surface brightness and angular sizes,

[TABLE]

where $S_{t}$ is the total flux density, $S_{p}$ is the peak surface brightness, $\theta_{b}$ is the FWHM of the circular beam ($0.75\arcsec$ in our observations), $\theta_{M}$ and $\theta_{m}$ are the major and minor FWHM intrinsic (deconvolved) angular sizes; see Bondi et al. 2008, where the same approach was used to derive a size estimate of sources affected by bandwidth smearing. In doing so we needed to make some assumptions on the morphology of the sources and in particular on how the sources are, eventually, resolved. We considered two limiting cases as follows:

Elongated geometry: sources are resolved in only one direction. This implies that $\theta_{m}=0$ and

[TABLE]

The simulated mock sources were accordingly generated as sources extended in one direction and eq. 2 is the appropriate relation between $S_{t}/S_{p}$ and the angular size. 2. 2.

Circular geometry: sources are uniformly resolved in all directions. This implies that $\theta_{M}=\theta_{m}$ and

[TABLE]

The simulated mock sources were accordingly generated as sources uniformly extended in all directions and eq. 3 is the appropriate relation between $S_{t}/S_{p}$ and the angular size.

Mock catalogs were generated for each combination of the 2 source count models (PL and MPL), the 8 different $n$ values (0.3–1.0 in steps of 0.1), and the 2 different source geometries ($\theta_{m}=0$ or $\theta_{m}=\theta_{M}$). For each of these 32 combinations we generated and merged 10 different mock samples. Then, we derived for each of the 32 different mock catalogs the $S_{t}/S_{p}$ distributions for sources with $S_{t}<100$ $\mu$Jy, splitting them into two subranges: $S_{t}\leq 40$ $\mu$Jy, and $40<S_{t}\leq 100$ $\mu$Jy. This range in total flux density is more affected by the choice on the intrinsic source size distribution and therefore is the best suited for a comparison between the $S_{t}/S_{p}$ distribution of the real cataloged sources and that derived from the mock samples reprocessed as the observed catalog. Using $S_{t}/S_{p}$ as a proxy for the angular size of the radio sources has the advantage that we do not need to deal with upper limits in the measured source sizes because of sources classified as unresolved.

The results of this comparison can be summarized as follows. No significant differences were found using the PL or MPL distributions for the source counts. For this reason, we were able to adopt either of the two models in the following analysis and we decided to use the MPL model, which provides a more realistic and detailed description of the observed source counts. However, none of the 16 combinations of $n$ value and source geometry provided a satisfying match between the $S_{t}/S_{p}$ distribution of the reprocessed mock catalog and that of the observed catalog, in the flux density range $S_{t}<100$ $\mu$Jy. While some combinations of parameters provided a reasonable match for sources with $S_{t}\gtrsim 40$ $\mu$Jy, they all failed to reproduce the observed distribution of $S_{t}/S_{p}$ below this threshold. In particular, the mock samples showed lower values of $S_{t}/S_{p}$ than the catalog for $S_{t}\lesssim 40$ $\mu$Jy. This is shown in Fig. 15 where we plot in the two panels the $S_{t}/S_{p}$ distribution for sources with $S_{t}<40$ $\mu$Jy and sources with $40<S_{t}\leq 100$ $\mu$Jy, respectively. Together with the observed distribution derived from the sources in the catalog we also plot the distribution obtained from our original simulation using $n=0.6$ and elongated geometry. The two distributions are clearly shifted and this effect is found in all the simulations.

This result is not completely unexpected, The extrapolation to very low flux density of our power-law relation between angular size and flux density, which has been previously tested only for sources more than one order of magnitude brighter, produces mock samples of radio sources dominated by extremely compact objects at the faint end of the total flux density distribution. For instance, for the simulations shown in Fig. 15, 45% of all the sources with $S_{t}\leq 40$ $\mu$Jy have $S_{t}/S_{p}<1$. This result is at odds with the distribution of the observed catalog, where only 26% of the observed sources fainter than 40 $\mu$Jy have $S_{t}/S_{p}<1$.

The simplest way to decrease the number of extremely compact objects at the faint end of the flux distribution in our mock sources, without modifying the adopted power-law relation between angular size and flux density, is to apply a minimum angular size to the faint mock sources. We tested the following expression for $\theta_{\mathrm{min}}$:

[TABLE]

in which the exponential part is motivated by the fact that, on the basis of the analysis shown in the lower panel of Fig. 15, no $\theta_{\mathrm{min}}$ is required at flux densities $\gtrsim 40$ microJy. We included the minimum angular size requirement in our procedure to generate the mock samples of radio sources and repeated the simulations, extraction process, and comparison of the $S_{t}/S_{p}$ distributions. By varying the parameters $k_{1}$ and $k_{2}$ we found that the best value for the $k_{2}$ parameter is $k_{2}=40$ $\mu$Jy, while for the normalization factor $k_{1}$ is equal to 0.3 (for the elongated geometry) and 0.2 (for the circular geometry). The different normalization is necessary because for a given intrinsic source size the area covered by a circular source is larger and derived from eq. 2 and 3. As shown in Fig. 15, this time we found a very good agreement between the observed and simulated distributions of $S_{t}/S_{p}$ at low flux densities as well. In particular, for the simulation shown in Fig. 15, introducing a minimum angular size as a function of the total flux density reduces the percentage of faint sources ($S_{t}\leq 40$ $\mu$Jy) with $S_{t}/S_{p}<1$ from 45% to 30% close to the observed value of 26%. Thus, we adopted the above-described size parametrization for our final simulations used to derive the completeness and bias corrections.

We further found that the geometry of the radio sources has some effects on the results we obtained. We obtained the best match between the $S_{t}/S_{p}$ distributions for $n=0.5-0.6$ for elongated geometry ($\theta_{m}=0$), and we obtained the best match for $n=0.6-0.7$ for circular geometry ($\theta_{m}=\theta_{M}$). We note that both the assumptions made on the source geometry are clearly simplistic and real sources will consist of a mix of elongated and circular sources. Thus, to compute our final completeness and bias corrections for the adopted MPL flux density distribution models we computed the completeness and bias correction factors using the average of those from the four best-matched simulations: i) elongated sources with $n=0.5$, and $n=0.6$, and ii) circular sources with $n=0.6$, and $n=0.7$, as described in more detail in the next section.

5.1.4 Derivation of completeness and bias corrections

We generated 60 mock catalogs using the parameterization as described above (Sect. 5.1.2 and Sect. 5.1.3; see also below). The mock sources were inserted into the mosaic and retrieved as described in Sect. 5.1.1. The retrieved mock sources were then cross-correlated with the input mock catalog and their measured total flux density chosen to be either their integrated flux density if resolved, or peak surface brightness if unresolved. For this, the same $S_{t}/S_{p}$ envelope was used as described in Sect. 3. Lastly, successfully extracted mock sources and original mock sources were binned separately in flux densities. The ratio of their numbers in each flux density bin represents the completeness and bias correction factor.

In Fig. 16 we show the net result of the above-described Monte Carlo simulations for the MPL model and best-matched simulations, i.e., i) elongated sources with $n=0.5$ and $0.6$, and ii) circular sources with $n=0.6$ and $0.7$. We take the average of these simulations as the completeness and bias correction with a confidence interval that takes into account the differences within the six sets of simulations. This is tabulated in Table 2. For reference, in Fig. 16 we also show the average completeness and bias corrections obtained using only the elongated and circular geometry approximations. The net curve yields values of about 55% at $12\leavevmode\nobreak\ \mu$Jy (S/N=5.2), and rather constant up to $20\leavevmode\nobreak\ \mu$Jy (S/N $<$9), beyond which they rise to a 94% completeness above $40\leavevmode\nobreak\ \mu$Jy (S/N$\geq$16). The mean error of the completeness and bias corrections is 5%. The two (elongated and circular geometry) approximations are consistent up to $\sim 30\leavevmode\nobreak\ \mu$Jy, beyond which they start diverging with the circular approximation being systematically lower at higher flux densities. However, beyond this limit both curves saturate at fairly constant values ($\sim 0.92$ for circular and $0.98$ for elongated morphology), implying average values of over 95% for fluxes higher than $\sim 40$ $\mu$Jy.

5.1.5 Biases addressed

There are several effects and biases that occur in the cataloging process that we addressed through our simulations. Firstly, an incompleteness in the extracted catalog exists as real sources on the sky are not be detected if i) their peak surface brightness falls below the chosen threshold of $5\sigma$ because of fluctuations in the local $rms$, or ii) they are extended enough for their peak surface brightness to fall below the detection threshold, even though their integrated flux density is well above it. Secondly, a contamination effect is also present. If a source is detected, its flux density might be wrongly computed because of the presence of a noise peak. Statistically, this effect is mostly symmetric around the mean flux density. However, when we set the total flux density of an unresolved source to its peak surface brightness we may introduce an asymmetric bias toward smaller flux densities. Some sources with $S_{t}>S_{p}$ within the envelope in Fig. 9 might truly be resolved, however noise variations do not allow us to determine this with sufficiently high accuracy leading to a potential bias. The final result is that a source can jump into a flux density bin where it does not belong, thus increasing its contamination. The combination of completeness, which always decreases with decreasing flux density, and the significant number of sources that move from their original flux density bin to another owing to errors in flux measurement at faint flux densities, produces the flat distribution of the completeness and bias correction factor seen at flux densitites of $\sim 12-20\leavevmode\nobreak\ \mu$Jy in Fig. 16.

In summary, the simulations we performed account for both the fraction of nondetected sources (incompleteness), and also the redistribution of sources between various flux density bins. Thus, in principle, its value can be larger than one if the contamination is high. These corrections, however, do not take into account the fraction of spurious sources as a function of flux density, which are derived separately in the next section.

5.2 False detection rate

To assess the false detection rate of our source extraction we ran blobcat on the inverted (i.e., multiplied by $-1$) continuum map with the same settings used for the main catalog. Since there is no negative emission on the sky, every source detected in the inverted map is per definition a noise peak (i.e., a false detection). The source extraction returned 414 negative detections with S/N $\geq 5$ across the entire observed field, 95 of which were outside the central two square degrees, demonstrating that 23% of false detections lies at the edge of the mosaic.

The highest S/N negative detections were predominantly located around true bright sources as they suffer from artifacts; up to six negative components could be found around a single bright object due to the VLA synthesized beam shape (see also Sect. 7.1.1. in Vernstrom et al. 2014 for an explanation of this effect). Since extraction of real emission does not exhibit this behavior, we removed all negative components that were less than 3$\arcsec$ away from a real source with S/N$>100$. This step removed further 40 components. We additionally removed four sources with catastrophic peak estimates, which increased their S/N by more than a factor of four owing to poor parabola fits. We note that there were no such sources in the catalog of real emission. The remaining 275 negative detections within the inner two square degrees were then classified into resolved and unresolved using the same envelope as was carried out for the real data. Finally, they were binned in S/N and flux densities alongside true detections to enable direct comparison. The results are shown in Fig 17 and also listed in Table 3. As expected, only the lowest S/N bins have a noticeable portion of false detections ($24\%$ for S/N$=5.0$–$5.1$), which quickly decreases to less than 3% for any S/N bin at S/N$>5.5$. The estimated fraction of spurious sources over the entire catalog above S/N$>5\,(5.5)$ drawn from the inner two square degrees is only 2.7 (0.4)% .

6 Radio source counts

In this section we present our 3 GHz radio source counts (Sect. 6.1), and compare them to 3 GHz and 1.4 GHz counts available in the literature (Sects. 6.2 and 6.3, respectively).

6.1 VLA-COSMOS 3 GHz radio source counts

We present our 3 GHz radio source counts normalized to Euclidean geometry, both corrected and uncorrected, in the top panel of Fig. 18. In Table 4 we list the counts, errors, number of sources, and radio source count corrections (i.e., completeness and bias and false detection fraction corrections) in each flux density bin. The source count errors take into account both the Poissonian errors as well as completeness and bias correction uncertainties. Most of our sources are located at low flux densities (below 0.5 mJy) with more than 500 sources in each flux density bin below 60 $\mu$Jy resulting in small Poissonian errors. As evident from the plot, our source counts at 3 GHz exhibit a flattening at about 0.3 mJy that continues one order of magnitude in flux densities down to 30 $\mu$Jy, steepening further at fainter flux densities.

6.2 Comparison with 3 GHz counts from the literature

In the middle panel of Fig. 18 we compare our 3 GHz source counts with other 3 GHz counts available in the literature (Condon et al., 2012; Vernstrom et al., 2014). Condon et al. (2012) performed a $P(D)$ analysis using 3 GHz confusion-limited data based on 50 hours of on-source C-array configuration observations of one VLA pointing targeting the Lockman hole and reaching an $rms$ of $1\leavevmode\nobreak\ \mu$Jy beam*-1*. Fitting single power-law models to the data the analysis allowed these investigators to constrain the counts of discrete sources in the $1-10\leavevmode\nobreak\ \mu$Jy range, also shown in Fig. 18. Vernstrom et al. (2014) performed a more complex $P(D)$ analysis on the same data fitting various (modified power-law, and node-based) models to the data allowing them to probe the counts down to $0.1\leavevmode\nobreak\ \mu$Jy. In Fig. 18 we show the counts based on the fit of a phenomenological parametric model of multiple joined power laws (their node-based model) applied to the inner circular area with a $5\arcmin$ radius (their Zone 1; see Vernstrom et al. 2014 for details).

The counts derived here are in very good agreement with those derived by Condon et al. (2012). Fitting the five faintest flux density bins using a power law, $dN/dS\propto S^{\gamma}$, we find that the slope $\gamma=-1.72$ is perfectly consistent with that inferred by Condon et al. (2012), while our normalization is slightly lower. Our comparison to the Vernstrom et al. (2014) results shows that the counts are in agreement down to $\sim 30\leavevmode\nobreak\ \mu$Jy with a discrepancy at fainter flux densitites as our counts are systematically lower than theirs.

In general, the strength of the $P(D)$ analysis is the ability to probe counts below the nominal noise in the data, while avoiding resolution biases as it is applied on confusion-limited (thus, low resolution) data. However, as the $P(D)$ analyses discussed above were performed on a single VLA pointing, the resulting counts may be subject to cosmic variance due to the small area covered. This could potentially explain the observed discrepancy between the VLA-COSMOS 3 GHz Large Project counts based on a two square degree area and the Vernstrom et al. (2014) results based on a 0.022 square degree area (their Zone 1). To test this we subdivided the two square degree COSMOS field into 100 square and nonoverlapping subfields, each with an area of 0.020 deg2 roughly corresponding to a circle with a radius of $5\arcmin$. In the middle panel of Fig. 18, we show the range of such obtained counts (corrected for completeness and bias, and false detection fractions, calculated on the full two square degrees and described in Sect. 5). We find that sample variance that is quantified in this way can introduce a (1$\sigma$) scatter of ${}^{+0.1}_{-0.2}$ dex in the source counts. The distribution in counts in the 100 subfields are likely to be an underestimate of the true cosmic variance, which is dominated by cosmic large-scale structure, rather than sample variance, because these fields are likely not fully independent from each other. Thus, cosmic variance may explain the observed discrepancy.

6.3 Comparison with 1.4 GHz counts from the literature

To compare our result with more abundant 1.4 GHz observations and models (e.g., Condon 1984; Bondi et al. 2008; Owen & Morrison 2008; Wilman et al. 2008; de Zotti et al. 2010; Condon et al. 2012) we scale our flux densities to the 1.4 GHz observed frame using a spectral index of -0.7. This value, which is also in agreement with the spectral index survival analysis described in Sect. 4, is commonly used and provides the easiest comparison (e.g., Condon et al. 2012). We show the 1.4 GHz source count comparison in the bottom panel of Fig. 18.

The large spread of the 1.4 GHz source counts available in the literature at submillijansky levels (see, e.g., Fig. 1 in Smolčić et al. 2015b) is usually attributed to a combination of i) cosmic variance as often the observed fields are rather small (see Fig. 1 and middle panel of Fig. 18), and ii) resolution bias leading to a loss of sources in radio continuum surveys conducted at intermediate to high ($\lesssim 2\arcsec$) angular resolution (as described in more detail in Sect. 5.1.5). The large, two square degree area of the COSMOS field minimizes the effect of cosmic variance, and in Sect. 5 we performed extensive Monte Carlo simulations to account for potential resolution biases. Our source counts agree well with those derived by Condon et al. (2012) based on the $P(D)$ analysis at the faint end (see previous section).

The counts derived here are in good agreement with those derived from the VLA-COSMOS 1.4 GHz Large Project (Schinnerer et al. 2007; Bondi et al. 2008; red diamonds in Fig. 18) at flux densities higher than $\sim 200\leavevmode\nobreak\ \mu$Jy, but are slightly lower in the flux density range of $100-200\leavevmode\nobreak\ \mu$Jy. As this is the same field, cosmic variance cannot explain the discrepancy. The uncorrected counts from the two surveys are in very good agreement; the difference is the largest in the flux density range where the 1.4 GHz survey is the least complete (about 60%), and the corrections, thus, are the largest. In the same flux density range the corrections for the 3 GHz survey are not as severe given the much higher sensitivity of the 3 GHz survey. Further reasons that could explain part of the discrepancy are i) the effect of BWS on the radio source count corrections that are present in the 1.4 GHz data, but are not present in the 3 GHz data (see Bondi et al. 2008), and ii) a possibly overly simplistic scaling of the 3 GHz counts to 1.4 GHz using just one spectral index value. Source counts at 1.4 GHz depend on the steepness of the counts at 3 GHz and the spread of the spectral indices. We leave the analysis of the potential bias in source counts due to this effect to an upcoming paper (Novak et al., in prep.).

The largest discrepancy between the counts derived here and those in the literature is observed relative to the Owen & Morrison (2008) results. Owen & Morrison (2008) have observed the Lockman hole at 1.4 GHz in A-, B-, C-, and D-array configurations with the VLA reaching an angular resolution of $\sim 1.6\arcsec$ and $rms\approx 2.7\leavevmode\nobreak\ \mu$Jy beam*-1*. To correct for the resolution bias, they assumed a source size distribution with an extended tail at the high end (see their Fig. 8) and that distribution remains constant as a function of flux density. The source count corrections are significant under these assumptions and result in a flat source count distribution at flux densities fainter than $\sim 200\leavevmode\nobreak\ \mu$Jy (green points in the bottom panel of Fig. 18). As already discussed by Condon et al. (2012) and Vernstrom et al. (2014), these corrections are most likely overestimated. In contrast, for the corrections applied to the data presented here we assumed a model for radio source sizes such that the radio size is a function of flux density with a limiting minimum size (see Sect. 5.1). The agreement between our source counts and those derived from confusion-limited data (Condon et al., 2012; Vernstrom et al., 2014) further strengthens the validity of this assumption.

In Fig. 18 we also compare our results with the models developed by Condon (1984), Wilman et al. (2008), and Béthermin et al. (2012). The faint end of our counts ($\lesssim 80\leavevmode\nobreak\ \mu$Jy), combined with the results from Condon et al. (2012) that appear as an extrapolation of our data, agree the best with the Condon (1984) model. The model was constrained by source counts, redshift, and spectral-index distributions for various 400 MHz to 5 GHz flux-limited samples as well as the local 1.4 GHz luminosity function for two dominant, spiral and elliptical galaxy populations. The adopted model is not a unique solution and evolves all sources, i.e., ellipticals and spirals, steep- and flat-spectrum sources, in the same way. At flux densities above $\sim 80\leavevmode\nobreak\ \mu$Jy the Condon (1984) model is slightly higher than our derived source counts and is consistent with the counts determined by Vernstrom et al. (2014).

Our derived source counts deviate from those predicted by the Wilman et al. (2008) and Béthermin et al. (2012) models. While they agree with the first down to $\sim 100\leavevmode\nobreak\ \mu$Jy, they are systematically higher at fainter flux densities. On the other hand, the Béthermin et al. (2012) model underpredicts our counts in the flux density range of $\sim 50-300\leavevmode\nobreak\ \mu$Jy, while it overpredicts the counts at flux densities $\lesssim 30\leavevmode\nobreak\ \mu$Jy. The discrepancies may possibly be understood when considering how AGN and star-forming galaxies were implemented in the models. Béthermin et al. (2012) implement only models for X-ray selected AGN ($L_{\mathrm{2-10\leavevmode\nobreak\ keV}}\sim 10^{42}-10^{44}$ erg s*-1*; see Mullaney et al. 2011, 2012; Aird et al. 2012), and thus ignore the population of radio-loud AGN hosted by red, quiescent galaxies, regularly not identified as X-ray AGN, yet still substantial (e.g., Best et al. 2006; Smolčić et al. 2008b; Bonzini et al. 2013; Smolčić et al, submitted; Delvecchio et al., accepted.). This could explain the lack of sources with flux densities in the range of $\sim 50-300\leavevmode\nobreak\ \mu$Jy in the model (see, e.g., Smolčić et al. 2008b; Padovani et al. 2015; Smolčić et al., submitted), compared to the observational results. On the other hand, Béthermin et al. (2012) model the star-forming galaxy population using the most recent results from Bouwens et al. 2007; Rodighiero et al. 2011; Magnelli et al. 2011; Karim et al. 2011; Sargent et al. 2012 by tracing the star-forming galaxy main sequence and the stellar mass function over cosmic time, while also taking into account main-sequence and starburst galaxy spectral energy distribution libraries. Therefore, the excess of Béthermin et al. (2012) model compared to that of Wilman et al. (2008) could suggest that the Wilman model carries potential for improvement in modeling the star-forming galaxy population.

7 Summary and conclusions

We presented the VLA-COSMOS 3 GHz Large Project based on 384 hours of observations with the Karl G. Jansky Very Large Array at 3 GHz (10 cm) toward the two square degree COSMOS field. Our final mosaic, imaged per pointing with the multiscale multifrequency algorithm and self-calibration, reaches a median $rms$ of 2.3$\leavevmode\nobreak\ \mu$Jy beam*-1* over the two square degrees, at an angular resolution of $0.75\arcsec$. We further presented a catalog of 10,830 radio sources. Combining our data with the 1.4 GHz VLA-COSMOS Joint Project data using survival analysis, we found the expected median spectral index $\alpha$ of -0.7. Comparing the positions of our 3 GHz sources with those from the high-resolution VLBA imaging at 1.4 GHz, we estimated that the astrometry is accurate to $0.01\arcsec$ at the bright end. Radio source count corrections were calculated for the central two square degrees and used to infer radio source counts. The radio angular size parametrization adopted based on the comparison of mock versus real source total over peak flux density ratios suggests that the angular sizes of radio sources at these flux density levels can be modeled as a power law in flux density ($\theta\propto S^{n}$) with a minimal, flux-dependent size cutoff (eq. 4). Our corrected radio counts with direct detections down to $20\leavevmode\nobreak\ \mu$Jy (at 1.4 GHz) are consistent with those derived based on P(D) analyses (Condon et al., 2012), and agree best with the Condon (1984) model, while they are systematically higher than those predicted by the SKADS (Square Kilometer Array Design Studies) simulations (Wilman et al., 2008).

The VLA-COSMOS 3 GHz Large Project simultaneously provides the largest and deepest radio continuum survey to date, bridging the gap between radio continuum surveys conducted with past generation and those planned with the next generation facilities. These radio data, in conjunction with the vast panchromatic COSMOS data sets, will allow for the exploration of various cosmologically relevant topics, such as i) the characteristics of the microJansky radio population, ii) radio-quiet QSOs by direct detection in the radio band, iii) modes of star formation at early cosmic epochs, and iv) studying stellar mass growth in typical galaxies since early cosmic epochs and star formation quenching via AGN feedback.

Acknowledgements.

We thank the referee, Jim Condon, for insightful comments that helped improve the manuscript. Based on observations with the National Radio Astronomy Observatory which is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This research was funded by the European Union’s Seventh Frame-work programs under grant agreements 333654 (CIG, ’AGN feedback’) and 337595 (ERC Starting Grant, ’CoSMass’). MB and PC acknowlege support from the PRIN-INAF 2014. AK and FB acknowledge support by the Collaborative Research Council 956, sub-project A1, funded by the Deutsche Forschungsgemeinschaft (DFG). BM and FB acknowledge support through DFG priority program 1573 founded by the DFG.