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

V. Smolčić¹, M. Novak¹, M. Bondi², P. Ciliegi³, K. P. Mooley⁴, E. Schinnerer⁵, G. Zamorani³, F. Navarrete⁶,
S. Bourke⁴, A. Karim⁶, E. Vardoulaki⁶, S. Leslie⁵, J. Delhaize¹, C. L. Carilli⁷, S. T. Myers⁷, N. Baran¹, I. Delvecchio¹,
O. Miettinen¹, J. Banfield⁸, ⁹, M. Baloković⁴, F. Bertoldi⁶, P. Capak¹⁰, D. A. Frail⁷, G. Hallinan⁴, H. Hao¹¹,
N. Herrera Ruiz¹², A. Horesh¹³, O. Ilbert¹⁴, H. Intema⁷, V. Jelić¹⁵, ¹⁶, ¹⁷, H.-R. Klöckner¹⁸, ¹⁹, J. Krpan¹, S. R. Kulkarni⁴,
H. McCracken²⁰, C. Laigle²⁰, E. Middleberg¹², E. J. Murphy²¹, M. Sargent²², N. Z. Scoville⁴, and K. Sheth²¹

¹ Department of Physics, University of Zagreb, Bijenička cesta 32, 10002 Zagreb, Croatia
e-mail: vs@phy.hr

² Istituto di Radioastronomia di Bologna – INAF, via P. Gobetti, 101, 40129 Bologna, Italy

³ INAF–Osservatorio Astronomico di Bologna, via Ranzani 1, 40127 Bologna, Italy

⁴ California Institute of Technology, MC 249-17, 1200 East California Boulevard, Pasadena, CA 91125, USA

⁵ Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany

⁶ Argelander Institut für Astronomie, Auf dem Hügel 71, 53121 Bonn, Germany

⁷ National Radio Astronomy Observatory, PO Box 0, Socorro, NM 87801, USA

⁸ CSIRO Australia Telescope National Facility, PO Box 76, Epping, NSW 1710, Australia

⁹ Research School of Astronomy and Astrophysics, Australian National University, Weston Creek, ACT 2611, Australia

¹⁰ Spitzer Science Center, 314-6 Caltech, Pasadena, CA 91125, USA

¹¹ Smithsonian Astrophysical Observatory, 60 Garden St., Cambridge, MA 02138, USA

¹² Astronomisches Institut, Ruhr-Universität Bochum, Universitätsstr. 150, 44801 Bochum, Germany

¹³ Benoziyo Center for Astrophysics, Weizmann Institute of Science, 76100 Rehovot, Israel

¹⁴ Aix Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille), UMR 7326, 13388 Marseille, France

¹⁵ Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV, Groningen, The Netherlands

¹⁶ ASTRON – The Netherlands Institute for Radio Astronomy, PO Box 2, 7990 AA, Dwingeloo, The Netherlands

¹⁷ Ruder Bošković Institute, Bijenička cesta 54, 10000 Zagreb, Croatia

¹⁸ Subdepartment of Astrophysics, University of Oxford, Denys-Wilkinson Building, Keble Road, Oxford OX1 3RH, UK

¹⁹ Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany

²⁰ Sorbonne Université, UPMC Univ. Paris 06, and CNRS, UMR 7095, IAP, 98bis Bd Arago, 75014 Paris, France

²¹ National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA

²² Astronomy Centre, Department of Physics and Astronomy, University of Sussex, Brighton, BN1 9QH, UK

Received 13 April 2016 / Accepted 20 June 2016

## ABSTRACT

We present the VLA-COSMOS 3 GHz Large Project based on 384 h 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 μJy 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σ, 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*<sub>3</sub> GHz > 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 α = −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.

**Key words.** catalogs – cosmology: observations – radio continuum: galaxies

## 1. Introduction

One of the main goals 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

<footnote>* The catalog is available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/602/A1</footnote>
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, 2015; Smolčić
2009; Smolčić & Riechers 2011).

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, 2003, 2012; Ciliegi et al. 1999; Georgakakis et al. 1999; Bock et al. 1999; Prandoni et al. 2001; 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. 2008, 2014; Owen & Morrison 2008; Miller et al. 2008, 2013; Owen et al. 2009; Hales et al. 2014). These past surveys have shown that deep observations at high angular resolution (≤1′′) 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. 2008, 2009b,a; 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, Australian Square Kilometre Array Pathfinder, MeerKAT, and Square Kilometre Array; 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 (Tisanić 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-h on-source with the VLA in C-array configuration. The ob-
servations resulted in a confusion-limited map with an rms of
1 µJy beam⁻¹. Based on this they constrained the counts of
discrete sources in the 1–10 µJy range via a *P*(D) analysis. A

**2. Observations and data reduction**

**2.1. Observations**

A total of 384 h 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 observa-
tions were taken from November 2012 to January 2013, June to
August 2013, and February to May 2014 in A-array (324 h) and

<figure><img src="image_1.png" /><figcaption>Fig. 1. Sensitivity (at the observed frequency of the given survey) vs. area for past, current, and future radio continuum surveys.</figcaption></figure>

more complex *P*(D) analysis using the same data was applied by
Vernstrom et al. (2014) who probed the counts down to 0.1 µJy.
Both results are qualitatively in agreement with the already well-
known flattening of the radio source counts (normalized to the
*N*(S) ∝ S⁻³/² of a static Euclidian space) below flux densities
of S₁.₄ GHz ≈ 1 mJy, and a further decrease of the counts with
decreasing flux density below S₁.₄ GHz ≈ 60 µ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 exper-
iment (Fixsen et al. 2011) 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 3 GHz
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, calibra-
tion, 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$.
<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Pointing pattern used for the 3 GHz VLA-COSMOS Large Project. The centers of the 192 pointings are indicated by the plus signs. Circles indicate the primary beam of each pointing, represented here by the HPBW at 3 GHz (15'; the primary beam HPBW is a function of frequency and varies by a factor 2 between the lower and upper edge of the S band).</figcaption></figure>

C-array configurations (60 h; Legacy ID AS1163). Sixty-four pointings, separated by 10' 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' in RA and Dec, while the third set is shifted by -5' in RA. Observing runs of 5 and 3 h in length were conducted. In each observing run J1331+3030 was observed for flux and bandpass calibration for about 3-5 min 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 min for 1m 40s on-source for gain and phase calibration, while the source J0713+4349 was observed for 5 min 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-h observing blocks was changed between the different observing runs to optimize the uv coverage. During the 3-h 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 Cookbook<fnref n="1" />.

In brief, the data are first loaded with the Obit<fnref n="2" /> 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 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 ~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

<footnote n="1">¹ http://www.aips.nrao.edu/cook.html</footnote>

<footnote n="2">² http://www.cv.nrao.edu/~bcotton/Obit.html</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Raw spectra of the gain calibrator source, i.e., phase vs. channel (top frame in each of the four panels) and amplitude vs. channel (bottom frame in each of the four panels) for the right-right (RR) and left-left (LL) polarizations (top panels and bottom panels, respectively). The panels to the left are for one night of observation in the A-array configuration and the panels to the right are for a C-array observation. No calibration was applied. All baselines and all pointings of the gain calibrator source were combined to produce these plots. We note the RFI in sub-bands 2 and 3.</figcaption></figure>

<table><caption>**Table 1.** Amplitude of the phase calibrator (J1024-0052) in each SPW for all observing blocks.</caption><thead><tr><th>SPW</th><th>Frequency (GHz)</th><th>Mean flux density (Jy)</th><th>Median flux density (Jy)</th><th>Standard deviation (Jy)</th></tr></thead><tbody><tr><td>1</td><td>2.060</td><td>0.739</td><td>0.735</td><td>0.029</td></tr><tr><td>4</td><td>2.444</td><td>0.707</td><td>0.704</td><td>0.025</td></tr><tr><td>5</td><td>2.572</td><td>0.700</td><td>0.696</td><td>0.023</td></tr><tr><td>6</td><td>2.700</td><td>0.684</td><td>0.680</td><td>0.023</td></tr><tr><td>7</td><td>2.828</td><td>0.668</td><td>0.665</td><td>0.024</td></tr><tr><td>8</td><td>2.956</td><td>0.652</td><td>0.648</td><td>0.026</td></tr><tr><td>9</td><td>3.084</td><td>0.645</td><td>0.642</td><td>0.020</td></tr><tr><td>10</td><td>3.212</td><td>0.635</td><td>0.632</td><td>0.020</td></tr><tr><td>11</td><td>3.340</td><td>0.625</td><td>0.622</td><td>0.021</td></tr><tr><td>12</td><td>3.468</td><td>0.615</td><td>0.611</td><td>0.020</td></tr><tr><td>13</td><td>3.596</td><td>0.603</td><td>0.600</td><td>0.021</td></tr><tr><td>14</td><td>3.724</td><td>0.539</td><td>0.579</td><td>0.129</td></tr><tr><td>15</td><td>3.852</td><td>0.525</td><td>0.569</td><td>0.153</td></tr><tr><td>16</td><td>3.980</td><td>0.535</td><td>0.566</td><td>0.124</td></tr></tbody></table>

Applications (CASA<fnref n="3" />)-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 *FLAGDATA*<fnref n="4" />, 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⁻¹) 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 side-lobes 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

<footnote n="3">³ CASA 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)</footnote>

