# The 1.28 GHz MeerKAT DEEP2 Image

T. Mauch, W. D. Cotton, J. J. Condon, A. M. Matthews, T. D. Abbott, R.

M. Adam, M. A. Aldera, K. M. B. Asad, E. F. Bauermeister, T. G. H.

Bennett, et al.

► To cite this version:

T. Mauch, W. D. Cotton, J. J. Condon, A. M. Matthews, T. D. Abbott, et al. The 1.28 GHz MeerKAT DEEP2 Image. The Astrophysical Journal, 2020, 888, 10.3847/1538-4357/ab5d2d. insu-03720945

HAL Id: insu-03720945

https://insu.hal.science/insu-03720945v1

Submitted on 7 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

The 1.28 GHz MeerKAT DEEP2 Image

T. Mauch<sup>1,10</sup>, W. D. Cotton<sup>1,2</sup>, J. J. Condon<sup>2,10</sup>, A. M. Matthews<sup>2,3,10</sup>, T. D. Abbott<sup>1</sup>, R. M. Adam<sup>1</sup>, M. A. Aldera<sup>4</sup>, K. M. B. Asad<sup>1,5,6</sup>, E. F. Bauermeister<sup>1</sup>, T. G. H. Bennett<sup>1</sup>, H. Bester<sup>1</sup>, D. H. Botha<sup>7</sup>, L. R. S. Brederode<sup>1</sup>, Z. B. Brits<sup>1</sup>, S. J. Buchner<sup>1</sup>, J. P. Burger<sup>1</sup>, F. Camilo<sup>1,10</sup>, J. M. Chalmers<sup>1</sup>, T. Cheetham<sup>1</sup>, D. de Villiers<sup>8</sup>, M. S. de Villiers<sup>1</sup>, M. A. Dikgale-Mahlakoana<sup>1</sup>, L. J. du Toit<sup>7</sup>, S. W. P. Esterhuyse<sup>1</sup>, G. Fadana<sup>1</sup>, B. L. Fanaroff<sup>1</sup>, S. Fataar<sup>1</sup>, S. February<sup>1</sup>, B. S. Frank<sup>1</sup>, R. R. G. Gamatham<sup>1</sup>, M. Geyer<sup>1</sup>, S. Goedhart<sup>1</sup>, S. Gounden<sup>1</sup>, S. C. Gumede<sup>1</sup>, I. Heywood<sup>5,9</sup>, M. J. Hlakola<sup>1</sup>, J. M. G. Horrell<sup>10</sup>, B. Hugo<sup>1,5</sup>, A. R. Isaacson<sup>1</sup>, G. I. G. Józsa<sup>1,5</sup>, J. L. Jonas<sup>1,5</sup>, R. P. M. Julie<sup>1</sup>, F. B. Kapp<sup>1</sup>, V. A. Kasper<sup>1</sup>, J. S. Kenyon<sup>5</sup>, P. P. A. Kotzé<sup>1</sup>, N. Kriek<sup>1</sup>, H. Kriel<sup>1</sup>, T. W. Kusel<sup>1</sup>, R. Lehmensiek<sup>7</sup>, A. Loots<sup>11</sup>, R. T. Lord<sup>1</sup>, B. M. Lunsky<sup>1</sup>, K. Madisa<sup>1</sup>, L. G. Magnus<sup>1</sup>, J. P. L. Main<sup>1</sup>, J. A. Malan<sup>1</sup>, J. R. Manley<sup>1</sup>, S. J. Marais<sup>7</sup>, A. Martens<sup>1</sup>, B. Merry<sup>1</sup>, R. Millenaar<sup>1</sup>, N. Mnyandu<sup>1</sup>, I. P. T. Moeng<sup>1</sup>, O. J. Mokone<sup>1</sup>, T. E. Monama<sup>1</sup>, M. C. Mphego<sup>1</sup>, W. S. New<sup>1</sup>, B. Ngcebetsha<sup>1,5</sup>, K. J. Ngoasheng<sup>1</sup>, M. T. O. Ockards<sup>1</sup>, N. Oozeer<sup>1</sup>, A. J. Otto<sup>1</sup>, A. A. Patel<sup>1</sup>, A. Peens-Hough<sup>1</sup>, S. J. Perkins<sup>1</sup>, A. J. T. Ramaila<sup>1,5</sup>, Z. R. Ramudzuli<sup>1</sup>, R. Renil<sup>1</sup>, L. L. Richter<sup>1</sup>, A. Robyntjies<sup>1</sup>, S. Salie<sup>1</sup>, C. T. G. Schollar<sup>1</sup>, L. C. Schwartz<sup>1</sup>, M. Serylak<sup>1</sup>, R. Siebrits<sup>1</sup>, S. K. Sirothia<sup>1,5</sup>, O. M. Smirnov<sup>1,5</sup>, L. Sofeya<sup>1</sup>, G. Stone<sup>1</sup>, B. Taljaard<sup>1</sup>, C. Tasse<sup>5,12</sup>, I. P. Theron<sup>7</sup>, A. J. Tiplady<sup>1</sup>, O. Toruvanda<sup>1</sup>, S. N. Twum<sup>1</sup>, T. J. van Balla<sup>1</sup>, A. van der Byl<sup>1</sup>, C. van der Merwe<sup>1</sup>, V. Van Tonder<sup>1</sup>, B. H. Wallace<sup>1</sup>, M. G. Welz<sup>1</sup>, L. P. Williams<sup>1</sup>, and B. Xaia<sup>1</sup>

<sup>1</sup>South African Radio Astronomy Observatory, 2 Fir Street, Black River Park, Observatory 7925, South Africa; tmauch@ska.ac.za

<sup>2</sup>National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA

<sup>3</sup>Department of Astronomy, University of Virginia, Charlottesville, VA 22904, USA

<sup>4</sup>Tellumat (Pty) Ltd. 64-74 White Road, Retreat 7945, South Africa

<sup>5</sup>Department of Physics and Electronics, Rhodes University, P.O. Box 94, Grahamstown 6140, South Africa

<sup>6</sup>Department of Physics and Astronomy, University of the Western Cape, Bellville, Cape Town 7535, South Africa

<sup>7</sup>EMSS Antennas, 18 Techno Avenue, Technopark, Stellenbosch 7600, South Africa

<sup>8</sup>Department of Electrical and Electronic Engineering, Stellenbosch University, Stellenbosch 7600, South Africa

<sup>9</sup>Oxford Astrophysics, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK

<sup>10</sup>Inter-University Institute for Data-Intensive Astronomy, University of Cape Town, Private Bag X3, Rondebosch 7701, South Africa

<sup>11</sup>Presidential Infrastructure Coordinating Commission, 77 Meintjies Street, Sunnyside, Pretoria 0001, South Africa

<sup>12</sup>GEPI, Observatoire de Paris, CNRS, PSL Research University, Université Paris Diderot, F-92190, Meudon, France
Received 2019 August 15; revised 2019 November 11; accepted 2019 November 26; published 2020 January 8

Abstract

We present the confusion-limited 1.28 GHz MeerKAT DEEP2 image covering one $\theta_b \approx 68'$ FWHM primary-beam area with $\theta = 7''$ 6 FWHM resolution and $\sigma_n = 0.55 \pm 0.01$ $\mu$Jy beam$^{-1}$ rms noise. Its J2000 center position $\alpha = 04^h 13^m 26.4^s$, $\delta = -80^\circ 00' 00''$ was selected to minimize artifacts caused by bright sources. We introduce the new 64-element MeerKAT array and describe commissioning observations to measure the primary-beam attenuation pattern, estimate telescope pointing errors, and pinpoint $(u, v)$ coordinate errors caused by offsets in frequency or time. We constructed a 1.4 GHz differential source count by combining a power-law count fit to the DEEP2 confusion $P(D)$ distribution from 0.25 to 10 $\mu$Jy with counts of individual DEEP2 sources between 10 $\mu$Jy and 2.5 mJy. Most sources fainter than $S \sim 100$ $\mu$Jy are distant star-forming galaxies (SFGs) obeying the far-IR/radio correlation, and sources stronger than 0.25 $\mu$Jy account for $\sim 93\%$ of the radio background produced by SFGs. For the first time, the DEEP2 source count has reached the depth needed to reveal the majority of the star formation history of the universe. A pure luminosity evolution of the 1.4 GHz local luminosity function consistent with the Madau & Dickinson model for the evolution of SFGs based on UV and infrared data underpredicts our 1.4 GHz source count in the range $-5 \lesssim \log[S(Jy)] \lesssim -4$.

Unified Astronomy Thesaurus concepts: Radio telescopes (1360); Galaxy counts (588); Star formation (1569)

1. Introduction

The extragalactic source population at 1.4 GHz is a mixture of galaxies with active galactic nuclei (AGNs) and star-forming galaxies (SFGs; Condon & Broderick 1988; Afonso et al. 2005; Simpson et al. 2006; Padovani et al. 2015). Radio sources powered by AGNs account for nearly all of the strong-source population, and SFGs whose radio emission primarily comes from synchrotron electrons accelerated by the supernova remnants of massive ($M > 8M_\odot$) stars (Condon 1992) dominate below $S \sim 100$ $\mu$Jy at 1.4 GHz (Simpson et al. 2006; Bonzini et al. 2013; Prandoni et al. 2018).

The far-infrared (FIR)/radio correlation obeyed by nearly all SFGs indicates that their 1.4 GHz luminosities are directly proportional to their star formation rates (SFRs; Condon 1992).

Dust is transparent at 1.4 GHz, so sufficiently sensitive radio continuum observations could trace the cosmic evolution of the mean star formation rate density (SFRD) unbiased by dust emission or absorption.

In the past two decades a number of wide-area redshift surveys (York et al. 2000; Colless et al. 2001; Jones et al. 2009) used in combination with mJy-sensitivity radio surveys such as the NRAO VLA Sky Survey (NVSS; Condon et al. 1998) and the Sydney University Molonglo Sky Survey (SUMSS; Mauch et al. 2003) have allowed accurate determinations of the local radio luminosity functions (RLFs) for both SFGs and AGNs (Sadler et al. 2002; Best et al. 2005; Mauch & Sadler 2007; Condon et al. 2019). In all cases the two populations were classified using available optical, mid-infrared, FIR, and radio data.
Measuring the evolving SFRD by directly determining RLFs at higher redshifts and comparing them with the well-determined local RLF is difficult. It requires deep multi-wavelength data to identify and classify the host galaxies of faint radio sources, plus photometric or spectroscopic redshifts. Recent studies of SFRD evolution by this method (e.g., Smolčić et al. 2009; Padovani et al. 2011; McAlpine et al. 2013) suggest pure luminosity evolution $\propto (1+z)^{2.5}$ for the most luminous SFGs at $z \lesssim 2.5$. However, current radio surveys are not sensitive enough to detect the fainter galaxies responsible for the bulk of star formation around “cosmic noon” at $z \sim 2$. Most current samples are further hampered by their reliance on deep multiwavelength data covering very small solid angles in fields selected using only optical/infrared criteria, which can be less than ideal for making deep radio images. A $5\sigma_n \approx 0.25$ $\mu$Jy beam$^{-1}$ sensitivity is needed to detect the individual SFGs accounting for most of the star formation history of the universe (SFHU), so making a survey with rms noise $\sigma_n \approx 0.05$ $\mu$Jy beam$^{-1}$ is one of the key continuum science goals of the Square Kilometre Array (SKA; Jarvis et al. 2015; Prandoni & Seymour 2015).