<footnote n="4">⁴ In total, about 30–35% of the data were flagged (using the tasks *RFLAG* and *FLAGDATA*).</footnote>
<figure><img src="image_4.png" /></figure>

**Fig. 4.** Artifacts around bright sources before (left panels) and after (right panels) applying self-calibration phase-gain solutions. The right column panels also have tapering applied that circularizes the beam shape. The top three rows show the same source with a peak surface brightness of around $S_p \approx 16$ mJy beam$^{-1}$, but located inside three different pointings that were observed during different time epochs. The fourth row shows artifacts around the brightest source in our data ($S_p \approx 18$ mJy beam$^{-1}$), which is also extended. The final row illustrates the improvement when applying self-calibration solutions only from the phase calibrator as this source with $S_p \approx 2$ mJy beam$^{-1}$ has insufficient S/N for self-calibration (see text for details).

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×, 5×, and 10× 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 8000 pixels on-the-side with a pixel size of 0.2 × 0.2 arcsec$^2$. 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 mosaic$^5$. Synthesized beam size variations between different pointings were about 0.03", which was small enough to allow restoration of every cleaned pointing to an average circular beam of 0.75". 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') using

<figure><img src="image_5.png" /></figure>

**Fig. 5.** Left: final (A- plus C-array configuration combined) dirty beam of one pointing, after tapering. This beam was used in the cleaning of that pointing (see text for details). Right: mean stack of all 192 dirty beams. The contribution of radial sidelobes is 10% at maximum.

the WIDEBANDPBCOR task. The noise level in the phase center of an individual pointing was usually around 4–5 µ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", with a median rms of 2.3 µ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_f$) 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

<footnote n="5">⁵ A 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" 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.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Final VLA-COSMOS 3 GHz MSMF mosaic with overlaid Gaussian fits to the pixel surface brightness distributions in various mosaic sectors. The rms obtained via the Gaussian fit (in units of µJy beam⁻¹) is indicated in each panel. The panels shown cover the full COSMOS two square degree field. The small-scale (~1′) rms variations due to the pointing layout are less than 2%.</figcaption></figure>

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 densities were computed by *blobcat* as described in the previous section. For a perfect Gaussian unresolved source, the peak surface brightness in Jy beam⁻¹ 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<sub>p</sub>/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<sub>t</sub> < S<sub>p</sub> and mirrored it above the S<sub>t</sub>/S<sub>p</sub> = 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<sub>t</sub>/S<sub>p</sub> = 1 + 6 × (S/N)⁻¹.⁴⁴. We consider the 3975 components above the upper envelope as resolved. These resolved components were flagged in the catalog. For the unresolved
<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Stamps from the VLA-COSMOS 3 GHz continuum mosaic imaged with the MSMF algorithm showing examples of extended and compact radio sources.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Visibility plot showing the total area covered down to a given noise level (black dashed line). Our data extend beyond the COSMOS two square degree field, which ensures more uniform noise inside it (red full line). The median noise level inside the COSMOS two square degrees is $\sigma = 2.3 \ \mu\text{Jy}\ \text{beam}^{-1}$.</figcaption></figure>

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

<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Ratio of total flux density to peak surface brightness as a function of S/N. Components below the upper envelope (gray points) indicated by red solid lines are considered unresolved, while those above the upper envelope (black points) are considered resolved (see text for details).</figcaption></figure>

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 2500 components. The inspected sample was a combination of the i) brightest 2500 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_{EST} > 1 + 30 \times (S/N)^{-1}$ (351 components). The $R_{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 identified 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.).
<figure><img src="image_10.png" /><figcaption>**Fig. 10.** Astrometry comparison between 3 GHz and 1.4 GHz VLBA data for 443 VLBA sources (PI: Middelberg; N. Herrera Ruiz et al., in prep.).</figcaption></figure>