Confusion by sources blended in the synthesized beams of deep continuum images also limits the ability of radio surveys to detect faint SFGs. For example, the MeerKAT MIGHTEE survey will have rms confusion $\sigma_c \approx 2$ $\mu$Jy beam$^{-1}$ in its $\theta \approx 8''$ synthesized beam (Jarvis et al. 2016). The detection threshold for individual sources will be about 17 $\mu$Jy, the level at which there are ~25 beam solid angles per source and below which the association of fitted components with individual galaxies is hampered by the increasing level of obscuration by stronger ones. Recent work by Condon et al. (2012) and Vernstrom et al. (2014) has shown that statistical analysis of the confusion brightness distribution, usually called the $P(D)$ distribution, can be used to estimate the source count at sub-$\mu$Jy levels. Such very deep source counts, combined with the already accurately determined local RLF, can be used to constrain the SFHU directly (Condon & Matthews 2018).

This paper presents the 1.28 GHz DEEP2 image observed with the South African Radio Astronomy Observatory's (SARAO) MeerKAT array, counts of radio sources in DEEP2, and a preliminary analysis of the SFHU constrained by those counts. Section 2 provides a brief introduction to the MeerKAT array and describes our early commissioning observations designed to measure and improve its performance. In order to maximize the depth we can reach with the finite dynamic range of the MeerKAT array, we chose the DEEP2 field to be as free as possible from bright radio sources in the MeerKAT primary beam, as described in Section 3. Our observing strategy for the DEEP2 field plus our method of calibration and imaging the raw data are outlined in Section 4. Section 5 presents source counts from 0.25 to 10 $\mu$Jy derived from the DEEP2 image $P(D)$ distribution and between 10 $\mu$Jy and 2.5 mJy based on discrete sources in DEEP2. Finally, for Section 6 we evolved the local 1.4 GHz luminosity function of SFGs (Condon et al. 2019) for pure luminosity evolution of the Madau & Dickinson (2014) fit to the UV/FIR SFHU and compared it with our 1.4 GHz source count.

## 2. The MeerKAT Radio Telescope

The MeerKAT array was used to observe the DEEP2 field. MeerKAT is a precursor to the SKA located in the Karoo region of South Africa's Northern Cape province. MeerKAT

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Distribution of all $N_{\text{bl}} = 2016$ MeerKAT baseline lengths, which range between 29 and 7698 m. Half of the baselines are between the 48 antennas in the densely packed 1 km diameter core.</figcaption></figure>

was inaugurated in 2018 July and, during the course of our DEEP2 observations, was in its early commissioning phase. MeerKAT is a new instrument that has little presence in the literature, so we describe the telescope features (Section 2.1) and results from our early commissioning observations (Section 2.2) relevant for deep continuum observations.

### 2.1. MeerKAT

MeerKAT is an array of 64 13.5 m diameter dish antennas spread out over roughly 8 km centered near latitude 30°42' S and longitude 21°23' E. Each dish has a 3.8 m offset Gregorian subreflector and a receiver indexer located just below the subreflector to ensure a completely unblocked aperture. A conical skirt extends below the subreflector to deflect radiation from the surrounding ground away from the receiver. The unblocked aperture is essential for deep continuum imaging at L band because it improves dynamic range by (1) lowering the sensitivity of primary-beam sidelobes to strong sources and radio frequency interference (RFI) outside the main beam and (2) reducing the system noise temperature by limiting pickup of ground radiation. The observations described in this paper were all carried out with the dual linear polarization (horizontal and vertical) L-band (856–1712 MHz) receivers (Lehmensiek & Theron 2012, 2014). Each antenna has a measured system noise temperature $T_{\text{sys}} \approx 20$ K and a remarkably low system equivalent flux density SEFD $\approx 430$ Jy on cold sky.

The MeerKAT antennas are named “m000” through “m063,” the order of which roughly follows their distances from the array center. The inner 48 antennas are located within a 1 km diameter central “core” region, and the remaining 16 are spread beyond this central area out to a radius of nearly 4 km. The distribution of baseline lengths between the 2016 antenna pairs (Figure 1) is peaked at lengths shorter than 1000 m, and roughly half of all MeerKAT baselines are between antennas in the core. MeerKAT and large proposed arrays such as the SKA and ngVLA have fixed antennas in centrally concentrated multiscale configurations, which are compromises designed to satisfy conflicting demands for high surface brightness sensitivity and high angular resolution. MeerKAT provides excellent L-band surface brightness sensitivity on angular scales $\gtrsim 1'$ at the cost of a very broad naturally weighted synthesized beam (Figure 2). To achieve $\theta \approx 7''6$ FWHM resolution in our DEEP2 image, we observed only when most of the outer antennas were working, and we gave up some
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Two slices through synthesized beams derived from a simulated 12 hr MeerKAT observation at $\delta = -80^{\circ}$. The dashed line shows a beam calculated using natural weights (i.e., constant per visibility) and has FWHM $\theta = 23''$.3. The solid line shows a beam calculated with uniform weights (i.e., inversely proportional to the number of visibilities in a grid sample) and has FWHM $\theta = 3''2$.</figcaption></figure>

sensitivity by heavily downweighting the (u, v) data from intracore baselines.

Visibilities are transported to the archive located in the Centre for High Performance Computing 600 km away in Cape Town, where they are converted to formats used by major radio astronomy imaging and analysis packages (e.g., Measurement Set or AIPS UV). Additional information about MeerKAT and its specifications can be found in Jonas & MeerKAT Team (2016) and Camilo et al. (2018).

## 2.2. Commissioning Observations of PKS B1934-638

To verify the accuracy of our early MeerKAT data, we made a series of snapshot images offset by 10', 20', 30', 40', and 50' to the north, south, east, and west of the calibration source PKS B1934-638. At the $\theta \approx 8''$ resolution of MeerKAT, PKS B1934-638 is a point source. It is also strong ($S \approx 15.1$ Jy) at 1284 MHz and has a well-established radio spectrum (Reynolds 1994). The usefulness of such offset snapshots is described in Condon et al. (1998): position errors can reveal incorrect (u, v) coordinates or frequency labeling in the data, and the variation in amplitude with position in the primary beam can be used to measure the primary-beam attenuation patterns and pointing errors of the MeerKAT dishes.

A 14-minute observation cycling through the offset pointings was made with a start time chosen to ensure that the scans were as close to the transit of PKS B1934-638 as possible (azimuth $\approx 180^{\circ}$, elevation $\approx 57^{\circ}$). Images made from individual pointings have rms noise $\sigma_n \sim 1$ mJy beam$^{-1}$. The flux density and position of the source were measured in each image by fitting an elliptical Gaussian using the AIPS task JMFIT. At all offsets the source is never attenuated to $S_a < 2.3$ Jy beam$^{-1}$ by the MeerKAT primary beam, so the source signal-to-noise ratio is always ($S/N) \gtrsim 2300:1$. The noise component of fitting errors should be <0.006 in position and <0.04% in flux density assuming a point source and a circular 8" beam (Condon 1997), so errors in the measured positions and flux densities are dominated by calibration errors.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Commanded position offsets (arcmin) for the observations of PKS B1934-638 described in Section 2.2 are shown as plus signs. The red points show the measured positions calculated from data with frequency and time errors; the open circles show the measured positions after correction for these errors. The differences between the measured and commanded offsets as shown have been magnified 600× to highlight the small errors in the measured positions. The inset in the upper right corner shows the distribution of corrected R.A. and decl. errors in units of arcsec for all 21 pointings.</figcaption></figure>

### 2.2.1. Positions

Geometric errors can be introduced into source positions on the image plane by incorrect calculations of the (u, v) coordinates associated with measured visibilities. The (u, v) coordinates used during imaging were calculated at a reference frequency of $\nu_0 = 886$ MHz. If $\nu_0$ differs from the actual array reference frequency $\nu_c$, this will rescale the image radial offsets of source positions from the phase center $\rho_m$ from their true offsets $\rho$ (Condon et al. 1998). Furthermore, if the time stamps used to calculate the (u, v) coordinates differ from the true time stamps of observation by $\Delta t$ s, an image near the celestial pole will be rotated about its phase center by an angle $\theta \approx 2\pi \Delta t/T_d$, where $T_d \approx 86,164$ s is one sidereal day.

Figure 3 compares the measured offsets (circles and red points) of the 21 pointings with their commanded offsets (plus signs). The red points indicate the measured offsets in our initial data set before any corrections; they clearly show both a rotation and a radial scaling. The rotation revealed a 2 s shift in the time stamps of the raw correlated visibilities. The radial scaling was caused by the frequency labeling being offset by 0.5 channels during conversion of the data to AIPS UV. Correcting these errors moved the red points to the black circles in Figure 3. To ensure that all MeerKAT data do not require these time and frequency corrections, the errors have now been fixed in the MeerKAT correlator and in the software.

The mean ratio of the corrected measured (circles) to the commanded (plus signs) radial offsets is consistent with unity:

$$ \langle \rho_m / \rho \rangle = \nu_c / \nu_o = 1 - (4.0 \pm 4.9) \times 10^{-5}. \quad (1) $$

Likewise, the circles in Figure 3 are consistent with no rotation about the center.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Blue curve shows a horizontal slice through holographic measurements of the 1.5 GHz Stokes *I* primary-beam power pattern of MeerKAT (M. de Villiers 2019, in preparation). The black curve is the attenuation pattern calculated from Equation (3) for FWHM $\theta_b = 57.5$.</figcaption></figure>

The upper right inset of Figure 3 shows the distribution of position errors in R.A. and decl. Their standard errors and means are

$$ 
\begin{aligned}
\sigma_{\Delta\alpha} &= 0''11, & \langle\Delta\alpha\rangle &= 0''001, \\
\sigma_{\Delta\delta} &= 0''06, & \langle\Delta\delta\rangle &= 0''02.
\end{aligned} 
\qquad (2) 
$$

Thus, we expect that individual strong sources within ~50' of the phase centers of *L*-band MeerKAT images will have rms position errors of ~0''1 in each coordinate, and the image frames will have astrometric uncertainties of ~0''01.

### 2.2.2. The MeerKAT Primary Beam

The attenuated flux densities of PKS B1934–638, observed at various pointing offsets $\rho$, divided by its flux density at the pointing center, measure the primary-beam attenuation pattern $a_b(\rho) \equiv S(\rho)/S(\rho=0)$ of the MeerKAT antennas. The attenuation pattern derived from these data and also from a horizontal slice through holographic measurements of the MeerKAT Stokes *I* beam at 1.5 GHz (M. de Villiers 2019, in preparation) is well matched by the attenuation pattern resulting from cosine-tapered field (or cosine-squared power) illumination (Condon & Ransom 2016):

$$ a_b(\rho/\theta_b) = \left[ \frac{\cos(1.189\pi\rho/\theta_b)}{1 - 4(1.189\rho/\theta_b)^2} \right]^2. \qquad (3) $$

For the purpose of comparing the observed attenuation pattern at 1.5 GHz with Equation (3), we set the FWHM of the horizontal slice through the primary power pattern to

$$ \theta_b = 57.5 \left( \frac{\nu}{1.5 \text{ GHz}} \right)^{-1}. \qquad (4) $$

Figure 4 shows that the attenuation pattern of the cosine illumination taper (Equation (3)) is a good match out to $\rho = 2.5$. Equation (3) also has the practical advantage of expressing $a_b(\rho/\theta_b)$ as an elementary function, so we used it to fit our PKS B1934–638 data and in all subsequent analyses requiring narrowband beam patterns at frequency $\nu$.