### 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" in ΔRA and 0.00" in ΔDec and a standard deviation of 0.01" for both. We note that we did not correct the catalog entries for the 0.01" offset in ΔRA. We took the standard deviation value (0.01") 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 (~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 \le S_t/S_p \le 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_p$) relative to the peak surface brightness retrieved from the mosaic ($S_M$) as a function of distance from the pointing center. If our data were affected by BWS, $S_p/S_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) Eq. (12) for the reduction of peak response $I/I_0 \approx 1/\sqrt{1+0.46\beta^2}$, where $\beta = (\Delta v/v_0) \times (\theta_0/\theta_{HPBW})$ equals fractional bandwidth times offset in synthesized beam-widths. Using

<figure><img src="image_11.png" /><figcaption>**Fig. 11.** Comparison of peak surface brightnesses (top) and MSMF-based spectral indices (bottom) determined inside the mosaic ($S_M$, $\alpha_M$) and individual pointings ($S_p$, $\alpha_p$) as a function of distance from the pointing center for 106 bright, point-like sources ($0.9 \le S_t/S_p \le 1.1$, $S/N > 200$) observed in up to 11 neighboring pointings at varying distances from the pointing center (gray points). In both panels the large red points and their corresponding errors indicate median values and interquartile ranges inside 5 equally spaced distance bins. In the bottom panel a linear fit is performed on the median values to obtain the needed correction of the systematic trend across all distances (see text for details).</figcaption></figure>

the VLA channel width $\Delta v = 2$ MHz, central frequency $v_0 = 3$ GHz, distance from the phase center $\theta_0 = 300''$, and beam size of $\theta_{HPBW} = 0.75''$ 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 $S/N > 200$ sources in different pointings. An offset of ~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 error<fnref n="6" />, 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 archive<fnref n="7" /> and at the CDS.

## 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_{MSMF}$ hereafter). The second method uses the cataloged monochromatic flux densities at 3 GHz in combination

<footnote n="6">⁶ 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).</footnote>

<footnote n="7">⁷ http://irsa.ipac.caltech.edu/frontpage/</footnote>
<figure><img src="image_12.png" /><figcaption>**Fig. 12.** Comparison of peak surface brightnesses over total flux densities for 106 bright, point-like sources ($0.9 \le S_t/S_p \le 1.1$, $S/N > 200$) observed in up to 11 neighboring pointings at varying distances from the pointing center (gray points). The large red points and their corresponding errors indicate median values and interquartile ranges inside 5 equally spaced distance bins. The theoretical prediction of the bandwidth smearing effect is shown by the blue curve (see text for details).</figcaption></figure>

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_{1.4-3\,\text{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 wide-band 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$^8$. 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",

<figure><img src="image_13.png" /><figcaption>**Fig. 13.** Comparison between MSMF-based and 1.4–3 GHz derived spectral indices, where the first were corrected for the observed systematic trend illustrated in Fig. 11. Red symbols and the corresponding errors denote median spectral indices and interquartile ranges, respectively, for sources in different S/N bins ($10 < S/N < 500$). The black dashed line indicates the median value of the red circles set at 0.02.</figcaption></figure>

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 2191 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 ~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", 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 10523 entries

<footnote>$^8$ The MSMF algorithm is still in active development and the upcoming software versions should correct for this.</footnote>
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_{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.

<figure><img src="image_14.png" /></figure>

**Fig. 14.** Top panel: distribution of 1.4–3 GHz spectral indices for sources detected at both frequencies (green line), and only at 3 GHz (lower limits, blue line) or 1.4 GHz (upper limits, red line). Middle panel: cumulative distribution function (CDF; black line) and error estimate (gray shaded area) of spectral indices calculated using the survival analysis also taking lower limits into account. The red dashed line shows a linear extrapolation of the distribution to zero assuming the maximum theoretically attainable spectral index of $\alpha = 2.5$ (see text for details). Bottom panel: probability density function (PDF) for spectral indices calculated as a first derivative of the CDF extrapolated to 0 at high end (middle panel). A Gaussian fit to the distribution is also shown (red curve) and its mean and standard deviation are indicated in the panel.

### 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 Sects. 5.1.2 and 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" 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 µ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 µ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 ~1 mJy from the VVDS 1.4 GHz observations with a resolution of 6" (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" 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,

$$ \frac{S_t}{S_p} = \frac{\sqrt{\theta_M^2 + \theta_b^2} \sqrt{\theta_m^2 + \theta_b^2}}{\theta_b^2}, \qquad (1) $$

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" 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:

1. **Elongated geometry**: sources are resolved in only one direction. This implies that $\theta_m = 0$ and

$$ \frac{S_t}{S_p} = \sqrt{\frac{\theta_M^2 + \theta_b^2}{\theta_b^2}}. \qquad (2) $$

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. **Circular geometry**: sources are uniformly resolved in all directions. This implies that $\theta_M = \theta_m$ and

$$ \frac{S_t}{S_p} = \frac{\theta_M^2 + \theta_b^2}{\theta_b^2}. \qquad (3) $$

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$ µJy, splitting them into two subranges: $S_t \le 40$ µJy, and $40 < S_t \le 100$ µ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\text{Jy}$. While some combinations of parameters provided a reasonable match for sources with $S_t \ge 40~\mu\text{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 \le 40~\mu\text{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\text{Jy}$ and sources with $40 < S_t \le 100~\mu\text{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 \le 40~\mu\text{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\text{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_{\min}$:

$$ \theta_{\min} = k_1 e^{-(S_t/k_2)^2}, \quad (4) $$

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_{\min}$ is required at flux densities $\ge 40~\mu\text{Jy}$. 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\text{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 Eqs. (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 \le 40~\mu\text{Jy}$) with $S_t/S_p < 1$ from 45% to 30% close to the

<figure><img src="image_15.png" /><figcaption>Fig. 15. Total-to-peak flux density ratio distributions in two total flux density ranges: $S_t \le 40~\mu\text{Jy}$ (top panel) and $40 < S_t \le 100~\mu\text{Jy}$ (bottom panel). Each panel shows the distribution of the observed sources (red histogram) that was derived from the 10 sets of simulations using an elongated geometry and $n = 0.6$, no minimum angular size (green histogram) and with a minimum angular size $\theta_{\min} = 0.3e^{-(S_t/40~\mu\text{Jy})^2}$ (blue histogram).</figcaption></figure>

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
<figure><img src="image_16.png" /><figcaption>**Fig. 16.** Completeness of our 3 GHz source catalog as a function of flux density and S/N. The mean completeness of all Monte Carlo runs (red line) and its standard deviation (gray shaded area) are shown. Also shown are the corrections when elongated (dash-dotted line) and circular (dashed line) geometries are assumed.</figcaption></figure>

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 (Sects. 5.1.2 and 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 \mu Jy$ ($S/N = 5.2$), and rather constant up to $20 \mu Jy$ ($S/N < 9$), beyond which they rise to a 94% completeness above $40 \mu Jy$ ($S/N \ge 16$). The mean error of the completeness and bias corrections is 5%. The two (elongated and circular geometry) approximations are consistent up to ~30 $\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 (~0.92 for circular and 0.98 for elongated morphology), implying average values of over 95% for fluxes higher than ~40 $\mu Jy$.

<table><caption>**Table 2.** Completeness and bias correction factors for the VLA-COSMOS 3 GHz catalog as a function of flux density.</caption><thead><tr><th>Flux density (μJy)</th><th>Completeness and bias correction factor (C<sub>compl</sub>)</th><th>Error</th></tr></thead><tbody><tr><td>&lt;10.4</td><td>0</td><td>—</td></tr><tr><td>11.0</td><td>0.08</td><td>0.01</td></tr><tr><td>11.6</td><td>0.40</td><td>0.02</td></tr><tr><td>12.3</td><td>0.55</td><td>0.03</td></tr><tr><td>13.0</td><td>0.58</td><td>0.03</td></tr><tr><td>13.8</td><td>0.56</td><td>0.04</td></tr><tr><td>14.6</td><td>0.57</td><td>0.03</td></tr><tr><td>15.5</td><td>0.57</td><td>0.04</td></tr><tr><td>16.4</td><td>0.57</td><td>0.04</td></tr><tr><td>17.3</td><td>0.57</td><td>0.04</td></tr><tr><td>18.4</td><td>0.59</td><td>0.04</td></tr><tr><td>19.4</td><td>0.56</td><td>0.04</td></tr><tr><td>21.1</td><td>0.68</td><td>0.05</td></tr><tr><td>23.3</td><td>0.73</td><td>0.04</td></tr><tr><td>25.8</td><td>0.82</td><td>0.05</td></tr><tr><td>28.6</td><td>0.85</td><td>0.05</td></tr><tr><td>31.7</td><td>0.89</td><td>0.06</td></tr><tr><td>35.1</td><td>0.91</td><td>0.07</td></tr><tr><td>38.8</td><td>0.94</td><td>0.06</td></tr><tr><td>43.0</td><td>0.92</td><td>0.06</td></tr><tr><td>47.6</td><td>0.95</td><td>0.07</td></tr><tr><td>53.9</td><td>0.94</td><td>0.05</td></tr><tr><td>62.4</td><td>0.95</td><td>0.07</td></tr><tr><td>72.2</td><td>0.93</td><td>0.06</td></tr><tr><td>83.5</td><td>0.92</td><td>0.08</td></tr><tr><td>96.7</td><td>0.97</td><td>0.08</td></tr><tr><td>&gt;100<sup>a</sup></td><td>1.00<sup>a</sup></td><td>0.05<sup>a</sup></td></tr></tbody></table>

Notes. (a) Assumed corrections for fluxes >100 $\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 densities of ~12–20 $\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 \ge 5$ across the entire observed field, 95 of which were outside the central two square degrees.

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''$ 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 µ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

<figure><img src="image_17.png" /><figcaption>**Fig. 17.** Fraction of false detections (red line) as a function of S/N (top panel) and flux density (bottom panel). The open (filled) histogram shows the number of components cataloged in the observed 3 GHz mosaic (detected in the inverted map), and limited to the central two square degrees. The data are also listed in Table 3.</figcaption></figure>

flux densities down to 30 µ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 h of on-source C-array configuration observations of one VLA pointing targeting the Lockman hole and reaching an rms of 1 µ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 µ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 µ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' 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 ~30 µJy with a discrepancy
<table><caption>**Table 3.** False detection probability as a function of S/N and flux density in the COSMOS two square degree field.</caption><thead><tr><th>S/N</th><th>Fraction</th><th>Flux density (μJy)</th><th>Fraction (F<sub>false-det</sub>)</th></tr></thead><tbody><tr><td>5.05</td><td>0.24</td><td>10.75</td><td>0.40</td></tr><tr><td>5.15</td><td>0.15</td><td>11.00</td><td>0.38</td></tr><tr><td>5.25</td><td>0.11</td><td>11.25</td><td>0.27</td></tr><tr><td>5.35</td><td>0.09</td><td>11.50</td><td>0.21</td></tr><tr><td>5.45</td><td>0.06</td><td>11.75</td><td>0.15</td></tr><tr><td>5.55</td><td>0.03</td><td>12.00</td><td>0.13</td></tr><tr><td>5.65</td><td>0.02</td><td>12.25</td><td>0.09</td></tr><tr><td>5.75</td><td>0.03</td><td>12.50</td><td>0.07</td></tr><tr><td>5.85</td><td>0.01</td><td>12.75</td><td>0.06</td></tr><tr><td>5.95</td><td>0.02</td><td>13.00</td><td>0.06</td></tr><tr><td>6.05</td><td>0.01</td><td>13.25</td><td>0.03</td></tr><tr><td>6.15</td><td>0.01</td><td>13.50</td><td>0.04</td></tr><tr><td>6.25</td><td>0.02</td><td>13.75</td><td>0.03</td></tr><tr><td>6.35</td><td>0.01</td><td>14.00</td><td>0.05</td></tr><tr><td>6.45</td><td>0.00</td><td>14.25</td><td>0.03</td></tr><tr><td>6.55</td><td>0.00</td><td>14.50</td><td>0.03</td></tr><tr><td>6.65</td><td>0.01</td><td>14.75</td><td>0.08</td></tr><tr><td>6.75</td><td>0.01</td><td>15.00</td><td>0.08</td></tr><tr><td>6.85</td><td>0.00</td><td>15.25</td><td>0.00</td></tr><tr><td>6.95</td><td>0.01</td><td>15.50</td><td>0.01</td></tr></tbody></table>

at fainter flux densities 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 deg² roughly corresponding to a circle with a radius of 5'. 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σ) 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

<figure><img src="image_18.png" /><figcaption>**Fig. 18.** Top panel: VLA-COSMOS 3 GHz Euclidean-normalized radio source counts, corrected using the completeness and bias and false-detection correction factors (black filled points) and without corrections (gray squares). Middle panel: VLA-COSMOS 3 GHz radio source counts compared to Condon et al. (2012) *P(D)* analysis with a single power law (dot-dashed red line) and Vernstrom et al. (2014) *P(D)* analysis with multiple power laws (green line) at 3 GHz. The yellow shaded area contains 95% of different source counts obtained from 100 square and nonoverlapping (8.5×8.5 arcmin²) subfields of the COSMOS field, thus demonstrating the effect of cosmic variance on fields with sizes similar to those analyzed by Condon et al. (2012) and Vernstrom et al. (2014). The dashed orange line shows the 68% interval of different source counts (obtained from 16th and 84th percentile in each flux density bin). Bottom panel: counts of the same sources, but shifted to the 1.4 GHz observed frame using a spectral index of $\alpha = -0.7$ prior to binning (black filled points). A selection of existing 1.4 GHz source counts in the literature is also shown, as indicated in the legend.</figcaption></figure>

(e.g., Condon et al. 2012). We show the 1.4 GHz source count comparison in the bottom panel of Fig. 18.
<table><caption>**Table 4.** Radio source counts at 3 GHz within the COSMOS two square degree field, normalized to Euclidean geometry.</caption><thead><tr><th>Flux density (mJy)</th><th>Counts<sup>a</sup> (Jy<sup>1.5</sup> sr<sup>-1</sup>)</th><th>Error<sup>b</sup> (Jy<sup>1.5</sup> sr<sup>-1</sup>)</th><th>N</th><th>Correction factor<sup>a</sup></th></tr></thead><tbody><tr><td>0.011</td><td>0.866</td><td>0.068</td><td>631</td><td>3.27</td></tr><tr><td>0.013</td><td>0.952</td><td>0.056</td><td>1109</td><td>1.64</td></tr><tr><td>0.015</td><td>1.10</td><td>0.078</td><td>991</td><td>1.70</td></tr><tr><td>0.018</td><td>1.16</td><td>0.094</td><td>849</td><td>1.67</td></tr><tr><td>0.020</td><td>1.40</td><td>0.11</td><td>888</td><td>1.54</td></tr><tr><td>0.024</td><td>1.38</td><td>0.086</td><td>811</td><td>1.33</td></tr><tr><td>0.028</td><td>1.47</td><td>0.10</td><td>780</td><td>1.18</td></tr><tr><td>0.032</td><td>1.57</td><td>0.12</td><td>702</td><td>1.12</td></tr><tr><td>0.037</td><td>1.59</td><td>0.13</td><td>587</td><td>1.08</td></tr><tr><td>0.045</td><td>1.58</td><td>0.13</td><td>753</td><td>1.07</td></tr><tr><td>0.059</td><td>1.53</td><td>0.12</td><td>505</td><td>1.05</td></tr><tr><td>0.076</td><td>1.67</td><td>0.14</td><td>366</td><td>1.08</td></tr><tr><td>0.098</td><td>1.60</td><td>0.17</td><td>250</td><td>1.03</td></tr><tr><td>0.13</td><td>1.65</td><td>0.15</td><td>181</td><td>1.00</td></tr><tr><td>0.16</td><td>1.62</td><td>0.17</td><td>121</td><td>1.00</td></tr><tr><td>0.21</td><td>1.82</td><td>0.21</td><td>93</td><td>1.00</td></tr><tr><td>0.27</td><td>1.93</td><td>0.25</td><td>67</td><td>1.00</td></tr><tr><td>0.35</td><td>2.03</td><td>0.31</td><td>48</td><td>1.00</td></tr><tr><td>0.78</td><td>3.95</td><td>0.37</td><td>159</td><td>1.00</td></tr><tr><td>3.0</td><td>10.4</td><td>1.5</td><td>56</td><td>1.00</td></tr><tr><td>11</td><td>37.4</td><td>7.4</td><td>27</td><td>1.00</td></tr><tr><td>44</td><td>72.3</td><td>28</td><td>7</td><td>1.00</td></tr></tbody></table>

**Notes.** (a) The listed counts were corrected for completeness and bias ($C_{compl}$), as well as false detection fractions ($F_{false-det}$), by multiplying the raw counts by the correction factor given in the last column, and equal to $(1-F_{false-det})/C_{compl}$ (see Tables 2 and 3). (b) The source count errors take into account only the Poissonian errors and completeness and bias correction uncertainties (see text for details).

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. 2015) 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 (≤2′′) 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 ~200 µJy, but are slightly lower in the flux density range of 100–200 µ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. 2017).

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 ~1.6′′ and rms ≈ 2.7 µJy beam⁻¹. 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 ~200 µ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 (≤80 µ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 ~80 µ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 ~100 µ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 ~50–300 µJy, while it overpredicts the counts at flux densities ≤30 µ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_{2-10 keV} \sim 10^{42}-10^{44}$ erg s⁻¹; 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. 2008; Bonzini et al. 2013; Smolčić et al. 2017; Delvecchio et al. 2017). This could explain the lack of sources with flux densities in the range of ~50–300 µJy in the model (see, e.g., Smolčić et al. 2008; Padovani et al. 2015; Smolčić et al. 2017), 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
<table><caption>**Table 5.** Catalog sample page.</caption><thead><tr><th>ID</th><th>NAME</th><th>RA<br>(J2000, deg)</th><th>RA_ERR<br>(arcsec)</th><th>DEC<br>(J2000, deg)</th><th>DEC_ERR<br>(arcsec)</th><th>FLUX<br>(µJy)</th><th>FLUX_ERR<br>(µJy)</th><th>RMS<br>(µJy beam<sup>-1</sup>)</th><th>SNR</th><th>NPIX</th><th>RES</th><th>MULTI</th></tr></thead><tbody><tr><td>78</td><td>COSMOSVLA3 J095709.33+020940.7</td><td>149.288886</td><td>0.01</td><td>2.161331</td><td>0.01</td><td>13400.0</td><td>670.0</td><td>28.7</td><td>385.0</td><td>104</td><td>1</td><td>0</td></tr><tr><td>1110</td><td>COSMOSVLA3 J095709.83+015457.4</td><td>149.290996</td><td>0.021</td><td>1.915946</td><td>0.021</td><td>1190.0</td><td>64.0</td><td>22.1</td><td>29.0</td><td>102</td><td>1</td><td>0</td></tr><tr><td>5144</td><td>COSMOSVLA3 J095710.49+013644.7</td><td>149.29372</td><td>0.065</td><td>1.612418</td><td>0.065</td><td>145.0</td><td>19.0</td><td>17.4</td><td>8.32</td><td>20</td><td>0</td><td>0</td></tr><tr><td>3749</td><td>COSMOSVLA3 J095710.53+025132.2</td><td>149.29389</td><td>0.053</td><td>2.858966</td><td>0.053</td><td>345.0</td><td>32.0</td><td>26.8</td><td>10.4</td><td>31</td><td>1</td><td>0</td></tr><tr><td>10099</td><td>COSMOSVLA3 J095710.57+022657.0</td><td>149.294072</td><td>0.1</td><td>2.449174</td><td>0.1</td><td>92.0</td><td>18.0</td><td>17.4</td><td>5.3</td><td>11</td><td>0</td><td>0</td></tr><tr><td>4979</td><td>COSMOSVLA3 J095710.88+020929.8</td><td>149.295339</td><td>0.076</td><td>2.158297</td><td>0.076</td><td>2820.0</td><td>140.0</td><td>21.4</td><td>7.13</td><td>332</td><td>1</td><td>0</td></tr><tr><td>3366</td><td>COSMOSVLA3 J095711.09+023031.3</td><td>149.296231</td><td>0.048</td><td>2.508722</td><td>0.048</td><td>318.0</td><td>27.0</td><td>21.5</td><td>11.4</td><td>35</td><td>1</td><td>0</td></tr><tr><td>8753</td><td>COSMOSVLA3 J095711.15+021104.2</td><td>149.296459</td><td>0.095</td><td>2.184514</td><td>0.095</td><td>121.0</td><td>22.0</td><td>21.2</td><td>5.68</td><td>18</td><td>0</td><td>0</td></tr><tr><td>4046</td><td>COSMOSVLA3 J095711.64+021236.8</td><td>149.298529</td><td>0.055</td><td>2.210248</td><td>0.055</td><td>223.0</td><td>26.0</td><td>22.6</td><td>9.83</td><td>22</td><td>0</td><td>0</td></tr><tr><td>6546</td><td>COSMOSVLA3 J095711.67+021401.3</td><td>149.298655</td><td>0.077</td><td>2.233706</td><td>0.077</td><td>135.0</td><td>21.0</td><td>19.3</td><td>7.02</td><td>18</td><td>0</td><td>0</td></tr><tr><td>10942</td><td>COSMOSVLA3 J110034.76+014635.7</td><td>150.144846</td><td>-99.0</td><td>1.776607</td><td>-99.0</td><td>374.0</td><td>-99.0</td><td>2.4</td><td>-99.0</td><td>828</td><td>1</td><td>1</td></tr><tr><td>9229</td><td>COSMOSVLA3 J110034.78+025027.4</td><td>150.144947</td><td>0.1</td><td>2.840947</td><td>0.1</td><td>25.6</td><td>2.7</td><td>2.4</td><td>5.18</td><td>25</td><td>1</td><td>0</td></tr><tr><td>8777</td><td>COSMOSVLA3 J110034.80+021421.1</td><td>150.145009</td><td>0.095</td><td>2.239195</td><td>0.095</td><td>13.0</td><td>2.4</td><td>2.27</td><td>5.7</td><td>14</td><td>0</td><td>0</td></tr><tr><td>489</td><td>COSMOSVLA3 J110034.81+025515.6</td><td>150.145042</td><td>0.013</td><td>2.921023</td><td>0.013</td><td>229.0</td><td>13.0</td><td>3.66</td><td>62.5</td><td>49</td><td>0</td><td>0</td></tr><tr><td>439</td><td>COSMOSVLA3 J110034.83+014247.1</td><td>150.145157</td><td>0.013</td><td>1.713091</td><td>0.013</td><td>162.0</td><td>9.0</td><td>2.27</td><td>71.3</td><td>51</td><td>0</td><td>0</td></tr><tr><td>2850</td><td>COSMOSVLA3 J110034.94+020234.9</td><td>150.145599</td><td>0.042</td><td>2.043054</td><td>0.042</td><td>29.9</td><td>2.8</td><td>2.29</td><td>13.1</td><td>29</td><td>0</td><td>0</td></tr><tr><td>2439</td><td>COSMOSVLA3 J110034.99+024524.5</td><td>150.145815</td><td>0.038</td><td>2.756821</td><td>0.038</td><td>32.7</td><td>2.8</td><td>2.22</td><td>14.7</td><td>29</td><td>0</td><td>0</td></tr><tr><td>2837</td><td>COSMOSVLA3 J110035.00+024614.3</td><td>150.145836</td><td>0.041</td><td>2.770643</td><td>0.041</td><td>39.6</td><td>2.9</td><td>2.18</td><td>13.3</td><td>41</td><td>1</td><td>0</td></tr><tr><td>4092</td><td>COSMOSVLA3 J110035.05+024154.6</td><td>150.146065</td><td>0.055</td><td>2.69852</td><td>0.055</td><td>28.3</td><td>2.6</td><td>2.23</td><td>9.93</td><td>31</td><td>1</td><td>0</td></tr><tr><td>4148</td><td>COSMOSVLA3 J110035.07+020350.5</td><td>150.14614</td><td>0.057</td><td>2.064044</td><td>0.057</td><td>23.4</td><td>2.7</td><td>2.43</td><td>9.64</td><td>23</td><td>0</td><td>0</td></tr><tr><td>...</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>5351</td><td>COSMOSVLA3 J110034.81+015420.5</td><td>150.940886</td><td>0.067</td><td>1.905719</td><td>0.067</td><td>117.0</td><td>16.0</td><td>14.6</td><td>8.05</td><td>18</td><td>0</td><td>0</td></tr><tr><td>5225</td><td>COSMOSVLA3 J110034.48+023458.3</td><td>150.943703</td><td>0.067</td><td>2.582882</td><td>0.067</td><td>127.0</td><td>17.0</td><td>15.5</td><td>8.14</td><td>21</td><td>0</td><td>0</td></tr><tr><td>8447</td><td>COSMOSVLA3 J110034.52+022031.2</td><td>150.943867</td><td>0.11</td><td>2.342001</td><td>0.11</td><td>290.0</td><td>20.0</td><td>14.4</td><td>5.04</td><td>48</td><td>1</td><td>0</td></tr><tr><td>10677</td><td>COSMOSVLA3 J110034.56+015500.1</td><td>150.944027</td><td>0.11</td><td>1.916707</td><td>0.11</td><td>75.7</td><td>15.0</td><td>14.9</td><td>5.1</td><td>11</td><td>0</td><td>0</td></tr><tr><td>6533</td><td>COSMOSVLA3 J110034.63+022415.8</td><td>150.944306</td><td>0.076</td><td>2.404402</td><td>0.076</td><td>143.0</td><td>16.0</td><td>14.2</td><td>7.09</td><td>28</td><td>1</td><td>0</td></tr><tr><td>2195</td><td>COSMOSVLA3 J110034.12+022510.6</td><td>150.946345</td><td>0.036</td><td>2.419631</td><td>0.036</td><td>548.0</td><td>31.0</td><td>15.3</td><td>15.6</td><td>83</td><td>1</td><td>0</td></tr><tr><td>10084</td><td>COSMOSVLA3 J110034.27+020117.6</td><td>150.94699</td><td>0.1</td><td>2.021575</td><td>0.1</td><td>73.1</td><td>15.0</td><td>14.1</td><td>5.18</td><td>16</td><td>0</td><td>0</td></tr><tr><td>9695</td><td>COSMOSVLA3 J110034.51+023539.7</td><td>150.947979</td><td>0.1</td><td>2.594384</td><td>0.1</td><td>93.2</td><td>18.0</td><td>17.6</td><td>5.3</td><td>10</td><td>0</td><td>0</td></tr><tr><td>10652</td><td>COSMOSVLA3 J110034.80+024951.9</td><td>150.949175</td><td>0.11</td><td>2.831096</td><td>0.11</td><td>80.7</td><td>17.0</td><td>16.1</td><td>5.01</td><td>14</td><td>0</td><td>0</td></tr><tr><td>10484</td><td>COSMOSVLA3 J110034.85+021102.3</td><td>150.952242</td><td>0.11</td><td>2.183987</td><td>0.11</td><td>82.3</td><td>17.0</td><td>16.3</td><td>5.05</td><td>13</td><td>0</td><td>0</td></tr></tbody></table>
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 h 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 \ \mu\text{Jy} \text{ beam}^{-1}$ over the two square degrees, at an angular resolution of 0.75". 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" 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 $\mu\text{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 quasi stellar objects 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 Framework programs under grant agreements 333654 (CIG, "AGN feedback") and 337595 (ERC Starting Grant, "CoSMass"). M.B. and P.C. acknowledge support from the PRIN-INAF 2014. A.K. and F.B. acknowledge support by the Collaborative Research Council 956, sub-project A1, funded by the Deutsche Forschungsgemeinschaft (DFG). B.M. and F.B. acknowledge support through DFG priority program 1573 founded by the DFG.

## References

* Afonso, J., Georgakakis, A., Almeida, C., et al. 2005, ApJ, 624, 135

* Aihara, H., Allende Prieto, C., An, D., et al. 2011, ApJS, 193, 29

* Aird, J., Coil, A. L., Moustakas, J., et al. 2012, ApJ, 746, 90

* Aretxaga, I., Wilson, G. W., Aguilar, E., et al. 2011, MNRAS, 415, 3831

* Becker, R. H., White, R. L., & Helfand, D. J. 1995, ApJ, 450, 559

* Bertoldi, F., Carilli, C., Aravena, M., et al. 2007, ApJS, 172, 132

* Best, P. N., Kaiser, C. R., Heckman, T. M., & Kauffmann, G. 2006, MNRAS, 368, L67

* Béthermin, M., Daddi, E., Magdis, G., et al. 2012, ApJ, 757, L23

* Bock, D. C.-J., Large, M. I., & Sadler, E. M. 1999, AJ, 117, 1578

* Bondi, M., Ciliegi, P., Zamorani, G., et al. 2003, A&A, 403, 857

* Bondi, M., Ciliegi, P., Venturi, T., et al. 2007, A&A, 463, 519

* Bondi, M., Ciliegi, P., Schinnerer, E., et al. 2008, ApJ, 681, 1129

* Bonzini, M., Mainieri, V., Padovani, P., et al. 2012, ApJS, 203, 15

* Bonzini, M., Padovani, P., Mainieri, V., et al. 2013, MNRAS, 436, 3759

* Bourke, S., Mooley, K., & Hallinan, G. 2014, in Astronomical Data Analysis Software and Systems XXIII, eds. N. Manset, & P. Forshay, ASP Conf. Ser., 485, 367

* Bouwens, R. J., Illingworth, G. D., Franx, M., & Ford, H. 2007, ApJ, 670, 928

* Bower, R. G., Benson, A. J., Malbon, R., et al. 2006, MNRAS, 370, 645

* Capak, P., Aussel, H., Ajiki, M., et al. 2007, ApJS, 172, 99

* Ciliegi, P., McMahon, R. G., Miley, G., et al. 1999, MNRAS, 302, 222

* Civano, F., Marchesi, S., Comastri, A., et al. 2016, ApJ, 819, 62

* Condon, J. J. 1984, ApJ, 287, 461

* Condon, J. J. 1992, ARA&A, 30, 575

* Condon, J. J. 2015, ArXiv e-prints [arXiv:1502.05616]

* Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693

* Condon, J. J., Cotton, W. D., Yin, Q. F., et al. 2003, AJ, 125, 2411

* Condon, J. J., Cotton, W. D., Fomalont, E. B., et al. 2012, ApJ, 758, 23

* Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11

* Delvecchio, I., Smolčić, V., Zamorani, G., et al. 2017, A&A, 602, A3 (VLA-COSMOS SI)

* de Zotti, G., Massardi, M., Negrello, M., & Wall, J. 2010, A&ARV, 18, 1

* Dickinson, M., Giavalisco, M., & GOODS Team 2003, in The Mass of Galaxies at Low and High Redshift, eds. R. Bender, & A. Renzini, Proc. ESO Astrophysics Symp., 324

* Driver, S. P., Norberg, P., Baldry, I. K., et al. 2009, Astron. Geophys., 50, 12

* Driver, S. P., Hill, D. T., Kelvin, L. S., et al. 2011, MNRAS, 413, 971

* Elvis, M., Civano, F., Vignali, C., et al. 2009, ApJS, 184, 158

* Evans, D. A., Worrall, D. M., Hardcastle, M. J., Kraft, R. P., & Birkinshaw, M. 2006, ApJ, 642, 96

* Feigelson, E. D., & Nelson, P. I. 1985, ApJ, 293, 192

* Fixsen, D. J., Kogut, A., Levin, S., et al. 2011, ApJ, 734, 5

* Georgakakis, A., Mobasher, B., Cram, L., et al. 1999, MNRAS, 306, 708

* Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35

* Haarsma, D. B., Partridge, R. B., Windhorst, R. A., & Richards, E. A. 2000, ApJ, 544, 641

* Hales, C. A., Murphy, T., Curran, J. R., et al. 2012, MNRAS, 425, 979

* Hales, C. A., Norris, R. P., Gaensler, B. M., et al. 2014, MNRAS, 441, 2555

* Hardcastle, M., Evans, D., & Croston, J. 2007, MNRAS, 376, 1849

* Hasinger, G., Cappelluti, N., Brunner, H., et al. 2007, ApJS, 172, 29

* Hopkins, A., Georgakakis, A., Cram, L., Afonso, J., & Mobasher, B. 2000, ApJS, 128, 469

* Hopkins, A. M., Afonso, J., Chan, B., et al. 2003, AJ, 125, 465

* Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55

* Jarvis, M. J. 2012, African Skies, 16, 44

* Karim, A., Schinnerer, E., Martínez-Sansigre, A., et al. 2011, ApJ, 730, 61

* Kimball, A. E., & Ivezić, Ž. 2008, AJ, 136, 684

* Koekemoer, A. M., Aussel, H., Calzetti, D., et al. 2007, ApJS, 172, 196

* Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36

* Krishna, G., Sirothia, S. K., Mhaskey, M., et al. 2014, MNRAS, 443, 2824

* Laigle, C., McCracken, H. J., Ilbert, O., et al. 2011, ApJS, 224, 24

* Le Fèvre, O., Tasca, L. A. M., Cassata, P., et al. 2015, A&A, 576, A79

* Lilly, S. J., Le Fèvre, O., Renzini, A., et al. 2007, ApJS, 172, 70

* Lilly, S. J., Le Brun, V., Maier, C., et al. 2009, ApJS, 184, 218

* Lisenfeld, U., & Völk, H. J. 2000, A&A, 354, 423

* Magnelli, B., Elbaz, D., Chary, R. R., et al. 2011, A&A, 528, A35

* McCracken, H. J., Milvang-Jensen, B., Dunlop, J., et al. 2012, A&A, 544, A156

* McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Data Analysis Software and Systems XVI, eds. R. A. Shaw, F. Hill, & D. J. Bell, ASP Conf. Ser., 376, 127

* Miettinen, O., Smolčić, V., Novak, M., et al. 2015, A&A, 577, A29

* Miller, N. A., Fomalont, E. B., Kellermann, K. I., et al. 2008, ApJS, 179, 114

* Miller, N. A., Bonzini, M., Fomalont, E. B., et al. 2013, ApJS, 205, 13

* Mooley, K. P., Hallinan, G., Bourke, S., et al. 2016, ApJ, 818, 105

* Mullaney, J. R., Alexander, D. M., Goulding, A. D., & Hickox, R. C. 2011, MNRAS, 414, 1082

* Mullaney, J. R., Daddi, E., Béthermin, M., et al. 2012, ApJ, 753, L30

* Norris, R. P., Huynh, M. T., Jackson, C. A., et al. 2005, AJ, 130, 1358
Norris, R. P., Hopkins, A. M., Afonso, J., et al. 2011, PASA, 28, 215

Norris, R. P., Afonso, J., Bacon, D., et al. 2013, PASA, 30, 20

Norris, R., Basu, K., Brown, M., et al. 2015, in Proc. Advancing Astrophysics with the Square Kilometre Array (AASKA14), 86

Novak, M., Smolčić, V., Civano, F., et al. 2015, MNRAS, 447, 1282

Novak, M., Smolčić, V., Delhaize, J., et al. 2017, A&A, 602, A5 (VLA-COSMOS SI)

Owen, F. N., & Morrison, G. E. 2008, AJ, 136, 1889

Owen, F. N., Morrison, G. E., Klimek, M. D., & Greisen, E. W. 2009, AJ, 137, 4846

Padovani, P. 2011, MNRAS, 411, 1547

Padovani, P., Mainieri, V., Tozzi, P., et al. 2009, ApJ, 694, 235

Padovani, P., Bonzini, M., Kellermann, K. I., et al. 2015, MNRAS, 452, 1263

Prandoni, I., & Seymour, N. 2015, in Proc. Advancing Astrophysics with the Square Kilometre Array (AASKA14), 67

Prandoni, I., Gregorini, L., Parma, P., et al. 2001, A&A, 365, 392

Prescott, M. K. M., Impey, C. D., Cool, R. J., & Scoville, N. Z. 2006, ApJ, 644, 100

Rau, U., & Cornwell, T. J. 2011, A&A, 532, A71

Rees, M. J. 1967, MNRAS, 136, 279

Rodighiero, G., Daddi, E., Baronchelli, I., et al. 2011, ApJ, 739, L40

Sanders, D. B., Salvato, M., Aussel, H., et al. 2007, ApJS, 172, 86

Sargent, M. T., Béthermin, M., Daddi, E., & Elbaz, D. 2012, ApJ, 747, L31

Schinnerer, E., Carilli, C. L., Scoville, N. Z., et al. 2004, AJ, 128, 1974

Schinnerer, E., Smolčić, V., Carilli, C. L., et al. 2007, ApJS, 172, 46

Schinnerer, E., Sargent, M. T., Bondi, M., et al. 2010, ApJS, 188, 384

Schmitt, J. H. M. M. 1985, ApJ, 293, 178

Scott, K. S., Austermann, J. E., Perera, T. A., et al. 2008, MNRAS, 385, 2225

Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1

Seymour, N., Dwelly, T., Moss, D., et al. 2008, MNRAS, 386, 1695

Smolčić, V. 2009, ApJ, 699, L43

Smolčić, V., & Riechers, D. A. 2011, ApJ, 730, 64

Smolčić, V., Schinnerer, E., Scodeggio, M., et al. 2008, ApJS, 177, 14

Smolčić, V., Zamorani, G., Schinnerer, E., et al. 2009a, ApJ, 696, 24

Smolčić, V., Schinnerer, E., Zamorani, G., et al. 2009b, ApJ, 690, 610

Smolčić, V., Aravena, M., Navarrete, F., et al. 2012, A&A, 548, A4

Smolčić, V., Ciliegi, P., Jelić, V., et al. 2014, MNRAS, 443, 2590

Smolčić, V., Padovani, P., Delhaize, J., et al. 2015, in Proc. Advancing Astrophysics with the Square Kilometre Array (AASKA14), 69

Smolčić, V., Delvecchio, G., Zamorani, G., et al. 2017, A&A, 602, A2 (VLA-COSMOS SI)

Tasse, C., Röttgering, H. J. A., Best, P. N., et al. 2007, A&A, 471, 1105

Trump, J. R., Impey, C. D., McCarthy, P. J., et al. 2007, ApJS, 172, 383

Vernstrom, T., Scott, D., & Wall, J. V. 2011, MNRAS, 415, 3641

Vernstrom, T., Scott, D., Wall, J. V., et al. 2014, MNRAS, 440, 2791

Wilman, R. J., Miller, L., Jarvis, M. J., et al. 2008, MNRAS, 388, 1335