Figure 5 shows fits of Equation (3) to the flux densities measured from our PKS B1934–638 observations. During fitting, we inserted the mean pointing offset $\Delta$ as a free parameter, replacing $(\rho/\theta_b)$ in Equation (3) by $[(\rho - \Delta)/\theta_b]$.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Fitted attenuation patterns of Equation (3) to vertical (top panel) and horizontal (bottom panel) slices through the Stokes *I* MeerKAT *L*-band primary beam. Plus signs show the measured attenuation in the flux density *S* of PKS B1934–638 derived from elliptical Gaussian fits to the source at each offset position from our observations. $\theta_b$ denotes the best-fitting FWHMs at $\nu = 1.28$ GHz, and $\Delta$ denotes the offsets of the peak in the fitted attenuation patterns from $\rho = 0$.</figcaption></figure>

<table><caption>**Table 1**
Frequency Dependence of the MeerKAT Primary Beamwidths $\theta_b$ and Typical Pointing Errors $\Delta$</caption><thead><tr><th rowspan="2">Subband Number</th><th rowspan="2">ν (MHz)</th><th colspan="2">Vertical</th><th colspan="2">Horizontal</th></tr><tr><th>θ<sub>b</sub> (arcmin)</th><th>Δ (arcsec)</th><th>θ<sub>b</sub> (arcmin)</th><th>Δ (arcsec)</th></tr></thead><tbody><tr><td>1</td><td>908.04</td><td>100.1</td><td>-17.1</td><td>96.2</td><td>-21.7</td></tr><tr><td>2</td><td>952.34</td><td>94.7</td><td>-5.1</td><td>91.1</td><td>-24.1</td></tr><tr><td>3</td><td>996.65</td><td>90.5</td><td>-2.4</td><td>87.1</td><td>-23.4</td></tr><tr><td>4</td><td>1043.46</td><td>86.1</td><td>6.2</td><td>82.8</td><td>-25.4</td></tr><tr><td>5</td><td>1092.78</td><td>81.7</td><td>14.0</td><td>78.6</td><td>-27.8</td></tr><tr><td>6</td><td>1144.61</td><td>78.2</td><td>17.4</td><td>75.2</td><td>-27.0</td></tr><tr><td>7</td><td>1198.94</td><td>73.4</td><td>21.8</td><td>70.5</td><td>-32.2</td></tr><tr><td>8</td><td>1255.79</td><td>70.0</td><td>21.1</td><td>67.3</td><td>-29.5</td></tr><tr><td>9</td><td>1317.23</td><td>65.7</td><td>19.7</td><td>63.2</td><td>-25.9</td></tr><tr><td>10</td><td>1381.18</td><td>63.1</td><td>12.1</td><td>60.6</td><td>-32.6</td></tr><tr><td>11</td><td>1448.05</td><td>60.6</td><td>-29.7</td><td>58.3</td><td>-68.3</td></tr><tr><td>12</td><td>1519.94</td><td>58.9</td><td>-30.2</td><td>56.8</td><td>-63.6</td></tr><tr><td>13</td><td>1593.92</td><td>56.2</td><td>-19.4</td><td>54.3</td><td>-39.4</td></tr><tr><td>14</td><td>1656.20</td><td>55.4</td><td>-50.4</td><td>53.6</td><td>-43.2</td></tr></tbody></table>

The horizontal and vertical fits of the Equation (3) attenuation power pattern are an excellent match to the data for all $\log_{10}(a_b) > -0.70$.

The primary-beam pattern is slightly elliptical; its vertical $\theta_b$ is larger than its horizontal $\theta_b$. Simulations of the MeerKAT antenna optics predict that asymmetries in the horizontal and vertical linearly polarized *L*-band feeds result in an ellipticity in the Stokes *I* primary-beam pattern. Our measurement of the ellipticity as a function of frequency (Table 1) agrees with these simulations to within 1%. We conclude that feed asymmetry is the reason for the ellipticity in the measured Stokes *I* MeerKAT *L*-band primary beam.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Variation of $\theta_b\nu$ with frequency $\nu$ in GHz from the values tabulated in Table 1. The lower line shows values from fits of Equation (3) to a horizontal slice through the primary beam, and the upper line shows fits to a vertical slice.</figcaption></figure>

We used wideband images of PKS B1934–638 covering the frequency range 886–1682 MHz centered on $\nu = 1.28$ GHz to fit the primary beam in Figure 5. We have further split the band into 14 narrow subbands and repeated the beam fitting described above for each. Table 1 summarizes the results. Note that the narrow subbands defined in this table are the same as the subbands in Table 5 used for the wideband imaging of the DEEP2 field.

Figure 6 shows that the product $\theta_b\nu$ varies by <±3% across the band for both the vertical and horizontal cuts through the beam. The best L-band approximations with fixed $\theta_b\nu$ are

$$ \theta_b = 89.5 \left( \frac{\nu}{\text{GHz}} \right)^{-1} \quad \text{(Vertical)} $$

$$ \theta_b = 86.2 \left( \frac{\nu}{\text{GHz}} \right)^{-1} \quad \text{(Horizontal).} \quad (5) $$

### 3. DEEP2 Field Selection

The finite dynamic range of MeerKAT limits the minimum rms fluctuation that can be achieved because every deep L-band image contains thousands of sources distributed over the primary-beam area. The dynamic range of such images is defined by the ratio of the effective source flux density $S_{\text{eff}}$ to the rms fluctuation $\sigma$ in regions devoid of sources:

$$ \text{DR} \equiv \frac{S_{\text{eff}}}{\sigma}, \quad (6) $$

where $S_{\text{eff}}$ is the quadratic sum of attenuated flux densities $S_a$ of all sources in the primary beam (Condon 2009)

$$ S_{\text{eff}} = \left( \sum_i S_{a,i}^2 \right)^{\frac{1}{2}}. \quad (7) $$

The attenuated flux density of a source with true flux density $S$ offset from the pointing center by an angle $\rho$ is

$$ S_a \equiv a_b(\rho)S, \quad (8) $$

where $a_b(\rho)$ is the primary-beam attenuation pattern approximated by a circular Gaussian. The quadratic sum over $S_a$ in Equation (7) is primarily determined by the strongest sources in

<table><caption>**Table 2**
Five Minimum-demerit Positions with $d < 1.4$ mJy</caption><thead><tr><th>α<br>(J2000)</th><th>δ<br>(mJy)</th><th>d<br>(mJy)</th><th>Ω (d &lt; 1.4 mJy)<br>(deg²)</th></tr></thead><tbody><tr><td>03<sup>h</sup>21<sup>m</sup></td><td>-18°53′</td><td>1.32</td><td>0.19</td></tr><tr><td>04<sup>h</sup>22<sup>m</sup></td><td>-80°15′</td><td>1.35</td><td>0.22</td></tr><tr><td>16<sup>h</sup>37<sup>m</sup></td><td>-70°46′</td><td>1.34</td><td>0.36</td></tr><tr><td>21<sup>h</sup>04<sup>m</sup></td><td>-54°25′</td><td>1.36</td><td>0.25</td></tr><tr><td>22<sup>h</sup>03<sup>m</sup></td><td>-35°43′</td><td>1.34</td><td>0.23</td></tr></tbody></table>

the field of view, so the best area for a deep field has the fewest and faintest bright sources in the primary beam.

The attenuated flux density $S_a$ of a source varies with telescope pointing errors and receiver gain fluctuations, both of which contribute to $\sigma$. Pointing errors contribute a flux density error

$$ \Delta S_p(\rho) \approx \epsilon(\rho)S_a, \quad (9) $$

where

$$ \epsilon(\rho) \approx \frac{8 \ln(2) \rho \sigma_p}{\theta_b^2} \quad (10) $$

is the fractional attenuation change at an angle $\rho$ from the pointing center caused by an rms pointing error $\sigma_p$ (Condon 2009). Receiver gain fluctuations cause a flux density change

$$ \Delta S_g(\rho) \approx \sigma_g S_a, \quad (11) $$

where $\sigma_g$ is the rms receiver gain error.

Each source in an image contributes an error flux density that is the quadratic sum of these flux density fluctuations:

$$ \Delta S = \sqrt{\Delta S_p^2 + \Delta S_g^2}. \quad (12) $$

Image artifacts produced by individual sources are independent, so we define an image “demerit” score $d$ equal to the quadratic sum of the independent source contributions to $\Delta S$ in the primary beam:

$$ d \equiv \left( \sum_i \Delta S_i^2 \right)^{\frac{1}{2}}. \quad (13) $$

To locate the best deep fields observable by MeerKAT, we searched the SUMSS (Mauch et al. 2003) catalog south of $\delta = -35^\circ$ and the NVSS (Condon et al. 1998) catalog from $\delta \ge -35^\circ$ to $\delta = +10^\circ$ over a fine grid of $\approx 6 \times 10^7$ potential pointings separated by 1.05 in the 18,185 deg² area defined by $\delta < 10^\circ$, $|b| > 10^\circ$. At each grid position we computed $d$ from flux densities shifted to 1.28 GHz assuming $S \propto \nu^{-0.7}$ and with parameters conservatively appropriate to MeerKAT: $\theta_b = 68'$, $\sigma_p = 30''$, and $\sigma_g = 0.01$ (Jonas & MeerKAT Team 2016) out to the radius $\rho = 3^\circ$ extending beyond the second sidelobe of the MeerKAT primary beam.

The five pointing centers with the smallest demerit flux densities $d$ in the southern hemisphere are listed in Table 2, along with their solid angles $\Omega$ in which $d < 1.4$ mJy. We chose the southernmost field at J2000 $\alpha = 04^h 22^m$, $\delta = -80^\circ 15'$ to ensure that observations of the field could be easily scheduled at most times of the day during the early commissioning and engineering phase of the telescope. Also,
<table><caption>**Table 3**
DEEP2 Observation Summary</caption><thead><tr><th>Date</th><th>Start Time<br>UTC</th><th>τ<sub>Total</sub><br>(hr)</th><th>τ<sub>Target</sub></th><th>N<sub>Ants</sub></th></tr></thead><tbody><tr><td>2018 Apr 27</td><td>07:11</td><td>11.0</td><td>8.4</td><td>61</td></tr><tr><td>2018 Jun 30</td><td>23:00</td><td>16.2</td><td>12.6</td><td>60</td></tr><tr><td>2018 Jul 7</td><td>21:39</td><td>17.2</td><td>13.4</td><td>61</td></tr><tr><td>2018 Jul 16</td><td>21:37</td><td>8.0</td><td>6.0</td><td>61</td></tr><tr><td>2018 Jul 24</td><td>21:07</td><td>8.9</td><td>6.9</td><td>59</td></tr><tr><td>2018 Jul 25</td><td>21:01</td><td>9.0</td><td>7.6</td><td>58</td></tr><tr><td>2018 Jul 27</td><td>21:01</td><td>16.1</td><td>14.0</td><td>61</td></tr><tr><td>2018 Jul 28</td><td>20:51</td><td>16.2</td><td>14.1</td><td>60</td></tr><tr><td>2018 Oct 8</td><td>21:33</td><td>9.5</td><td>8.5</td><td>59</td></tr><tr><td>2018 Nov 4</td><td>14:37</td><td>16.2</td><td>14.2</td><td>62</td></tr><tr><td>2019 Jan 19</td><td>09:31</td><td>16.1</td><td>13.9</td><td>63</td></tr><tr><td>2019 Jan 20</td><td>09:41</td><td>10.8</td><td>9.2</td><td>63</td></tr></tbody><tfoot><tr><td>Total</td><td></td><td>155.2</td><td>128.8</td><td></td></tr></tfoot></table>

at ecliptic latitude $\beta \approx -75^{\circ}$, the DEEP2 field is easily observed by orbiting telescopes and is minimally affected by zodiacal dust. We inspected a mosaic image made in 2017 with a 16-antenna MeerKAT subarray and covering a 2 deg$^2$ region surrounding our selected position, and we finally chose the field centered at J2000 $\alpha = 04^{\mathrm{h}}13^{\mathrm{m}}26^{\mathrm{s}}4$, $\delta = -80^{\circ}00'00''$ (which we call DEEP2) in order to move a few moderately bright extended sources farther from the pointing center. Our DEEP2 field has a demerit score $d = 1.4$ mJy, only slightly higher than the 1.35 mJy minimum at J2000 $\alpha = 04^{\mathrm{h}}22^{\mathrm{m}}$, $\delta = -80^{\circ}15'$.

# 4. Observations and Imaging

## 4.1. DEEP2 Observations

The DEEP2 field centered on J2000 $\alpha = 04^{\mathrm{h}}13^{\mathrm{m}}26^{\mathrm{s}}4$, $\delta = -80^{\circ}00'00''$ was observed at L band in 12 separate sessions between 2018 April 27 and 2019 January 20 for a total of 155.2 hr (Table 3). We always required that at least 58 of the 64 antennas and at least seven of the nine outer-ring antennas (those providing the longest baselines) be available. This ensured sufficient long- and intermediate-baseline coverage to produce a fairly clean “dirty beam” point-spread function (PSF) with $\theta \lesssim 8''$ FWHM resolution. We preferentially observed during the night, though sessions with longer duration and other scheduling constraints meant that ~30% of the observations occurred in the daytime. The -80° decl. of the DEEP2 field ensures that it is never <55° from the Sun.

Scans on the DEEP2 target lasted 15 minutes and were interleaved with scans on the S(1284 MHz) ≈ 6.1 Jy phase and secondary gain calibrator PKS J0252-7104 located ≈10° from the DEEP2 field center. Scans on PKS J0252-7104 were 2 minutes long before July 25, after which we shortened them to 1 minute because we found that we were getting sufficient S/ Ns on the calibrator gain solutions. The primary flux density and bandpass calibrator PKS B1934-638 was observed for 10 minutes at the start of each observation and then subsequently every 3 hr until it set. Its assumed flux densities from Reynolds (1994) are listed as a function of frequency in Table 4.

Table 3 summarizes our observations of the DEEP2 field. From the total 155.2 hr, calibration/slewing overheads took 17% and left 128.8 hr on the DEEP2 target. Our observations

<table><caption>**Table 4**
The Eight Calibration Subbands</caption><thead><tr><th>Subband<br>Number</th><th>ν<br>(MHz)</th><th>S(PKS B1934-638)<sup>a</sup><br>(Jy)</th></tr></thead><tbody><tr><td>1</td><td>935.728</td><td>14.516</td></tr><tr><td>2</td><td>1035.204</td><td>14.921</td></tr><tr><td>3</td><td>1134.680</td><td>15.108</td></tr><tr><td>4</td><td>1234.157</td><td>15.129</td></tr><tr><td>5</td><td>1333.634</td><td>15.028</td></tr><tr><td>6</td><td>1433.110</td><td>14.835</td></tr><tr><td>7</td><td>1532.586</td><td>14.577</td></tr><tr><td>8</td><td>1632.063</td><td>14.428</td></tr></tbody></table>

**Note.**

had an integration period of 8 s except for the initial April 27 observation that we averaged from 4 to 8 s prior to any calibration. Observations were all carried out in the 4096 × 208.984 kHz channel L-band continuum mode. The raw data volume was typically ≈2 TB for a 12 hr observation. Each data set was calibrated separately prior to imaging.

The uncalibrated visibilities from our observations are publicly available and were obtained from the SARAO archive at https://archive.sarao.ac.za. Readers interested in obtaining these data can do so by searching the archive for Proposal ID: SCI-20180426-TM-01.

## 4.2. Editing and Calibration

The full MeerKAT L band covers the frequency range 856-1712 MHz with 4096 spectral channels. We trimmed 144 channels each from the lower and upper ends of the full band because the receiver response is too weak at the band edges to be useful. The remaining frequency range of our data is 886-1682 MHz. This band is contaminated by various sources of strong RFI, the broadest of which are difficult to detect automatically. We therefore developed an empirical mask that flags at all times those channels most contaminated by RFI. Figure 7 shows an example of the raw MeerKAT bandpass on two baselines, one short (319 m) and one long (7566 m). The long baselines are typically not as badly affected by RFI as the short ones, so we chose to apply the mask only to baselines shorter than 1000 m. The mask rejects 34% of the trimmed L band, or 17% of the full data set when applied only to the subset of short baselines.

The unmasked data were searched in time and frequency for deviations and flagged. Our flagging method is similar to the SumThreshold technique described by Offringa et al. (2010). A smooth background was fitted to the data by convolving them with a Gaussian whose widths in frequency and time are larger-than-expected RFI spike widths and are smaller than any variations in the bandpass or changes in amplitude with time. The smoothed background was then subtracted from the data, and outliers in the residuals were detected in increasing averaged widths in both time and frequency.

After the data were masked and initially edited, we used the Obit package (Cotton 2008) for further editing and calibration. Prior to calibration, the raw data were smoothed with a Hanning filter to prevent Gibbs ringing. This filter combines adjacent channels with weights (1/4, 1/2, 1/4) and effectively doubles the channel width to $2 \times 208.984$ kHz = 417.968 kHz. It makes neighboring

<footnote>a Flux densities from Reynolds (1994).</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Average amplitude of cross-hand polarization visibilities (Horizontal × Vertical) from a 10-minute scan on PKS B1934−638 on two baselines. We chose to plot data from a cross-hand polarization, as this is more sensitive to polarized RFI signals. A short baseline (m000 × m010; 319 m) is plotted in the top panel, and a long one (m059 × m063; 7566 m) is shown in the bottom panel. The gray shaded areas in the top panel show the regions masked for all times on baselines shorter than 1000 m.</figcaption></figure>

frequency channels degenerate, so we excised every second channel from the smoothed data prior to calibration.

Our calibration consisted of the following steps:

1. The variations in group delays were calculated from the observations of PKS B1934−638 and PKS J0252−7104 and were interpolated to all of the data.

2. A bandpass calibration was determined to remove residual variations in gain and phase as a function of frequency. This was based on a CLEAN component model within 1° of PKS B1934−638. The average correction from all scans on PKS B1934−638 was applied to the data.

3. The observations of PKS J0252−7104 were used to correct the phases and amplitudes on the target for each calibration subband as a function of time. The amplitudes of PKS B1934−638 were derived from the model of Reynolds (1994) in each subband (see Column (3) of Table 4) and used to determine the amplitude spectrum of PKS J0252−7104. Amplitude and phase corrections were then interpolated in time and applied to the entire data set.

4. Some residual errors that were not detected during earlier editing were found by searching for gains with amplitudes that are discrepant by more than 20σ and flagging them.

5. The fully calibrated data were edited once more. At this stage the few visibilities with extremely discrepant Stokes *I* amplitudes (>200 Jy) and outliers from a running median in time and frequency were flagged.

6. The previous steps were repeated after resetting the calibration while retaining the flags on the calibrated data.

After the editing steps listed above, ~35% of the data were flagged. Prior to imaging, the calibrated visibilities were

spectrally averaged again, to 2 × 417.968 kHz = 835.936 kHz. The calibrated, flagged, and averaged data volume for the target was ≈200 GB in a typical 12 hr observation.

### 4.3. Self-calibration and Final Editing

The DEEP2 data were collected during 12 observing sessions spread over 9 months (Table 3) and were individually phase-referenced to PKS J0252−7104, which is ~10° from the DEEP2 field center at J2000 $\alpha = 04^h 13^m 26.4^s 4$, $\delta = -80^\circ 00' 00''$. The resulting phases were corrected to the direction and times of the target observations and astrometrically aligned with each other before imaging. The data from one day (2018 July 27) were imaged using two iterations of phase-only self-calibration with a 30 s averaging time. The resulting model of the DEEP2 field was used as the initial model to start the phase self-calibration of all data sets. Starting the self-calibrations of other days from one model ensures that images from all days are astrometrically aligned. We avoided amplitude self-calibration for two reasons: (1) amplitude self-calibration does not work well in a field containing many faint sources and no dominant point source, and (2) the external gain calibration based on observations of PKS B1934−638 worked very well. We measured the flux density of the $S \approx 12$ mJy point source at J2000 $\alpha = 04^h 15^m 08.2^s 21$, $\delta = -79^\circ 59'41.0''$ on the 12 daily DEEP2 images (Table 3); its rms variation is only 2% over the full 9-month observation period.

The data from each session were then imaged and deconvolved without further self-calibration to ensure good data quality. As a final editing step, the models derived from these preliminary images were Fourier transformed and subtracted from the calibrated (u, v) data. Residual amplitudes >0.5 Jy in the difference data were flagged in the calibrated session data. Next, the data were time averaged in a baseline-dependent fashion using Obit/UVBIAvg with the constraints of <1% amplitude loss within 1°.5 of the field center and averaging time <30 s. Finally, the averaged data sets were concatenated by Obit/UVAppend into a single data set for imaging.

### 4.4. Imaging

The concatenated data set was imaged using the Obit task MFImage (Cotton et al. 2018) without further self-calibration. MFImage used small planar facets to cover the wide field of view and made separate images in the 14 imaging subbands (Table 5) having fractional bandwidths $\Delta\nu/\nu < 0.05$ small enough to accommodate the frequency dependence of sky brightness and primary-beam attenuation. Dirty/residual images were formed in each subband, but a weighted average image was used to drive the CLEAN process. CLEAN components included the pixel flux densities from each subband, and the subtraction during the major cycles used a spectral index fitted to each component to interpolate between the subband center frequencies listed in Table 5. The joint deconvolution of the subband images requires that the width of the dirty PSF be nearly independent of frequency. This was accomplished by a frequency-dependent (u, v) taper that downweighted the longer baselines at the higher frequencies.

The facet images were reprojected during gridding to form a coherent grid of pixels on the plane tangent to the celestial
<table><caption>**Table 5**
DEEP2 Imaging Subband Frequencies and Weights</caption><thead><tr><th>Subband Number</th><th>ν<sub>i</sub> (MHz)</th><th>σ<sub>n,i</sub> (μJy beam<sup>-1</sup>)</th><th>w<sub>i</sub> for Min σ<sub>n</sub><sup>2</sup></th><th>w<sub>i</sub> for Max S/N</th></tr></thead><tbody><tr><td>i = 1</td><td>908.04</td><td>4.22</td><td>0.0225</td><td>0.0378</td></tr><tr><td>2</td><td>952.34</td><td>5.04</td><td>0.0158</td><td>0.0248</td></tr><tr><td>3</td><td>996.65</td><td>3.20</td><td>0.0393</td><td>0.0580</td></tr><tr><td>4</td><td>1043.46</td><td>2.88</td><td>0.0483</td><td>0.0669</td></tr><tr><td>5</td><td>1092.78</td><td>2.76</td><td>0.0526</td><td>0.0683</td></tr><tr><td>6</td><td>1144.61</td><td>2.58</td><td>0.0603</td><td>0.0733</td></tr><tr><td>7</td><td>1198.94</td><td>4.20</td><td>0.0227</td><td>0.0259</td></tr><tr><td>8</td><td>1255.79</td><td>3.98</td><td>0.0253</td><td>0.0271</td></tr><tr><td>9</td><td>1317.23</td><td>1.85</td><td>0.1171</td><td>0.1170</td></tr><tr><td>10</td><td>1381.18</td><td>1.64</td><td>0.1486</td><td>0.1389</td></tr><tr><td>11</td><td>1448.05</td><td>1.55</td><td>0.1672</td><td>0.1463</td></tr><tr><td>12</td><td>1519.94</td><td>1.87</td><td>0.1147</td><td>0.0938</td></tr><tr><td>13</td><td>1593.92</td><td>2.89</td><td>0.0481</td><td>0.0368</td></tr><tr><td>14</td><td>1656.20</td><td>1.85</td><td>0.1173</td><td>0.0850</td></tr></tbody></table>

sphere at the pointing center. This allows a joint CLEAN of many facets in a given major cycle.

The MeerKAT antenna array is centrally concentrated, so a relatively strong ROBUST weighting (ROBUST = -1.3 in AIPS/Obit usage) was used to downweight the shortest baselines to give a $\theta = 7.6''$ FWHM synthesized beamwidth. The DEEP2 field was completely imaged out to a radius of 1.5, and facets out to 2° were added as needed to cover outlying strong sources selected from the SUMSS (Mauch et al. 2003) catalog. The source density in our DEEP2 image is so high that no CLEAN windowing was used.

CLEAN used a loop gain of 0.15 and found 250,000 point components stronger than 7 µJy for a total CLEAN flux density of $S = 1.301$ Jy. Prior to restoring the CLEAN components with a circular Gaussian of FWHM $\theta = 7.6''$, the residuals in each subband and facet were convolved with a Gaussian with widths calculated to give a dirty PSF having approximately the same FWHM as the restoring beam. First, the CLEAN components appearing in each facet of each subband image were restored, and then the subband facets were collected into a single subband plane. The output of this imaging process is an image cube containing the 14 CLEANed and restored subband images.

### 4.5. Wideband Images

We took weighted averages of the subband images, all of which have well-defined center frequencies $\nu_i$ and small fractional bandwidths $\Delta\nu/\nu_i \lesssim 0.05$, to produce single-plane wideband images using two different weighting schemes. If the rms noise in each subband image is $\sigma_i$, then subband weights $w_i \propto \sigma_i^{-2}$ minimize the wideband image noise variance

$$ \sigma_n^2 = \frac{\sum_{i=1}^{14} w_i \sigma_i^2}{\sum_{i=1}^{14} w_i}. \qquad (14) $$

More generally, subband weights $w_i \propto (\nu_i^\alpha/\sigma_i)^2$ maximize the wideband image S/N for sources with spectral index $\alpha \equiv +d \ln(S)/d \ln(\nu)$:

$$ (S/N)^2 \propto \frac{\sum_{i=1}^{14} w_i (\nu_i^\alpha / \sigma_i)^2}{\sum_{i=1}^{14} w_i}. \qquad (15) $$

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Distribution of peak flux densities of the S/N-weighted wideband image in five apparently source-free regions covering 14,373 synthesized beam solid angles of $\sim 1.6''$ from the pointing center is shown by the thick curve. It is well matched by the Gaussian with rms $\sigma_n = 0.55$ µJy beam$^{-1}$ shown plotted as a dashed curve. The Gaussian is a parabola on this semilogarithmic plot.</figcaption></figure>

Minimizing $\sigma_n$ (Equation (14)) is equivalent to choosing $\alpha = 0$ in Equation (15). The median spectral index of faint sources selected at frequencies $\nu \sim 1.4$ GHz is $\langle \alpha \rangle \approx -0.7$ (Donald 1984), so this is the best choice of $\alpha$ for maximizing the S/N. The first three columns of Table 5 list the subband numbers, center frequencies $\nu_i$ (MHz), and rms noise values $\sigma_i$ (µJy beam$^{-1}$). Column (4) tabulates the weights $w_i$ that minimize $\sigma_n^2$, and Column (5) shows the different weights $w_i$ that maximize the S/N in the wideband DEEP2 image. Both sets of weights have been normalized to make $\sum_{i=1}^{14} w_i = 1$ for convenience.

If a source near the pointing center has flux density $S \propto \nu^\alpha$, its flux density in the weighted wideband image will be

$$ S \propto \frac{\sum_{i=1}^{14} w_i \nu_i^\alpha}{\sum_{i=1}^{14} w_i}. \qquad (16) $$

We define the “effective” frequency $\nu_e$ of a weighted wideband image as the frequency at which the image flux density equals the source flux density $\nu_e^\alpha$:

$$ \nu_e = \left( \frac{\sum_{i=1}^{14} w_i \nu_i^\alpha}{\sum_{i=1}^{14} w_i} \right)^{1/\alpha}. \qquad (17) $$

Thus, different weighting schemes yield slightly different effective frequencies for the wideband images. The effective frequencies of our DEEP2 images are $\nu_e \approx 1329$ MHz for minimum $\sigma_n^2$ and $\nu_e \approx 1278$ MHz for maximum S/N weighting.

The rms noise in the S/N-weighted DEEP2 image was estimated in five widely separated and apparently source-free regions covering 14,373 CLEAN beam solid angles at offsets $\rho \sim 1.6''$ from the pointing center, where the primary attenuation is $a_b < 0.01$. The distribution of their peak flux densities $S_p$ is nearly Gaussian (Figure 8) with rms $\sigma_n = 0.55 \pm 0.01$ µJy beam$^{-1}$. The noise in a synthesis image not corrected for primary-beam attenuation should be uniform across the whole image.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Circularized effective primary attenuation pattern for the wideband DEEP2 image made with subband weights that maximize the S/N for sources with $\alpha = -0.7$ (solid curve) has slightly broader wings than the narrowband attenuation pattern at $\nu_e \approx 1278$ MHz (dotted curve).</figcaption></figure>

In each narrow subband the primary-beam attenuation pattern $a_{b,i}(\rho)$ is well defined and was measured accurately in the horizontal and vertical planes (Section 2.2.2) of the alt-az mounted MeerKAT dishes. The primary beam is slightly elliptical and rotates with parallactic angle on the sky. For the long DEEP2 tracks we approximated the ellipse by a circle whose diameter is the geometric mean of the horizontal and vertical diameters. The primary beamwidth is inversely proportional to frequency, so the effective primary pattern $a_b(\rho)$ of the weighted wideband image must be calculated from

$$a_b(\rho) = \frac{\sum_{i=1}^{14} w_i a_{b,i}(\rho) \nu_i^{\alpha}}{\sum_{i=1}^{14} w_i}. \quad (18)$$

The circularized attenuation pattern for the wideband DEEP2 image weighted for maximum S/N is shown in Figure 9. Its FWHM is $\theta_b \approx 68'$. Numerically it can be approximated within 0.1% for all $a_b > 0.25$ by the polynomial

$$a_b(\rho) \approx 1.0 - 0.3514x/10^3 + 0.5600x^2/10^7 - 0.0474x^3/10^{10} + 0.00078x^4/10^{13} + 0.00019x^5/10^{16}, \quad (19)$$

where $x \equiv [\rho(\text{arcmin})\nu_e(\text{GHz})]^2$ and $\nu_e \approx 1.278$ GHz. This primary-beam attenuation can be used by the AIPS task PBCOR via the parameter PBPARM = 0.250, 1.0, -0.3514, 0.5600, -0.0474, 0.00078, and 0.00019. The wideband attenuation pattern is close to the narrowband attenuation pattern at $\nu = \nu_e \approx 1278$ MHz shown by the dotted curve in Figure 9, but it has slightly broader wings.

The wideband DEEP2 images are so sensitive ($\sigma_n \approx 0.55$ µJy beam$^{-1}$) that they are densely covered by sources with flux densities $S \gg \sigma_n$ (Figures 10 and 11). The total flux density in a dirty interferometer image is always zero because there are no zero-spacing $(u, v)$ data, only sinusoidal fringes with zero means. The dirty image of a single strong point source contains a negative “bowl” whose angular size is inversely proportional to the smallest $(u^2 + v^2)^{1/2}$ sampled. For our DEEP2 data, the bowl surrounding each source is much wider than the source spacing but much narrower than the primary attenuation pattern $a_b(\rho)$ (Figure 9). Faint sources have a fairly uniform random distribution on the sky, and their attenuated brightness distribution in the DEEP2 images is

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Central $5' \times 5'$ of the wideband DEEP2 image made with subband weights that maximize the S/N for sources with $\alpha = -0.7$ is covered by sources unresolved by the $\theta = 7''$ beam and brighter than the rms noise $\sigma_n \approx 0.55$ µJy beam$^{-1}$. The $-1.4$ µJy beam$^{-1}$ “bowl” (as described in Section 4.5) has been removed from this image.</figcaption></figure>

multiplied by the primary attenuation pattern. Consequently, each dirty DEEP2 image has a negative bowl whose size and shape closely match the primary beam and whose depth exceeds $\sigma_n$. Partial CLEANing reduces the depth of the bowl but does not completely eliminate it.

We estimated the central depth of the bowl of our S/N-optimized DEEP2 image using the mode of the brightness distribution in the $6' \times 6'$ square at the center of the wideband CLEAN image, where the primary attenuation is confined to the narrow range $0.98 < a_b < 1.00$ (Figure 9); it is $-1.4 \pm 0.1$ µJy beam$^{-1}$. To fill in the bowl and flatten the image baseline level, we added the weighted wideband primary attenuation pattern $a_b(\rho)$ (Equation (19)) multiplied by $1.4 \times 10^{-6}$ to the wideband image, whose brightness units are Jy beam$^{-1}$.

The central square cutout 4640 pixels $\times$ 1.25 pixel$^{-1} = 5800''$ on a side from the S/N-weighted 1.28 GHz wideband DEEP2 image with the bowl removed, but not corrected for primary-beam attenuation, is available in FITS format at https://archive.sarao.ac.za.

## 5. Source Counts at 1.4 GHz

We used the DEEP2 confusion $P(D)$ distribution to make the best power-law approximation to the sky density of sources fainter than 10 µJy, and we counted individual DEEP2 sources with $-5.0 < \log[S(\text{Jy})] < -2.6$.

### 5.1. Deep Power-law P(D) Counts

The differential number $n(S)$ of sources per unit flux density per steradian is a statistical quantity that can be calculated directly from the distribution of image brightnesses expressed in units of flux density per beam solid angle. The term “confusion” describes the brightness fluctuations caused by radio sources, and an image is said to be confusion limited if the confusion fluctuations are larger than the rms noise
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Central 1° × 1° of the wideband DEEP2 sky image made with subband weights that maximize the S/N for sources with $\alpha = -0.7$. The $-1.4~\mu\text{Jy beam}^{-1}$ bowl described in Section 4.5 has been removed from this image, and it has also been corrected for the primary-beam attenuation using Equation (19). The gray scale is stretched by an exponent of 1.3 between -15 and 30 $\mu\text{Jy}$ as indicated by the bar at the top. The dashed square in the center of the image bounds the $1250'' \times 1250''$ region whose $P(D)$ distribution we used to calculate the power-law source count described in Section 5.1.</figcaption></figure>

fluctuations $\sigma_n$. The normalized probability distribution of confusion brightness is traditionally called the $P(D)$ distribution, where $D$ originally stood for the pen deflection on a chart recording but now refers to image peak flux density $S_p$, that is, flux density per beam solid angle.

The noiseless $P(D)$ distribution for point sources can be derived analytically only for power-law differential source counts $n(S) = kS^{-\gamma}$, where $k$ is the overall source density parameter and $\gamma$ is the power-law exponent (Condon 1974). Power-law distributions are scale-free, so the shape of a power-law $P(D)$ distribution depends only on $\gamma$. The amplitudes and widths of power-law $P(D)$ distributions obey the scaling relation

$$P[(k\Omega_e)^{1/(\gamma-1)}D] = (k\Omega_e)^{-1/(\gamma-1)}P(D), \quad (20)$$

where

$$\Omega_e \equiv \int [a(\rho, \phi)]^{\gamma-1} d\Omega \quad (21)$$

is the “effective” solid angle $\Omega_e$ of a beam whose polar attenuation pattern is $a(\rho, \phi)$. For power-law counts, the $P(D)$ distribution depends on $\Omega_e$ but not on the form of the beam attenuation pattern. The effective solid angle of an elliptical

Gaussian beam with FWHM axes $\theta_1$ and $\theta_2$ is

$$\Omega_e = \left( \frac{\pi\theta_1\theta_2}{4\ln 2} \right) \left( \frac{1}{\gamma - 1} \right) = \frac{\Omega_b}{\gamma - 1}, \quad (22)$$

where $\Omega_b$ is the Gaussian beam solid angle

$$\Omega_b \equiv \int a(\rho, \phi) d\Omega = \frac{\pi\theta_1\theta_2}{4\ln 2}. \quad (23)$$

A convenient normalization for displaying analytic $P(D)$ distributions is $k\Omega_e = \eta_1^{-1}$, where

$$\eta_1^{-1} = \frac{2\Gamma(\gamma/2)\Gamma[(\gamma+1)/2]\sin[\pi(\gamma-1)/2]}{\pi^{\gamma+1/2}} \quad (24)$$

and $\Gamma(x)$ is the factorial function (Condon 1974).

The shape of the $P(D)$ distribution varies significantly with the count slope $\gamma$. Four examples of $P(D)$ distributions with $k\Omega_e = \eta_1^{-1}$ are shown in Figure 12. The super-Euclidean slope $\gamma = 2.8$ is close to the actual slope observed above $S \sim 1$ Jy at 1.4 GHz, so the $\gamma = 2.8$ curve in the top panel of Figure 12 represents the troublesome confusion that affected early radio surveys such as 2C (Shakeshaft et al. 1955). In the limit $\gamma \to 3-$, the $P(D)$ distribution is Gaussian and would appear as a parabola in the semilogarithmic Figure 12. The $\gamma = 2.8$ curve is dominated by its nearly parabolic core and has only a weak tail extending to the right. The sky brightness contributed by point sources diverges for all $\gamma \ge 2$ (Olbers’s paradox), so the
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Noiseless *P(D)* distributions for power-law source counts *n(S)* ~ S⁻⁷. The top panel shows γ > 2 counts for which the sky brightness diverges (Olbers's paradox), so the deflections are shown relative to the mean deflection ⟨D⟩. The bottom panel shows *P(D)* distributions for γ < 2 relative to the absolute zero of the source sky.</figcaption></figure>

calculated values of *D* are relative to the mean deflection ⟨*D*⟩. The γ = 2.1 curve shows how the *P(D)* peak shifts to the left and the tail becomes more prominent as γ → 2+.

The bottom panel of Figure 12 shows sample *P(D)* distributions when γ < 2, and *D* represents the deflection above the absolute zero of sky brightness contributed by radio sources. The γ = 1.9 curve is similar to the γ = 2.1 curve in form, and its peak *D* is still significantly offset above zero by contributions from faint sources so numerous that multiple sources are blended together in each beam. This is characteristic of confusion seen at levels 10 µJy ~ S₁.₄ GHz ~ 0.1 Jy.

At the lower value γ ≈ 1.5 observed when S₁.₄ GHz ≪ 10 µJy, the nature of confusion changes again. There are so few sub-µJy sources that the *P(D)* peak deflection approaches *D* → 0+, indicating that the extragalactic background has been largely resolved into discrete sources. The long tail of the *P(D)* distribution so completely dominates the narrow core that the rms confusion σ_c is ill-defined and should not be used to describe the amount of confusion when S₁.₄ GHz ≪ 10 µJy. Confusion in the traditional sense “melts away” at the sub-µJy levels reached by DEEP2. Instead of numerous even fainter sources tending to boost the flux densities of faint sources, faint sources are more likely to suffer obscuration by stronger sources.

The central 1250'' × 1250'' square covering about 2.4 × 10⁴ restoring beam solid angles was extracted from the DEEP2 sky image (Figure 11) and rescaled to 1.4 GHz via the spectral index α = −0.7 typical of faint sources. Its 1.4 GHz *P(D)* distribution is shown by the red curve in Figure 13. The best least-squares power-law fit with fixed rms noise σ_n = 0.55 ± 0.01 µJy beam⁻¹ but free *k* and γ is *n(S)* = 1.07 × 10⁵ S⁻¹.⁵² Jy⁻¹ sr⁻¹ between *S* = 0.25 and

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Solid black curve is the noiseless *P(D)* distribution rescaled for a θ = 7''6 Gaussian beam and 1.4 GHz source count *n(S)* = kS⁻⁷. The dotted parabola is the normalized probability distribution of the σ_n = 0.55 µJy beam⁻¹ Gaussian noise, and the blue curve is the convolution of the noiseless calculated *P(D)* distribution with the noise. The observed *P(D)* distribution in the central 1250'' × 1250'' (~2.4 × 10⁴ Ω_b) of DEEP2 is shown as the red curve.</figcaption></figure>

10 µJy. Letting σ_n be a free parameter and varying σ_n by 2Δσ_n = ±0.02 µJy beam⁻¹ has a negligible effect on this fit. The fit and its rms uncertainties (68% confidence region) are bounded by the thick box at the left in Figure 14. The actual source count is not a perfect power law, and the DEEP2 dirty beam is not perfectly Gaussian, so getting more accurate source counts will require numerical simulations based on more realistic non-power-law source counts and non-Gaussian dirty beams (A. M. Matthews et al. 2019, in preparation).

Smoothly extrapolating our power-law 1.4 GHz count below S₀ ≈ 0.25 µJy is reasonable because fainter sources are evolved from the power-law region below the “knee” in the local RLF (Condon et al. 2019), so simple evolutionary models (Condon 1984) predict nearly power-law sub-µJy source counts. This extrapolation yields an estimate of the Rayleigh–Jeans brightness temperature T_b contributed by all fainter SFGs:

$$ \Delta T_b(<S_0) = \left[ \frac{\ln(10)c^2}{2k_B\nu^2} \right] \int_0^{S_0} S^2 n(S) d[\log S], \quad (25) $$

where k_B ≈ 1.38 × 10⁻²³ J K⁻¹ is the Boltzmann constant (Condon et al. 2012). For *n(S)* = 1.07 × 10⁵ S⁻¹.⁵² Jy⁻¹ sr⁻¹ at 1.4 GHz, the background contributed by sources fainter than S₀ = 0.25 µJy is ΔT_b ≈ 2.5 mK, which is only ≈7% of the T_b ≈ 37 mK total background (Condon et al. 2012) contributed by SFGs and <3% of the background contributed by all extragalactic sources.

## 5.2. Confusion and Obscuration Sensitivity Limits for Individual Sources

An image is said to be confusion limited if the errors caused by confusion exceed the errors caused by Gaussian noise. Both confusion and noise set lower limits to the flux density *S* of the faintest individual source that can be reliably detected. This section extends earlier calculations of the confusion limit below
<figure><img src="image_16.png" /><figcaption>**Figure 14.** The 1.4 GHz differential source count has been plotted with the traditional static Euclidean normalization $S^{5/2}n(S)$ (top panel) and the brightness-weighted normalization $S^2n(S)$ (bottom panel). The black data points show the DEEP2 discrete source counts from Table 6, the red data points are from Prandoni et al. (2018), the green data points are from Smolčić et al. (2017), and the blue data points are from Hopkins et al. (2003). The box covering $-5.0 < \log[S(Jy)] < -4.1$ bounds the Mitchell & Condon (1985) 1.4 GHz $P(D)$ count, and the box covering $-5.8 < \log[S(Jy)] < -4.8$ bounds the 3 GHz $P(D)$ count from Condon et al. (2012) converted to 1.4 GHz via spectral index $\alpha = -0.7$. The heavy box spanning $-6.6 < \log[S(Jy)] < -5$ indicates the best fit to the DEEP2 confusion $P(D)$ distribution discussed in Section 5.1.</figcaption></figure>

$S \sim 10$ µJy, where $\gamma \ll 2$ and obscuration dominates confusion. This flux density range affects the confusion/obscuration corrections to DEEP2 direct source counts (Section 5.3) and directly impacts the design of future arrays such as the SKA and ngVLA.

A common way to compare confusion and noise is to calculate their variances $\sigma_c^2$ and $\sigma_n^2$. This calculation must be done carefully because the confusion variance $\sigma_c^2 = \int_0^\infty D^2P(D) dD$ formally diverges for all power-law source counts $n(S) = k S^{-\gamma}$ and is finite only if the $P(D)$ distribution is truncated above some cutoff deflection $D_c$. Then (Condon 1974),

$$ \sigma_c = \left( \frac{k\Omega_e}{3-\gamma} \right)^{1/2} D_c^{(3-\gamma)/2}, \quad 1 < \gamma < 3. \qquad (26) $$

The cutoff should be proportional to $\sigma_c$: $D_c = q\sigma_c$, where the constant $q \sim 5$ is the usual cutoff signal-to-confusion ratio. Then,

$$ \sigma_c = \left( \frac{k\Omega_e}{3-\gamma} \right)^{1/(\gamma-1)} q^{(3-\gamma)/(\gamma-1)}. \qquad (27) $$

Note that the rms confusion $\sigma_c$ still depends on the choice of $q$. Only as $\gamma \to 3-$ is $\sigma_c$ nearly independent of $q$. When $\gamma = 2$, $\sigma_c$

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Solid curve is the minimum number $\beta_{min}$ of beam solid angles per source for a confusion $S/N q = 5$ (Equation (31)), and the dotted curve shows an alternative $\beta_{min}$ for which the probability of confusion is $P_c < 0.159$ (the probability that Gaussian noise exceeds $+\sigma_n$) for confusion $> S/5$ (Equation (34)). The dashed curve is the minimum number of sources per beam solid angle at which 10% of sources are obscured by stronger sources (Equation (36)).</figcaption></figure>

$\propto q$. In the sub-µJy regime where $\gamma \sim 1.5$, $\sigma_c \propto q^3$ depends so sensitively on $q$ that the very concept of rms confusion is nearly useless.

For imaging with a Gaussian beam, $\Omega_e = \Omega_b/(\gamma - 1)$ and

$$ \sigma_c^{\gamma-1} = \frac{k\Omega_b q^{3-\gamma}}{(\gamma-1)(3-\gamma)}. \qquad (28) $$

Solving the cumulative source count

$$ N(>S) = \frac{kS^{1-\gamma}}{\gamma-1} \qquad (29) $$

for $k = (\gamma - 1) N(>S)$ $S^{\gamma-1}$ and substituting the detection limit $S \approx q\sigma_c$ results in

$$ \left(\frac{\sigma_c}{S}\right)^{\gamma-1} \approx \left(\frac{q^{3-\gamma}}{3-\gamma}\right) N(>S) \Omega_b = \left(\frac{q^{3-\gamma}}{3-\gamma}\right)^{\beta-1}, \qquad (30) $$

which can be solved for $\beta_{min}$, the minimum number of beam solid angles per reliably detectable source, in terms of the confusion $S/N q$:

$$ \beta_{min} \approx \frac{q^2}{3-\gamma}. \qquad (31) $$

The number of beams per source at the $q = 5$ confusion limit is shown as a function of $\gamma$ by the solid curve in Figure 15. In the limit $\gamma \to 3-$, the $P(D)$ distribution is dominated by the very faintest sources and the *fluctuations* in sky brightness diverge, just as the total sky brightness diverges as $\gamma \to 2-$ (Olbers's paradox). The $P(D)$ distribution becomes nearly Gaussian, the same as the noise distribution, so sources can never be distinguished from noise and $\beta \to \infty$. For sources stronger than $S \sim 1$ Jy at 1.4 GHz, a super-Euclidean differential count slope $\gamma \approx 2.7$ implies $\beta_{min} \approx 80$ beam solid angles per reliable $q = 5$ source detection. This very severe requirement on $\beta_{min}$ was not obvious when the first extragalactic radio surveys were being made, so sources with smaller $\beta$ were often cataloged as real and later found to be spurious.
An alternative approach to calculating the confusion $\beta_{min}$ uses the probability $P_c$ that a source of flux density $S$ will be confused by a weaker source of flux density between $fS$ and $S$, where $f < 1$. For a cumulative source count $N(>S) \propto S^{1-\gamma}$ and a top-hat beam ($a(\theta, \phi) = 1$ over solid angle $\Omega_e$ in the notation of Equation (21)),

$$P_c = \Omega_e[N(>fS) - N(>S)] = \Omega_e[N(>S)(f^{1-\gamma} - 1)]. \quad (32)$$

A Gaussian beam with beam solid angle $\Omega_b = \Omega_e(\gamma - 1)$ yields exactly the same $P(D)$ distribution, so

$$P_c = N(>S)\Omega_b\left(\frac{f^{1-\gamma}-1}{\gamma-1}\right) = \frac{1}{\beta}\left(\frac{f^{1-\gamma}-1}{\gamma-1}\right) \quad (33)$$

and

$$\beta_{min} = \frac{1}{P_c}\left(\frac{f^{1-\gamma}-1}{\gamma-1}\right) \quad (34)$$

is the minimum number of Gaussian beams per source consistent with this confusion requirement. A reasonable choice for $f$ would be $f \approx 0.2$, so the confusion is at least $S/5$. A reasonable choice for the probability of confusion at this level is $P_c \approx 0.159$, the probability that Gaussian noise exceeds $+\sigma_n$. The dotted curve in Figure 15 shows the resulting $\beta_{min}$ as a function of $\gamma$.

In the broad flux density range $10 \mu Jy < S < 0.1$ Jy covered by most 1.4 GHz surveys, $\gamma \sim 2$ and $\beta_{min} \sim 25$. Below $S \sim 10 \mu Jy$, the differential count slope falls again, to $\gamma \lesssim 1.5$. While $\beta_{min}$ for avoiding confusion by fainter sources continues to fall, the probability of obscuration by stronger sources grows. For a top-hat beam the probability that a source of flux density $S$ will be obscured by a stronger source is $P_o = N(>S) \Omega_e$. For power-law source counts, the $P(D)$ distribution is independent of beam shape and depends only on $\Omega_e$, so for a Gaussian beam with beam solid angle $\Omega_b = (\gamma - 1) \Omega_e$ (Equation (22)),

$$P_o = N(>S) \Omega_e = \frac{N(>S) \Omega_b}{\gamma - 1}, \quad (35)$$

and the minimum number of beam solid angles per source to minimize obscuration is

$$\beta_{min} \approx [P_o(\gamma - 1)]^{-1}. \quad (36)$$

The dashed curve in Figure 15 shows the number $\beta_{min}$ of beam solid angles per source corresponding to $P_o = 0.1$.

Choosing $\beta_{min} \approx 25$ is a good rule of thumb for detecting individual sources valid in the flux density range below $S_{1.4 GHz} \sim 0.1$ Jy, where $1.4 \lesssim \gamma \lesssim 2$. The solid curve in Figure 16 shows the 1.4 GHz flux density $S$ at which $\beta = 25$. For example, the 1.4 GHz EMU survey (Norris et al. 2011) has $\theta = 10''$ FWHM resolution, rms confusion $S/\Omega_b \approx 240$ nJy arcsec$^{-2}$, and beam solid angle $\Omega_b = \pi\theta^2/(4\ln2) \approx 113$ arcsec$^2$. The weakest reliably detectable sources at that resolution have flux densities $S \approx 27$ $\mu$Jy.

The sharp decline of the minimum $S$ caused by the low $\gamma \sim 1.5$ at $\theta \ll 10''$ means that even the most sensitive SKA images will not be confusion limited at beamwidths $\theta \gtrsim 1''$. This is fortunate because the median angular size of faint SFGs is $\langle\phi\rangle = 0''3 \pm 0''1$ with an rms scatter $\sigma_{\phi} \lesssim 0''3$ (Cotton

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Solid curve shows the individual-source detection limit $S$ divided by the beam solid angle $\Omega_b$ in units of nJy arcsec$^{-2}$ (left ordinate) and Rayleigh-Jeans brightness temperature $\sigma_T$ in units of mK (right ordinate) at 1.4 GHz at which $\beta = 25$. The dotted line is $5\sigma_c$ calculated for the Loi et al. (2019) 1.4 GHz SKA sky simulation.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 17.** Solid curve based on Equation (37) shows the flux density $\langle S \rangle$ of the strongest source in half of the Gaussian beam areas of an image made with resolution $\theta$.</figcaption></figure>

et al. 2018), and sources with $\phi \gtrsim \theta$ are more difficult to detect because their peak flux densities are reduced by factors $\gtrsim 2$. Although the Loi et al. (2019) 1.4 GHz simulation of the radio sky for the SKA and its precursors also predicts $\gamma \lesssim 1.5$ below $S = 10 \mu Jy$, their estimated rms confusion $\sigma_{mJy/bm} = (0.237 \pm 0.001)(\nu_{GHz})^{-0.8}(\theta)^{2.149 \pm 0.001}$ (the dotted line in Figure 16 shows their $5\sigma_{mJy/bm}$) does not take the lower $\gamma$ into account and overpredicts the rms confusion at nJy flux densities.

Statistical counts using the $P(D)$ distribution can usefully reach much lower $\beta$ values and hence much lower flux densities. For point sources randomly placed on the sky, the Poisson probability that any beam solid angle will contain no sources stronger than $S$ is $P_P = \exp(-1/\beta)$. Half of all beam solid angles ($P_P = 0.5$) must satisfy $\beta = 1/\ln(2) \approx 1.44$, and the flux density of the strongest source in them, $\langle S \rangle$, is the solution of

$$N(>\langle S \rangle) = \frac{\ln 2}{\Omega_b} = \frac{1}{\pi} \left( \frac{2 \ln 2}{\theta} \right)^2. \quad (37)$$

For the $\theta = 7''6$ FWHM DEEP2 beam and Condon (1984) model 1.4 GHz source count, Figure 17 shows that half of all
<table><caption>**Table 6**
DEEP2 1.278 GHz Source Counts</caption><thead><tr><th>⟨log[S(Jy)]⟩</th><th>n<sub>bin</sub></th><th>log[S<sup>2</sup>n(S)(Jy sr<sup>-1</sup>)]</th></tr></thead><tbody><tr><td>-4.9</td><td>5572</td><td>2.682 ± 0.031</td></tr><tr><td>-4.7</td><td>4521</td><td>2.754 ± 0.031</td></tr><tr><td>-4.5</td><td>3025</td><td>2.762 ± 0.032</td></tr><tr><td>-4.3</td><td>2008</td><td>2.770 ± 0.032</td></tr><tr><td>-4.1</td><td>1070</td><td>2.740 ± 0.033</td></tr><tr><td>-3.9</td><td>576</td><td>2.690 ± 0.036</td></tr><tr><td>-3.7</td><td>289</td><td>2.578 ± 0.040</td></tr><tr><td>-3.5</td><td>139</td><td>2.444 ± 0.047</td></tr><tr><td>-3.3</td><td>85</td><td>2.450 ± 0.055</td></tr><tr><td>-3.1</td><td>49</td><td>2.428 ± 0.066</td></tr><tr><td>-2.9</td><td>30</td><td>2.411 ± 0.080</td></tr><tr><td>-2.7</td><td>24</td><td>2.532 ± 0.087</td></tr></tbody></table>

beam solid angles contain no sources stronger than ⟨S⟩ ≈ 0.25 µJy. Although the observed DEEP2 P(D) distribution (Figure 13) is broadened by σn ≈ 0.55 µJy beam⁻¹ noise, it samples ~1.2 × 10⁴ independent effective beam areas Ωe, so it should be able to constrain the 1.4 GHz source count down to the S ≈ 0.25 µJy SKA1 goal for studying the SFHU (Prandoni & Seymour 2015). At 1.4 GHz the ratio of the S ~ 16 µJy detection limit βmin = 25 for discrete sources to the confusion sensitivity limit ⟨S⟩ ≈ 0.25 µJy is about 64! Another advantage of P(D) counts is that the low minimum β ~ 1.44 allows relatively large beamwidths, while the larger minimum β ≈ 25 for deep counts of individual sources requires beamwidths small enough that large and difficult corrections for partial source resolution become necessary (Owen 2018).

### 5.3. Direct Counts of DEEP2 Sources

Figure 16 indicates that the β = 25 limit for reliably detecting individual sources with the DEEP2 θ₁/₂ = 7.6' beam is S ≈ 17 µJy at 1.278 GHz (≈16 µJy at 1.4 GHz). Below that limit, a significant fraction of sources will be confused or obscured. We produced preliminary counts of DEEP2 sources stronger than 10 µJy inside the Ω ≈ 1.026 deg² DEEP2 half-power circle. To estimate confusion and obscuration corrections, we simulated a DEEP2 image, counted sources in the simulated image, and compared those counts with the actual counts used in the simulation.

For the simulation, point sources with approximately the correct source count were randomly placed on the image and convolved with the DEEP2 dirty beam. Then, the simulated image was multiplied by the primary attenuation, convolved with 0.55 µJy beam⁻¹ rms noise, and CLEANed. A source catalog was extracted from the simulated image, and the source counts from this catalog were compared with the input counts to derive count corrections. These corrections are typically ~10%, and we estimate that their rms uncertainties are about half of the corrections themselves. The brightness-weighted 1.278 GHz DEEP2 source counts log[S²n(S)(Jy sr⁻¹)] in 12 bins of logarithmic width Δ = 0.2 in log(S) centered on ⟨log[S(Jy)]⟩ = -4.9, -4.7, ..., -2.7 and their rms errors are listed in Table 6.

We converted the 1.278 GHz counts to 1.4 GHz via a median spectral index -0.7, and these 1.4 GHz counts are plotted as the black points in Figure 14. The top panel of Figure 14 presents the traditional static Euclidean normalization S⁵/²n(S), and the bottom panel shows the brightness-

weighted normalization S²n(S). The red points are from Prandoni et al. (2018), the blue points are from Hopkins et al. (2003), and the green points are from Smolčić et al. (2017). The DEEP2, Prandoni et al. (2018), and Hopkins et al. (2003) counts agree within the errors. The counts from Smolčić et al. (2017) are somewhat lower in the range -4 ≲ log[S(Jy)] ≲ -3, perhaps because some extended sources were partially resolved by their 0.75 beam. Figure 14 also shows that the DEEP2 direct count is slightly higher than the Mitchell & Condon (1985) 68% confidence P(D) box, possibly because DEEP2 is significantly more sensitive (rms noise σn = 0.55 µJy beam⁻¹ ≈ 7.2 mK) to low-brightness SFGs.

## 6. DEEP2 and the SFHU

The 1.4 GHz continuum emission from an SFG is a combination of synchrotron radiation from relativistic electrons accelerated in the supernova remnants of short-lived (τ ≲ 30 Myr) massive stars and free-free radiation from thermal electrons in H II regions ionized by stars that are even more massive and short-lived. It is uniquely unbiased by dust or older stars. The tight and fairly linear FIR/radio correlation (Condon et al. 1991) indicates that the radio luminosity of an SFG is directly proportional to its recent SFR. Thus,

$$ \frac{\text{SFR}(M > 0.1 M_{\odot})}{M_{\odot} \text{ yr}^{-1}} = \kappa \left( \frac{L_{1.4 \text{ GHz}}}{\text{W Hz}^{-1}} \right), \quad (38) $$

where the constant of proportionality κ is in the range of (0.6–1.2) × 10⁻²¹ (Condon 1992; Sullivan et al. 2001; Bell 2003; Kennicutt et al. 2009; Murphy et al. 2011) and depends on the unknown initial mass function below M ~ 8 M☉. However, the only consequence of the uncertainty in κ is a global scaling in the total mass of stars in the universe. It does not alter the determination of the evolutionary models of the SFRD through comparisons of local RLFs with radio source counts. Therefore, the local RLF and the DEEP2 confusion probability distribution will constrain the luminosity and density evolution of SFGs independent of dust, older stars, and SFR calibration errors.

The universe is homogeneous on large scales, so its spatially averaged SFRD ψ depends only on cosmic time or a proxy for time such as redshift z. It is usually written in units of M☉ yr⁻¹ Mpc⁻³. Combining the 1.4 GHz local luminosity function of SFGs and the source count of distant SFGs yields an independent extinction-free means of measuring the SFHU. The main obstacle has been the fact that SFGs are intrinsically weak radio sources, so tracing the formation of most stars, and not just the tip of the iceberg in ultraluminous starburst galaxies, requires counting sources fainter than S ~ 1 µJy.

Our method for calculating the SFRD in the standard ΛCDM universe with Ωm = 0.3 and H₀ = 70 km s⁻¹ Mpc⁻¹ from radio data follows Appendix C of Condon & Matthews (2018). Let ρ(L_ν|z)dL be the number of sources with spectral luminosities L_ν to L_ν + dL_ν in comoving volume element dVc and let η(S, z)dS dz be the number of sources per steradian with flux densities S to S + dS in the redshift range z to z + dz. Then,

$$ \eta(S, z) dS dz = \rho(L_{\nu}|z) dL_{\nu} dV_{c}, \quad (39) $$
where the comoving volume element per steradian is

$$dV_C = \frac{D_C^2 D_{H_0}}{E(z)} dz, \qquad (40)$$

$D_{H_0} \equiv c/H_0$, and $E(z) \equiv H/H_0$ specifies the expansion history of the universe. For sources with spectral index $\alpha$,

$$L_{\nu} = 4\pi D_{C}^{2}(1 + z)^{1-\alpha} S. \qquad (41)$$

Multiplying both sides of Equation (39) by $S^2$ converts this redshift-dependent differential source count into a brightness-weighted source count

$$S^2\eta(S, z) = L_{\nu}^2 \rho(L_{\nu}|z) \left[ \frac{(1+z)^{\alpha-1} D_{H_0}}{4\pi E(z)} \right]. \qquad (42)$$

Similarly, multiplying the spectral luminosity function $\rho(L_\nu|z)$ by luminosity emphasizes the luminosity ranges contributing the most to the spectral luminosity density $U(L_\nu|z) \equiv L_\nu \rho(L_\nu|z)$. SFGs span several decades of luminosity, so it is convenient to replace $U$ by the spectral luminosity density per decade of luminosity,

$$U_{\text{dex}}(L_{\nu}|z) = L\rho(L_{\nu}|z) \frac{dL}{d\log L} = L_{\nu}^2 \rho(L_{\nu}|z) \ln(10). \qquad (43)$$

The local energy density functions $U_{\text{dex}}(L_{1.4\,\text{GHz}}|z=0)$ for SFGs and AGN-powered radio galaxies were separately derived from a spectroscopically complete sample containing 9517 NVSS sources cross-identified with Two Micron All Sky Survey (2MASS) galaxies (Condon et al. 2019). In terms of $U_{\text{dex}}$,

$$S^2\eta(S, z) = U_{\text{dex}}(L_{\nu}|z) \left[ \frac{(1+z)^{\alpha-1} D_{H_0}}{4\pi \ln(10) E(z)} \right]. \qquad (44)$$

The brightness-weighted source count is obtained by integrating over redshift:

$$S^2 n(S) = \frac{D_{H_0}}{4\pi \ln(10)} \int_0^\infty U_{\text{dex}}(L_\nu|z) \left[ \frac{(1+z)^{\alpha-1}}{E(z)} \right] dz. \qquad (45)$$

Thus, the evolving SFRD $\psi(z)$ can be constrained by comparing the observed brightness-weighted 1.4 GHz source count $S^2n(S)$ with counts predicted by evolving the local energy density function $U_{\text{dex}}(L_{1.4\,\text{GHz}}|z=0)$ with redshift $z$ using various evolutionary models (e.g., Hopkins & Beacom 2006; Madau & Dickinson 2014). For example, Madau & Dickinson (2014) derived the SFRD evolutionary model

$$\frac{\psi(z)}{\psi(0)} = \frac{(1+z)^{2.7}}{1 + [(1+z)/2.9]^{5.6}} \qquad (46)$$

by fitting available UV and infrared data.

Pure luminosity evolution consistent with Equation (46) for SFGs predicts the brightness-weighted source counts shown in Figure 18 by the blue curve. The red curve is an estimate of the AGN contribution (Condon 1984). The black curve is their sum, and it is significantly lower than the actual 1.4 GHz source count in the range $-5 \lesssim \log[S(\text{Jy})] \lesssim -4$. The DEEP2 source count extends $\sim 7 \times$ deeper than the Condon et al. (2012) source count (Figure 18), constraining for the first time star formation in all galaxies with SFRs as low as $5 M_{\odot} \text{yr}^{-1}$ at “cosmic noon” ($z \approx 2$), not just the rare starbursts with SFR $> 35 M_{\odot} \text{yr}^{-1}$. DEEP2 also reduces the statistical

<figure><img src="image_20.png" /><figcaption>Figure 18. Blue and red curves represent the brightness-weighted counts $S^2n(S)$ of SFGs (from Madau & Dickinson 2014) and AGNs (from Condon 1984), respectively. Their sum is shown by the black curve, which lies significantly below the observed counts from Figure 14 in the range $-5 \lesssim \log[S(\text{Jy})] \lesssim -4$.</figcaption></figure>

uncertainty in $S^2n(S)$ by a factor of four because it samples a much larger solid angle of sky.

While this simple preliminary analysis suggests that SFGs evolve more strongly, it uses a power-law approximation for the faint-source $P(D)$ counts and assumes a perfectly Gaussian beam. To fully exploit our DEEP2 image, we are now developing numerical simulations of synthetic images populated with sources having arbitrary counts and beams that take the limitations of CLEAN into account (A. M. Matthews et al. 2019, in preparation).

The MeerKAT telescope is operated by the South African Radio Astronomy Observatory, which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated by Associated Universities, Inc. This material is based on work supported by the National Science Foundation Graduate Research Fellowship under grant No. DDGE-1315231. We thank our anonymous referee for a detailed and constructive review.

## ORCID iDs

T. Mauch https://orcid.org/0000-0003-2716-9589

J. J. Condon https://orcid.org/0000-0003-4724-1939

A. M. Matthews https://orcid.org/0000-0002-6479-6242

F. Camilo https://orcid.org/0000-0002-1873-3718

## References

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

Bell, E. F. 2003, ApJ, 586, 794

Best, P. N., Kauffmann, G., Heckman, T. M., & Ivezić, Ž. 2005, MNRAS, 362, 9

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

Camilo, F., Scholz, P., Serylak, M., et al. 2018, ApJ, 856, 180

Colless, M., Dalton, G., Maddox, S., et al. 2001, MNRAS, 328, 1039

Condon, J. J. 1974, ApJ, 188, 279

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

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

Condon, J. J. 1997, PASP, 109, 166

Condon, J. J. 2009, Sensitive Continuum Surveys with the SKA: Goals and Challenges Memo 114, (Cheshire: SKA), https://www.skatelescope.org/uploaded/12336_114_Memo_Condon.pdf

Condon, J. J., Anderson, M. L., & Helou, G. 1991, ApJ, 376, 95

Condon, J. J., & Broderick, J. J. 1988, AJ, 96, 30

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

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

Condon, J. J., & Matthews, A. M. 2018, PASP, 130, 073001

Condon, J. J., Matthews, A. M., & Broderick, J. J. 2019, ApJ, 872, 148
Condon, J. J., & Ransom, S. M. 2016, Essential Radio Astronomy (Princeton, NJ: Princeton Univ.)

Cotton, W. D. 2008, PASP, 120, 439

Cotton, W. D., Condon, J. J., Kellermann, K. I., et al. 2018, ApJ, 856, 67

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

Hopkins, A. M., & Beacom, J. F. 2006, ApJ, 651, 142

Jarvis, M., Seymour, N., Afonso, J., et al. 2015, in Proc. Advancing Astrophysics with the Square Kilometre Array (AASKA14) Conf., ed. T. L. Bourke et al. (Trieste: PoS), 68

Jarvis, M., Taylor, R., Agudo, I., et al. 2016, in Proc. MeerKAT Science: On the Pathway to the SKA Conf., ed. R. Taylor et al. (Trieste: PoS), 51

Jonas, J. & MeerKAT Team 2016, in Proc. MeerKAT Science: On the Pathway to the SKA Conf., ed. R. Taylor (Trieste: PoS), 1

Jones, D. H., Read, M. A., Saunders, W., et al. 2009, MNRAS, 399, 683

Kennicutt, R. C., Jr., Hao, C.-N., Calzetti, D., et al. 2009, ApJ, 703, 1672

Lehmensiek, R., & Theron, I. P. 2012, in Proc. Int. Conf. Electromagn. Adv. Appl. (ICEAA) Conf., ed. Zich R. S. et al. (Piscataway, NJ: IEEE), 321

Lehmensiek, R., & Theron, I. P. 2014, in Proc. 31st URSI General Assembly and Scientific Symp. (GASS) (Piscataway, NJ: IEEE), 1

Loi, F., Murgia, M., Govoni, F., et al. 2019, MNRAS, 485, 5285

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Mauch, T., Murphy, T., Buttery, H. J., et al. 2003, MNRAS, 342, 1117

Mauch, T., & Sadler, E. M. 2007, MNRAS, 375, 931

McAlpine, K., Jarvis, M. J., & Bonfield, D. G. 2013, MNRAS, 436, 1084

Mitchell, K. J., & Condon, J. J. 1985, AJ, 90, 1957

Murphy, E. J., Condon, J. J., Schinnerer, E., et al. 2011, ApJ, 737, 67

Norris, R. P., Hopkins, A. M., Afonso, J., et al. 2011, PASA, 28, 215

Offringa, A. R., de Bruyn, A. G., Biehl, M., et al. 2010, MNRAS, 405, 155

Owen, F. N. 2018, ApJS, 235, 34

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

Padovani, P., Miller, N., Kellermann, K. I., et al. 2011, ApJ, 740, 20

Prandoni, I., Guglielmino, G., Morganti, R., et al. 2018, MNRAS, 481, 4548

Prandoni, I., & Seymour, N. 2015, in Proc. Advancing Astrophysics with the Square Kilometre Array (AASKA14) Conf., ed. T. L. Bourke et al. (Trieste: PoS), 67

Reynolds, J. E. 1994, A Revised Flux Scale for the AT Compact Array, ATNF Memo. AT/39.3/040, (Epping: ATNF)

Sadler, E. M., Jackson, C. A., Cannon, R. D., et al. 2002, MNRAS, 329, 227

Shakeshaft, J. R., Ryle, M., Baldwin, J. E., Elsmore, B., & Thomson, J. H. 1955, MmRAS, 67, 106

Simpson, C., Rawlings, S., & Martínez-Sansigre, A. 2006, AN, 327, 270

Smolčić, V., Novak, M., Bondi, M., et al. 2017, A&A, 602, A1

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

Sullivan, M., Mobasher, B., Chan, B., et al. 2001, ApJ, 558, 72

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

York, D. G., Adelman, J., Anderson, J. E., Jr., et al. 2000, AJ, 120, 1579