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

# CHILES VERDES: Radio Variability at an Unprecedented Depth and Cadence in the COSMOS Field

Sumit K. Sarbadhicary<fnref n="1" />, Evangelia Tremou<fnref n="2" />, Adam J. Stewart<fnref n="3" />, Laura Chomiuk<fnref n="1" />, Charee Peters<fnref n="4" />, Chris Hales<fnref n="5" />, Jay Strader<fnref n="1" />, Emmanuel Momjian<fnref n="5" />, Rob Fender<fnref n="6" />,<fnref n="7" />, and Eric M. Wilcots<fnref n="4" />

Received 2020 September 12; revised 2021 April 19; accepted 2021 May 21; published 2021 December 8

## Abstract

Although it is well established that some extragalactic radio sources are time-variable, the properties of this radio variability, and its connection with host galaxy properties, remain to be explored—particularly for faint sources. Here we present an analysis of radio variable sources from the CHILES Variable and Explosive Radio Dynamic Evolution Survey (CHILES VERDES)—a partner project of the 1.4 GHz COSMOS H I Large Extragalactic Survey. CHILES VERDES provides an unprecedented combination of survey depth, duration, and cadence, with 960 hr of 1–2 GHz continuum VLA data obtained over 209 epochs between 2013 and 2019 in a 0.44 deg<fnref n="2" /> section of the well-studied extragalactic deep field, COSMOS. We identified 18 moderate-variability sources (showing 10%–30% flux density variation) and 40 lower-variability sources (2%–10% flux density variation). They are mainly active galactic nuclei (AGNs) with radio luminosities in the range of $10^{22}–10^{27}$ W Hz⁻<fnref n="1" />, based on cross-matching with COSMOS multiwavelength catalogs. The moderate-variability sources span redshifts $z = 0.22–1.56$, have mostly flat radio spectra ($\alpha > −0.5$), and vary on timescales ranging from days to years. The lower-variability sources have similar properties, but generally have higher radio luminosities than the moderate-variability sources, extending to $z = 2.8$, and have steeper radio spectra ($\alpha < −0.5$). No star-forming galaxy showed statistically significant variability in our analysis. The observed variability likely originates from scintillation on short (~week) timescales, and Doppler-boosted intrinsic AGN variability on long (month–year) timescales.

*Unified Astronomy Thesaurus concepts: Radio continuum emission (1340); Radio active galactic nuclei (2134); Time domain astronomy (2109); Time series analysis (1916)*

*Supporting material: data behind figures, machine-readable tables*

## 1. Introduction

The time domain radio sky hosts a diversity of astrophysical phenomena, including Galactic events such as novae, flare stars, X-ray binaries, cataclysmic variables, and magnetar flares; and extragalactic phenomena like supernovae, gamma-ray bursts, tidal disruption events, neutron star mergers, and active galactic nucleus (AGN) flares (Metzger et al. 2015; Fender et al. 2017). The most common phenomenon in time domain radio surveys is AGN variability (Thyagarajan et al. 2011; Mooley et al. 2016), which is hypothesized to be caused by shocks within a relativistic jet powered by an accreting supermassive black hole (Marscher & Gear 1985; Hughes et al. 1989), as well as due to extrinsic causes in our galaxy, such as scintillation (Walker 1998; Lovell et al. 2003; Koay et al. 2018).

Radio observations complement other wavelengths in the characterization of these events, for several reasons. Radio directly probes shocks and relativistic jets via synchrotron emission (Chevalier & Fransson 2017; Panessa et al. 2019), it is sensitive to dust-obscured phenomena missed in optical wave bands (Gal-Yam et al. 2006; Brunthaler et al. 2009), and it can provide fast, accurate localizations of events compared to wavelengths like X-ray and gamma rays. Monitoring the variable radio sky is therefore one of the main science drivers of future radio facilities, such as the Square Kilometer Array (SKA; Fender et al. 2015) and

the Next Generation Very Large Array (ngVLA; Murphy et al. 2018).

The last decade has seen a rise in the number of blind surveys attempting to systematically discover and characterize the rates of radio variability in the sky.<fnref n="8" /> These efforts were motivated by the emergence of powerful radio facilities with increasingly sensitive wide-field and wideband observing capabilities, such as the Karl G. Jansky Very Large Array (VLA; Perley et al. 2011), MeerKAT (Booth & Jonas 2012), and the Australian SKA Pathfinder Telescope (Murphy et al. 2013); and low-frequency facilities like the Low Frequency Array (LOFAR; van Haarlem et al. 2013) and Murchison Widefield Array (Bell et al. 2019). These studies have converged on a common result: the time-variable radio sky is relatively quiet. Only a handful of transients have been discovered in blind surveys (e.g., Gregory & Taylor 1986; Levinson et al. 2002; Bower et al. 2007; Bannister et al. 2011; Frail et al. 2012; Jaeger et al. 2012; Mooley et al. 2016; Radcliffe et al. 2019; Driessen et al. 2020). Most radio transient surveys have yielded nondetections, and found that only about 1% of the radio source population is variable (Carilli et al. 2003; Frail et al. 2003; de Vries et al. 2003).

<footnote n="1">¹Center for Data Intensive and Time Domain Astronomy, Department of Physics and Astronomy, Michigan State University, East Lansing, MI 48824, USA
sarbadhi@msu.edu</footnote>

<footnote n="2">²LESIA, Observatoire de Paris, CNRS, PSL, SU/UPD, Meudon, France</footnote>

<footnote n="3">³Sydney Institute for Astronomy, School of Physics, University of Sydney, Sydney, New South Wales 2006, Australia</footnote>

<footnote n="4">⁴Department of Astronomy, University of Wisconsin-Madison, Madison, WI 53706, USA</footnote>

<footnote n="5">⁵National Radio Astronomy Observatory, P.O. Box O, Socorro, NM 87801, USA</footnote>

<footnote n="6">⁶Astrophysics, Department of Physics, University of Oxford, Keble Road, Oxford, OX1 3RH, UK</footnote>

<footnote n="7">⁷Department of Astronomy, University of Cape Town, Private Bag X3, Rondenbosch, 7701, South Africa</footnote>

<footnote n="8">⁸For an exhaustive list, see http://www.tauceti.caltech.edu/kunal/radio-transient-surveys/.</footnote>
2004; Bower et al. 2010; Croft et al. 2010; Lazio et al. 2010; Bell et al. 2011; Bower et al. 2011; Croft et al. 2011; Ofek et al. 2011; Thyagarajan et al. 2011; Croft et al. 2013).

Time domain radio surveys face several challenges. First-generation radio surveys mainly probed radio sources brighter than 1 mJy, and often relied on archival all-sky data sets, such as the NRAO VLA Sky Survey (Condon et al. 1998) and the Faint Images of the Radio Sky at Twenty Centimeters Survey (Becker et al. 1995). The mismatched sensitivities, spatial resolutions, and absolute flux density calibrations encountered when comparing these surveys made it challenging to produce a complete sample of variables and transients. Imaging artifacts, as a result of improper data cleaning, calibration, and sidelobe contamination, were also identified as a potential source of false positive detections (Frail et al. 2012). More recent studies have utilized the newer generation of sensitive radio facilities, deep field extragalactic observations, and/or improvements in mosaicking and imaging techniques (Mooley et al. 2013; Hancock et al. 2016; Mooley et al. 2016; Bhandari et al. 2018; Radcliffe et al. 2019), but have still been forced to work with just a few epochs, as a compromise between epoch sensitivity, total time allocation, and survey area. The small number of epochs makes it difficult to characterize the temporal behavior of transients and variables occurring on a wide variety of timescales in the radio sky (Metzger et al. 2015).

Here, we fill in a portion of this under-explored parameter space with the CHILES Variable and Explosive Radio Dynamic Evolution Survey (CHILES VERDES), designed to look for variable phenomena within a single VLA pointing observed at L-band (1–2 GHz), using nearly 1000 hr of observations obtained between 2013 and 2019. The data were taken as part of the COSMOS H I Legacy Survey (CHILES), a survey of neutral hydrogen traced by the 21 cm emission line in the COSMOS field out to an unprecedented depth of $z \approx 0.5$ (Fernández et al. 2013; Dodson et al. 2016; Fernández et al. 2016; Hess et al. 2019; Luber et al. 2019; Blue Bird et al. 2020). The VLA Wideband Interferometric Digital Architecture correlator also provides simultaneous wideband continuum observations, which can be searched for synchrotron emission from energetic phenomena.

The CHILES VERDES survey is unique among GHz radio variable surveys, with its unprecedented combination of survey depth, cadence, and duration—sources are sampled with 172 epochs of observations every few days for 5.5 yr (except for the gaps between the B-configuration semesters), with each epoch reaching rms sensitivities of ~10 µJy beam⁻¹ per epoch. While the 22.5 radius field of view limits the recovered sample size, we are able to probe down to faint flux densities (<100 µJy) and sample the time domain characteristics of our sources on timescales spanning days to years. The survey also avoids the challenges associated with mosaics and multiple pointings (e.g., Mooley et al. 2013, 2016), and thus complements shallower, wide-field variability surveys. Finally, the COSMOS field is well supplemented with multiwavelength data, providing us with redshifts, multiwavelength luminosities, and detailed characterizations of the host galaxies of our sample (Scoville et al. 2007).

In this paper, we present a survey of variables carried out with the CHILES VERDES observations. A separate paper by A. J. Stewart et al. (2021, in preparation) will discuss the discovery and statistics of transient phenomena in CHILES VERDES. AGN variables are the dominant category of events discovered in time domain radio surveys (Thyagarajan et al. 2011), which vary on timescales similar to transients (Metzger

et al. 2015; Mooley et al. 2016), and can be misinterpreted as transients (e.g., Williams & Berger 2016). With more sensitive blind surveys, variables will likely dominate the background or confusion level against which transients are discovered (Carilli et al. 2003; Rowlinson et al. 2019). It is therefore important to understand the variability characteristics of the radio source population in blind surveys. Particularly important, and relatively less explored, is the variability in the submJy regime, where the source population is increasingly dominated by star-forming galaxies (Smolčić et al. 2017a) and the expected rate of all categories of transients is higher (Frail et al. 2012; Mooley et al. 2013; Metzger et al. 2015). The few surveys in this regime (e.g., Carilli et al. 2003; Radcliffe et al. 2019) show that the variable fraction in the submJy regime can be larger than that in the mJy regime. Extragalactic variability is also an important constraint on scintillation models and turbulence in the ionized interstellar medium (ISM), particularly off the Galactic plane, where pulsars are rare (Rickett 1990).

This paper is organized as follows. Section 2 describes the CHILES VERDES observations, data reduction, and imaging; Section 3 describes the search method and analysis of variables using an automated pipeline; and Section 4 describes the statistical properties, light curves, and structure functions (SFs) of the radio variables. Section 5 discusses the host galaxy properties, AGN properties, redshift distribution, and spectral indices of our sources from multiwavelength catalogs. In Section 6, we discuss the information learned about radio variability from our analyses, and we conclude in Section 7.

## 2. The Survey

### 2.1. The CHILES Survey

The COSMOS H I Legacy Extragalactic Survey, or CHILES, is an ambitious project using the VLA to obtain a deep field in the H I 21 cm spectral line (Fernández et al. 2013, 2016). The enormous correlator power of the upgraded (Jansky) VLA allows observers to obtain large bandwidths at high spectral resolutions (Perley et al. 2011), enabling a survey of neutral hydrogen from $z = 0.0–0.5$ in a single observation.

The CHILES pointing was chosen to be in the COSMOS field (Scoville et al. 2007), where it can enjoy significant multiwavelength supporting data (e.g., Davies et al. 2015; Civano et al. 2016; Andrews et al. 2017). It is centered at the J2000 position, R.A. = 10h01m24s.00, decl. = +02°21′00″.0. The single L-band pointing yields a field of view 0.5° in diameter (half-power beamwidth) at 1.5 GHz (we note that our search diameter for variables is slightly larger, see Section 4). The data were obtained in the VLA’s B-configuration, which yields a resolution of 4.3 at 1.5 GHz.

The CHILES data were obtained as both spectral line and continuum mode, the latter observed in full polarization mode and comprised of four spectral windows, each with 128 MHz bandwidth sampled by 64 channels. These spectral windows span the L-band and cover frequencies relatively free of radio frequency interference (RFI), centered at 1032, 1416, 1672, and 1800 MHz. The CHILES data were saved to disk with 8 s correlator integration time. The data were calibrated using 3C286 as a bandpass and absolute flux density calibrator, and J0943–0819 as a complex gain calibrator.
<table><caption>**Table 1**
Log of CHILES VERDES Observations</caption><thead><tr><th>Semester</th><th>Date Range</th><th>Number of Epochs</th><th>Total Time (hr)</th><th>Total Time on Source (hr)</th></tr></thead><tbody><tr><td>1</td><td>2013 Oct 25–2014 Jan 21</td><td>45</td><td>169.55</td><td>128.56</td></tr><tr><td>2</td><td>2015 Feb 25–2015 May 04</td><td>39</td><td>207.11</td><td>155.78</td></tr><tr><td>3</td><td>2016 May 16–2016 Sep 27</td><td>49</td><td>177.52</td><td>133.36</td></tr><tr><td>4</td><td>2017 Nov 11–2018 Jan 29</td><td>48</td><td>227.51</td><td>170.77</td></tr><tr><td>5</td><td>2019 Mar 01–2019 Apr 11</td><td>28</td><td>176.78</td><td>134.05</td></tr></tbody></table>

## 2.2. The CHILES VERDES Project

The CHILES VERDES project is the study of transient and variable sources within the CHILES deep field, using the continuum data set produced by CHILES. The data were observed over the time frames listed in Table 1. Each VLA B-configuration semester lasts approximately 4 months, and is followed by a 12-month break before the VLA reenters the B-configuration. During a given B-configuration semester, the CHILES observations were obtained, on average, every few days. The CHILES observations were obtained in scheduling blocks of durations ranging from 1–8 hr, with a median duration of 4.5 hr. We call each of these scheduling blocks an “epoch” throughout this paper.

The continuum data were reduced using a custom-developed scripted pipeline for the CHILES Continuum Polarization (or CHILES Con Pol) project, operated in CASA release version 4.7.2. Unlike CHILES VERDES, which uses total intensity (Stokes I) continuum maps, CHILES Con Pol is a partner project dedicated to making full polarization maps of the CHILES data in order to probe the cosmological evolution of polarized radio sources and intergalactic magnetic fields (Hales & Chiles Con Pol Collaboration 2014). The pipeline was run in a semiautomated manner by a single person (C. Hales), and each calibration step was manually inspected for quality assurance. High-quality flagging was performed using *pieflag* (Hales & Middelberg 2014), supplemented by CASA functionality provided by the custom tasks *antintflag3* (Hales 2016a), *plot3d* (Hales 2016b), and *interpgain* (Hales 2016c). An additional assessment of the calibration quality of each observation was carried out in Hales & Stephenson (2019) by applying the calibration solutions to 3C286, measuring the spectrum of fractional linear polarization, and comparing this with the known spectrum from the D-configuration observations by Perley & Butler (2013a). The pipeline was run twice for each observation, first to obtain calibrated data and a quick-look image within a few hours of coming off the telescope (for transient searches and potential optical follow-up), and the second at least two weeks after the observation, so as to incorporate GPS-derived ionospheric total electron content (TEC) data obtained from the International Global Navigation Satellite System Service, which has a latency period of up to two weeks. Our variable search described in Section 3 was carried out on the TEC-corrected total intensity data.

We created a time-averaged Stokes I image for each epoch using *tclean* in CASA (McMullin et al. 2007). We used multifrequency synthesis with *nterms = 2* and a reference frequency of 1.45 GHz. Each image is substantially larger than the primary beam half-power point, and calculated with w-projection and *wprojplanes = 128*, but is truncated at a radius corresponding to 20% of the sensitivity at the image’s

center. The images were created with a Briggs weighting of robust = 0.7, and the pixels are 1" across. Image deconvolution was carried out with the multi-scale multifrequency synthesis algorithm on scales of [0.5,15] pixels (Rau & Cornwell 2011), and a clean threshold was set at 6 × the image’s theoretical noise.

A total of 209 CHILES continuum images (totaling 960 hr of observations) were obtained between 2013 and 2019 from all the B-configuration scheduling blocks, of which we decided to use 172 for our CHILES VERDES search. The images have a typical spatial resolution of 4"5, with a range of 3"8–5"7 from image to image. The typical nonprimary beam-corrected image rms noise per epoch is about 10 µJy beam⁻¹, varying from image to image within the range 7–27 µJy beam⁻¹, depending on the on-source integration time per epoch. The rms noise was calculated inside a relatively uncrowded region of the CHILES field using the CASA task *imstat*. Of the rejected images, nine were obtained during the VLA B → A move time at the end of the 2016 semester, resulting in the systematic decrease in intensity of many CHILES VERDES sources, due to the changing uv coverage. The rest of the rejected images had various issues, from high rms noise to residual artifacts in the final images, that made them unsuitable for probing variability. The rejected images were affected either by corrupted data, high RFI levels due to daytime observations (particularly in 2016), or suboptimal calibration solutions obtained by the pipeline on short (~1 hr) epochs. A more detailed investigation will be carried out in the future in order to recover these epochs. For our purposes, the current set of 172 images has stable image quality, with the image rms being roughly a factor of 2 of the theoretical rms noise values in Figure 1. The theoretical rms noise was calculated using the VLA sensitivity equation,<fnref n="9" /> given the number of antennas used, the total time on source, and an assumed bandwidth of 0.5 GHz.

## 3. Finding Variables with the LOFAR Transients Pipeline

We identified potential sources and extracted their light curves from our radio images using the LOFAR Transients Pipeline (TraP; Swinbank et al. 2015). TraP can systematically extract radio sources across multiple images using a set of uniform and quantifiable selection criteria, and it has been used by a variety of contemporary radio variability surveys (e.g., Hobbs et al. 2016; Stewart et al. 2016; Rowlinson et al. 2019; Driessen et al. 2020).

The details of the TraP process are described in Swinbank et al. (2015), but for the convenience of the reader, we provide a brief description of the process. We ran TraP on our nonprimary beam-corrected images (excluding the rejected

<footnote n="9">⁹ https://science.nrao.edu/facilities/vla/docs/manuals/oss/performance/sensitivity</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** The ratio of the image rms noise ($\sigma_{rms,obs}$) to the theoretical rms noise ($\sigma_{rms,th}$) as a function of time for all 209 CHILES VERDES epochs, divided by the B-configuration semester. The 172 included epochs and the 38 excluded epochs are shown as black and red markers, respectively. The time-axis is the CHILES epoch in units of days, starting from our first observation, $t_0 = 2,456,588.5$ JD (2013 October 23).</figcaption></figure>

images described in Section 2.2), which have more uniform noise properties than the primary beam-corrected images and provide the opportunity to discover interesting variable sources farther away from the phase center. Once the final source catalog was obtained, primary beam correction was applied to the flux densities of the sources (Section 3.1). The first image in the time series searched by TraP was a deep reference image obtained by coadding the 2013 semester data. This was done so TraP could reliably identify sources with accurate positions, and to prevent spurious transient or variable detections from sources near the detection thresholds. TraP identified sources in each image by selecting pixels with flux densities above a $5\sigma$ threshold compared to the local rms, and then “growing” the source region to adjacent pixels above a $3\sigma$ threshold value. For each of these sources, TraP recorded the position, flux densities, and their respective uncertainties per image. TraP then calculated the flux densities of the sources by fitting an elliptical Gaussian at the recovered source location, keeping the width of the fitting beam fixed to the shape of the synthesized beam. It then matched each detected source per image with its counterpart in adjacent images, using the dimensionless de Ruiter radius ($r_{ij}$), which is the angular separation between source i and its counterpart j, normalized by their position uncertainties (de Ruiter et al. 1977). A source is said to have a counterpart in the subsequent image if their angular separation is less than the semimajor axis of the restoring beam, and if $r_{ij} \le r_s$ (where $r_s$ is a threshold set by the deruiter\_radius = 5.68, so that the probability of missing an associated source <10$^{-7}$). TraP is able to unambiguously match pointlike sources across images, but a small fraction of sources (mostly resolved) were identified as multiple entries in the final source catalog. These entries were automatically filtered out in the variable selection criteria further applied by us in Section 4.1.

### 3.1. Primary Beam Correction

After TraP has searched through all the epochs and recovered a source catalog, we corrected the integrated flux density of each source for the decreasing sensitivity of the primary beam with the increasing distance from the phase center. We followed the formalism of Perley (2016) by dividing the nonprimary beam-corrected flux densities by a factor $f(x)$, defined as

$$f(x) = 1 + a_1x + a_2x^2 + a_3x^3, \quad (1)$$

where $x = (R \nu_{obs})^2$, R is the distance of the source from the beam center in arcmin, and $\nu_{obs}$ is the observing frequency in GHz. In this paper, we use $a_1 = -1.343 \times 10^{-3}$, $a_2 = 6.579 \times 10^{-7}$, and $a_3 = -1.186 \times 10^{-10}$, which are the values for an average observing frequency of 1.465 GHz.<fnref n="10" /> All statistical calculations in the subsequent sections use the primary beam-corrected flux densities.

### 3.2. Identifying Variables

The sources recovered by TraP comprised a mixture of steady (nonvariable) and variable sources, and we differentiated them using statistics calculated from their radio light curves. A radio light curve is composed of a series of flux density measurements $F_i$, each with uncertainty $\sigma_{F,i}$; N is the number of measurements (=172 for most sources). We used these data to calculate the flux coefficient of variation (V) and the weighted reduced $\chi^2$ statistic ($\eta$), calculated as:

$$V = \frac{s}{\bar{F}} = \frac{1}{\bar{F}} \sqrt{\frac{N}{N-1} (\bar{F}^2 - \bar{F})^2} \quad (2)$$

and

$$\eta = \frac{1}{N-1} \sum_{i=1}^{N} \frac{(F_i - \xi_F)^2}{\sigma_{F,i}^2}. \quad (3)$$

Here, s refers to the standard deviation of the flux density measurements $F_i$, and $\bar{F}$ is the arithmetic mean of $F_i$. $\xi_F$ is the weighted mean of the flux density measurements, defined as:

$$\xi_F = \frac{\sum_{i=1}^{N} F_i / \sigma_{F,i}^2}{\sum_{i=1}^{N} 1 / \sigma_{F,i}^2}. \quad (4)$$

The coefficient V is equivalent to the fractional variability or modulation index parameter (m) used in previous transient surveys (e.g., Jenet & Gil 2003; Bell et al. 2014; Moooley et al. 2016). The $\eta$ parameter is similar to the reduced $\chi^2$ statistic, and quantifies the significance of the variability.

Together, V and $\eta$ have been shown to effectively separate the parameter space of variables from steady sources (Swinbank et al. 2015; Rowlinson et al. 2019). A common practice for identifying potential variable sources is to select those with V and $\eta$ above some threshold, defined using the

<footnote n="10">¹⁰ The primary beam correction formula and coefficients are adopted from http://www.aips.nrao.edu/cgi-bin/ZXHLP2.PL?PBCOR.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The threshold determination for $\eta$ as described in Section 3.2. The histogram shows the normalized log $\eta$-distribution of our 370 pointlike sources. The curve shows the best fit normal function. Sources with $\eta$ within the shaded 4$\sigma$ region are designated steady sources, while those with larger $\eta$ values are designated variable candidates.</figcaption></figure>

means and standard deviations of the $V$ and $\eta$ distributions (e.g., Rowlinson et al. 2016; Stewart et al. 2016). We adopt a similar strategy here, first determining a threshold on $\eta$ by fitting the log $\eta$ distribution with a normal function, with the mean ($\mu_{\log\eta}$) and standard deviation ($\sigma_{\log\eta}$) as free parameters. We follow the strategy of Rowlinson et al. (2019) by making a histogram of the log $\eta$ data, with the binning determined by using the Bayesian Blocks algorithm (Scargle et al. 2013) implemented in `astropy`.<fnref n="11" /> Bayesian Blocks is part of a family of algorithms that chooses the optimal binning to minimize the error of the histogram’s approximation of the data, with the added advantage of allowing variable bin widths. The histogram is then fitted with a normal function using a nonlinear least-squares method<fnref n="12" /> to obtain the best fit values for $\mu_{\log\eta}$ and $\sigma_{\log\eta}$ (Figure 2).

For this paper, we choose a $4\sigma$ threshold in log $\eta$. Sources with $\log\eta > (\mu_{\log\eta} + 4\sigma_{\log\eta})$ are considered variable candidates, while sources below this threshold are considered steady sources (Figure 2).

### 3.3. Structure Functions

The well-sampled light curves in CHILES VERDES, with their high signal-to-noise ratio (S/N), provide an excellent opportunity to examine the power spectra of radio variability in extragalactic sources, and to do this—possibly for the first time —in a detailed manner with a blind radio survey. In an attempt to understand these power spectra, we calculate SFs, which have commonly been used to quantify variability power spectra for blazars and quasars (e.g., Hughes et al. 1992), X-ray binaries (e.g., Plotkin et al. 2019), and radio scintillation (Simonetti et al. 1985). We apply the SF formalism of (Plotkin et al. 2019; hereafter P19) to our sources.

For a light curve with flux density measurements $F(t)$, the first-order SF is defined as

$$V(\tau) = \langle (F(t + \tau) - F(t))^2 \rangle, \qquad (5)$$

where $\tau$ is defined as the lagtime, or the time difference between two different epochs. For our discretely sampled light curves, we calculate $V(\tau)$ by taking each pair of flux density measurements in the light curve ($F_i, F_j$) at the times ($t_i, t_j$) and calculating $V_{ij}(\tau_{ij}) = (F_j - F_i)^2$, where $\tau_{ij} = t_j - t_i$. Similar to P19, we bin $V_{ij}$ into bins of $\tau_{ij}$, and calculate the mean $V_{ij}$ in that bin. The uncertainty in $V_{ij}$ is given by the standard error $=\sigma/\sqrt{N}$, where $\sigma$ is the standard deviation of the $V_{ij}$ values in the bin and $N$ is the number of points per bin of $\tau_{ij}$. The bin sizes were chosen to be approximately uniformly spaced in log lagtime, although note that because of the gaps between semesters, we do not have information on $V(\tau)$ for lagtimes between 123 and 363 days.

As explained in Hughes et al. (1992), the shape of the SF can reveal information about the underlying physical processes driving variability. At the shortest lagtimes, the SF plateaus to roughly twice the variance of the flux density uncertainties, and at the longest lagtimes to roughly twice the variance of the flux density measurements. The shape of the SF in between these extremes depends on the noise regime (e.g., white noise, flicker noise, shot noise; Hughes et al. 1992). The shortest lagtime at which the SF breaks from the white noise regime (where the SF flattens) corresponds to the shortest timescale of uncorrelated behavior (in the case of white noise) or the distribution of impulse response times (in the case of flicker noise). It is generally interpreted as the characteristic timescale of the underlying process driving variability.

We acknowledge, however, that the shapes of SFs are affected by the finite length of light curves, the large gaps in observations, and the random nature of variability, and therefore they may not accurately represent the true variability power spectrum of the source (Emmanoulopoulos et al. 2010). Furthermore, in many of the variables in Figures 8, 9, and 10, the SFs have complex shapes with no well-defined plateaus or breaks, and they cannot be easily fitted with simple functional forms to derive characteristic timescales. We therefore defer a more detailed analysis of the variability spectrum for future papers, and simply present the SFs alongside the light curves as a rough guide for the potential timescales present in the light curves.

## 4. CHILES VERDES Variable Candidates

### 4.1. Selecting Candidates

TraP searched a total of 172 images obtained between 2013 and 2019 that passed quality control, and yielded a catalog of 6821 sources with light curves. The distribution of these sources as a function of distance from the phase center is shown in Figure 3. The concentration of sources detected is highest near the phase center, and drops off rapidly as the sensitivity decreases with increasing distance from the phase center. For the subsequent analysis of the variable population in CHILES, we restrict ourselves to sources within a radius of 22′.5, which is $1.5 \times$ the primary beam half-width at half-maximum (HWHM) at 1.5 GHz ($HWHM_{PB} = 15'$, corresponding to about 17% primary beam sensitivity based on

<footnote n="11">¹¹ Implemented via `astropyvisualization` in Python 3+; see https://docs.astropy.org/en/stable/visualization/histogram.html.</footnote>

<footnote n="12">¹² Implemented via `scipyoptimize.curve_fit()`.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** The surface density of the sources found by TraP within the nonprimary beam-corrected images vs. their separation from the phase center. The vertical lines show the different distance cuts from the L-band primary beam phase center in units of the HWHM, $\theta_{PB} = 15'$, at 1.5 GHz. Since the surface density begins to flatten roughly at 1.5 HWHM<sub>PB</sub>, we include all sources within this radius for our variability analysis.</figcaption></figure>

Equation 1), beyond which the source density undergoes a noticeable flattening.

We find a total of 2713 sources inside the 1.5 HWHM<sub>PB</sub> region. Before investigating which of these sources are variable, we carry out two other selection cuts based on the positional uncertainties ($\Delta \theta_{pos}$) and average flux density uncertainties ($\langle\sigma_I\rangle$) returned by TraP, shown in Figure 4. The purpose of these cuts is to select for compact or pointlike sources, which are the only sources expected to vary on a ~6 yr timescale, and the sources which will have reliable photometry returned from TraP.

1. $\Delta\theta_{pos} \le 0''9$—We refer to the uncertainties of the R.A. and decl. of the sources returned by TraP as $\sigma_{R.A.}$ and $\sigma_{decl.}$, respectively. We define $\Delta\theta_{pos} = \sqrt{\sigma_{R.A.}\sigma_{decl.}}$, the geometric mean of the R.A. and decl. uncertainties. The average synthesized beam FWHM of our images is about 4''5, so for sources observed at S/N = 5, their positional accuracy should be ~4''5/5 ≈ 0''9. Nearly 85% of the sources in Figure 4 have $\Delta \theta_{pos} > 0''9$. Visual inspection revealed that the majority of these sources are faint (low S/N), while the rest are associated with bright extended sources with lobes, jets, and filaments, whose spatial scales would be too large to produce variability observable on our timescales. We therefore restrict ourselves to sources with $\Delta \theta_{pos} \le 0''9$.

2. $\langle\sigma_F\rangle \le 25~\mu\text{Jy}$—The average flux density measurement uncertainties of our sources are narrowly distributed around $\log\langle\sigma_F/\mu\text{Jy}\rangle = 1.2 \pm 0.04$ (or $\langle\sigma_F\rangle \approx 16~\mu\text{Jy}$), as seen in Figure 4, which is a typical rms noise in the individual epoch images. Visual inspection revealed that sources with high $\langle\sigma_F\rangle$ are mainly associated with bright extended features. We exclude these sources by setting an upper limit of $\langle\sigma_F\rangle < 25~\mu\text{Jy}$ (roughly five times the standard deviation of the log $\langle\sigma_F\rangle$ distribution) for our final catalog.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Mean flux density uncertainties ($\langle\sigma_F\rangle$) vs. position uncertainties ($\Delta\theta_{pos}$; as defined in Section 4.1) for all the CHILES sources identified by TraP inside the 1.5 × HWHM<sub>PB</sub> area. Sources with $\Delta\theta_{pos} < 0''9$ and $\langle\sigma_I\rangle < 0.025~\mu\text{Jy}$ (the gray region) are considered variable candidates for the remainder of the paper. The colorbar shows the fraction of flux density measurements >5× the local rms noise.</figcaption></figure>

We note some caveats here in our selection criteria. It is possible that some fraction of the faint sources that are excluded with $\Delta \theta_{pos} \le 0''9$ criteria may in fact exhibit variability that could be measured more accurately and reliably in deeper images obtained by coadding multiple epochs. Such deep images are being prepared as part of our follow-up transients paper, by A. J. Stewart et al. (2021, in preparation), in which we will comment on them further. Additionally, there may be some compact components of extended sources that could exhibit variability, but the surrounding diffuse emission leads them to having a lower S/N over the local background, meaning that they were therefore skipped by TraP's selection criteria. Making a catalog of such sources might be possible by reimaging the individual epochs, excluding the shorter baselines to filter out extended emission (e.g., from jets and filaments), and rerunning TraP on this image series. We defer additional analysis to a future paper. We note, however, that such extended sources form a small fraction of the source population in the CHILES field, and are unlikely to change the statistical results we infer for the variable population.

## 4.2. V–η Statistics

Our selection criteria in the previous section yielded a catalog of 370 sources inside the CHILES primary beam (Figure 5), in which we searched for variables. As described in Section 4.1, we selected for sources that appear compact in our images, so we will henceforth refer to this sample of 370 sources as “pointlike radio sources.” Based on V and η (Section 3.2), these sources are classified as either “moderate-variability” sources (showing significant, moderate-amplitude variability); “low-variability” sources (showing significant variability, but of relatively low amplitude); and nonvariable sources (showing no significant variability; in this paper, we will call them “steady” sources).

Figure 6 shows the V–η distribution of the 370 pointlike radio sources. The sample has means of $\bar{V} = 0.115 \pm 0.055$ and $\bar{\eta} = 0.43 \pm 0.13$. Most sources fall below the threshold
<figure><img src="image_6.png" /><figcaption>**Figure 5.** The CHILES field showing the 370 variable candidates in our paper (Section 4.1), from one of our images from 2017 November 16. The large dashed circle has a radius of 22.5, and is roughly 1.5 times the primary beam HWHM at 1.5 GHz. The red, blue, and yellow circles are pointlike radio sources identified by TraP that were designated moderate-variability, low-variability, and steady sources, respectively, as described in Section 4.2.</figcaption></figure>

$\eta = 1$ determined in Section 3.2, as expected for sources with no statistically significant variation (Rowlinson et al. 2016; Stewart et al. 2016; Rowlinson et al. 2019); these are called “steady” sources. We are left with 58 sources with $\eta > 1$, which are divided into “moderate-variability” and “low-variability” categories at the threshold of $V = 0.1$, corresponding to 10% variability in integrated flux density. We note that this threshold was chosen purely for discussion purposes, to separate sources that seemed to show more variability than others; the threshold is not a physical cutoff (in fact, as we show in Figure 10, some of the low-variability sources have light curves that are quite similar to those of the moderate-variability ones). Of the 58 sources, 18 are classified as moderate-variability, with $V \approx 0.1-0.3$, while 40 are classified as low-variability, with $V = 0.02-0.1$.

In Figure 7, we further assess the $V-\eta$ parameter space, by comparing the $V$ and $\eta$ values with the brightest per-epoch flux density measured in our CHILES VERDES light curves (called maximum flux density), and then comparing the ratio of this maximum flux density to the median flux density for that source, measured over its light curve. Not surprisingly (and as

also seen in Figure 6), the high $\eta$ sources are also the brightest sources; a source needs to be detected at high significance to make high-significance measurements of its variability. Variability ($V$) appears to be anticorrelated with flux density, with the low $V$ sources being some of the brightest objects. Again, this is not surprising—small-amplitude variations can be securely measured if a source is bright and detected at high significance. Faintier sources need to have higher $V$ values in order to reach a given $\eta$ value, compared to the brighter sources. High $V$, high $\eta$ sources would be easily detectable in our sample, but are not present.

The steady sources have flux densities $\lesssim 1.5$ mJy, and their $V$ values are inversely correlated with the maximum flux density, but positively correlated with the maximum-to-median flux density ratio. Both these trends can be explained if the epoch-by-epoch variability in steady sources is random noise on the order of the image rms. For example, consider a case where the deviation in flux density from the median, as well as the uncertainty in the flux density measurements, are both similar to the average rms of the images ($\sigma_{rms}$); then, $I_i - \xi_i \approx \sigma_{rms}$ and $\sigma_i \approx C\sigma_{rms}$, where $C$ is a constant ($\gtrsim 1$ to account for the fact
<figure><img src="image_7.png" /><figcaption>**Figure 6.** The V–η distribution of the 370 pointlike sources, from which we identified variables, color-coded by their median flux density. The sources at the left of the plot, with $η \lesssim 1$, do not have statistically significant variability and are designated steady. The higher $η$ sources are classified as moderate-variability and low-variability, differentiated at $V = 0.1$. The blue star shows the V–η position of the complex gain (phase) calibrator, J0943–0819.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** The values of V and η plotted against the maximum flux density attained during the observations and the ratio of maximum-to-median flux density. The points are color-coded so that steady sources are gold, low-variability sources are blue, and moderate-variability sources are magenta.</figcaption></figure>

that the local rms might be higher than the image rms). In that case, $η \lesssim 1$ from Equation (3), as observed for steady sources. Similarly for V as defined in Equation (2), if we assume $s \approx σ_{rms}$ and $I_{max} \approx (1-2)ξ_I$, then $\bar{V} \propto I_{max}^{-1}$. In other words, steady sources whose fluctuations are mainly driven by the rms noise of images will have lower V for brighter (and therefore higher S/N-detected) objects. The moderate- and low-variability sources, on the other hand, are as bright as ~10 mJy, implying that, contrary to the steady sources, both moderate-

and low-variability sources undergo flux density changes that are significantly higher than their measurement uncertainties (a fact also reflected in the SFs shown later).

The bottom panels of Figure 7 show that the sources span almost two orders of magnitude of maximum flux density (0.1–10 mJy), but that their maximum flux density is within a factor of 2 of their average brightness, consistent with previous deep-field variability surveys (e.g., Mooley et al. 2013; Bhandari et al. 2018; Radcliffe et al. 2019). Moderate-variability sources show a higher maximum-to-median flux density ratio than low-variability sources (they undergo up to a factor of 2 in brightening), explaining their higher V. Unlike the steady sources, both the moderate- and low-variability sources show V values that are uncorrelated with flux density, at least for maximum flux densities above 1 mJy, further suggesting that their variability is distinct from noise fluctuations in images. Faintier variable sources, however, show some correlation between V and maximum flux density, like the steady sources, and this is likely due to their lower S/N detections, which make them more vulnerable to image noise.

### 4.3. Individual Light Curves

In this section, we review the basic features of the individual light curves and SFs of our moderate-variability and low-variability sources, as shown in Figures 8, 9, and 10, and with measurements summarized in Tables 2 and 3. We only show 6 of the 40 low-variability sources, due to considerations of space (the full light curve tables will be available in electronic form), although the majority of them are similar to the first two sources in Figure 10. Hereafter, we give the prefix “CV” to our CHILES VERDES variables. The moderate-variability sources are numbered from 1–18 in order of increasing $η$, and the low-variability sources are numbered from 19–58, also in order of increasing $η$.

The variability patterns in the radio light curves are quite complex, with various episodes of brightening, dimming, and jittering on a variety of timescales, ranging from days to years. The majority (15 of 18) of the moderate-variability sources are fainter than 1 mJy in average brightness, with variability on the order of 10%–30% over the full course of the CHILES observations. Evidence of monthly variability is seen clearly in cases such as CV13, as well as CV10 and CV6. More jittery, short-timescale variability, however, is seen in cases like CV14. CV7 shows a rare example of a brightening event by more than a factor of 2 (in semester 2, between $t-t_0 = 520–540$ days) that is well-resolved by our observing cadence (i.e., longer than the timescale between subsequent epochs, but much shorter than a typical semester). Only three moderate-variability sources—CV16, CV17, and CV18—are brighter than 1 mJy, but they are among the most prominent moderate-variability sources in our sample. CV18, as we show later, is also the most luminous and distant radio variable source in our sample. CV18 shows variability on monthly timescales (e.g., the fluctuations in the first semester), but also a brightening and dimming on a longer, several-year timescale, as seen across the rest of the semesters. CV17 shows shorter-term variability at the level of ~5% on timescales of a few days, as compared to CV18, but also a systematic dimming by about 2 mJy (or 26%) over the full CHILES observation time. The SF does not plateau on the longest timescales, indicating that any characteristic timescale driving the long-term variability must be >6 yr. CV16 is the other >mJy source, and it mostly shows
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Light curves and SFs for sources designated "moderate-variability" (Section 4.2), in increasing order of $\eta$. Each light curve (left panel) shows the integrated flux density and their uncertainties over our five observing semesters, as measured by TraP, per CHILES epoch in units of days, starting from our first observation, $t_0 = 2,456,588.5$ JD (2013 October 23). The right panel shows the SF, calculated from the light curve, as described in Section 3.3. The upper and lower horizontal black dashed lines show the variance of the flux densities ($2\sigma_f^2$) and the variance of the flux density measurement uncertainties ($2\sigma_{err}^2$), respectively. The vertical orange dashed line shows the estimated scintillation timescale, and the green dashed line shows the minimum variability timescale needed to not exceed the inverse Compton limit on brightness temperatures (see Section 6.2).</figcaption></figure>

(The data used to create this figure are available.)
<figure><img src="image_10.png" /><figcaption>Figure 9. Light curves (continued), as in Figure 8.</figcaption></figure>

(The data used to create this figure are available.)

variability on timescales of about a week (11.5% in flux density), without any discernible long-term variability as observed in the other two cases.

The majority of the low-variability sources, true to their classification, have low-amplitude fluctuations at the few percent level, as evidenced by their light curves and SFs
<figure><img src="image_11.png" /><figcaption>**Figure 10.** The light curves of 6 of the 40 low-variability sources in our sample. The first two light curves are representative of the majority of the sample (most of the light curves look like these). The last four are sources that show interesting variability by eye, even though they do not pass the selection cut for being “variable.” The plots are as described as in Figure 8, and again $t_0 = 2,456,588.5$ JD.</figcaption></figure>

(e.g., CV22, Figure 10). Aside from their low-amplitude fluctuations, their variability characteristics are not all that dissimilar from those of the moderate-variability sources. For example, CV55 and CV57 are similar to CV18, one of the prominent variables discussed in the previous paragraph. CV55 shows monthly variability, seen prominently in semesters 1 and 4, but also a slower variability across the semesters, on a rough timescale of 5.4 yr, based on the SF. CV57 similarly shows variability on both short and long timescales. More short-term jittering variation in flux density is seen in CV50's light curve, with about 8% variation in flux density occurring on a timescale of about 10 days. CV46 is an example of a quiescent source that underwent a sudden short brightening. The source had a median flux density of 0.67 mJy and only 5% flux density variation until the middle of semester 3, then it brightened by about 44% to 1 mJy in about 2.5 weeks. It is unknown how far the source actually brightened, since we do not have reliable observations afterwards.

In summary, we find that the moderate-variability (and some of the low-variability) sources show a range of time series behaviors, with variations on both short and slow timescales,

and the occasional brightening during a quiescent phase. In contrast, the steady sources have light curves that do not show statistically significant variability over the time period in which we observed them, and an example of such a source is shown in Figure 12.

4.4. *Variability Due to Calibration*

One pertinent question is the amount of observed variability induced by the calibration process of the CHILES VERDES data. The flux densities in each CHILES epoch depend on the quality of the calibration solutions, as well as any intrinsic variability in the flux calibrator 3C286. We can get an idea of these effects from the changes in the flux density scale of the complex gain calibrator J0943-0819, whose flux density is measured in CHILES VERDES epochs, by comparison with 3C286, shown in Figure 11. J0943-0819 shows a low-level scatter in its flux density values of roughly 0.2%-0.3% per semester, except for the first semester, where it shows a slightly elevated level of variability (~0.6%). There is also a hint of a slow rise of 0.6% in the median flux density over 5.5 yr. While
<table><caption>**Table 2**
The Multiband and Multiwavelength Properties of Our 18 CHILES VERDES Moderate-variability Sources</caption><thead><tr><th>ID</th><th>R.A.</th><th>Decl.</th><th>⟨F⟩<br>(mJy)</th><th>V</th><th>η</th><th>Log L<sup>S10</sup><sub>1.4</sub><br>(W Hz<sup>-1</sup>)</th><th>Log L<sup>M16</sup><sub>X</sub><br>(erg s<sup>-1</sup>)</th><th>α</th><th>z</th><th>Classification</th></tr></thead><tbody><tr><td>CV1</td><td>10:00:06.6</td><td>+02:22:26.1</td><td>0.22</td><td>0.24</td><td>1.00</td><td>22.17</td><td>&lt;41.24</td><td>0.19 ± 0.47</td><td>0.222</td><td>...</td></tr><tr><td>CV2</td><td>10:02:15.8</td><td>+02:29:47.0</td><td>0.18</td><td>0.22</td><td>1.05</td><td>23.58</td><td>&lt;42.90</td><td>0.65 ± 0.43</td><td>1.146</td><td>MLAGN</td></tr><tr><td>CV3</td><td>10:00:21.8</td><td>+02:12:20.2</td><td>0.41</td><td>0.11</td><td>1.12</td><td>23.37</td><td>&lt;41.89</td><td>-0.49 ± 0.22</td><td>0.426</td><td>MLAGN</td></tr><tr><td>CV4</td><td>10:00:37.6</td><td>+02:29:49.0</td><td>0.14</td><td>0.27</td><td>1.22</td><td>23.33</td><td>42.3</td><td>0.06 ± 0.5</td><td>0.671</td><td>HLAGN</td></tr><tr><td>CV5</td><td>10:01:33.6</td><td>+02:27:49.5</td><td>0.09</td><td>0.30</td><td>1.23</td><td>22.22</td><td>&lt;41.21</td><td>-0.83 ± 0.58</td><td>0.215</td><td>MLAGN</td></tr><tr><td>CV6</td><td>10:01:10.8</td><td>+02:02:04.1</td><td>0.48</td><td>0.14</td><td>1.37</td><td>23.84</td><td>42.87</td><td>0.99 ± 0.27</td><td>0.972</td><td>HLAGN</td></tr><tr><td>CV7</td><td>10:00:34.0</td><td>+02:26:45.7</td><td>0.14</td><td>0.26</td><td>1.47</td><td>23.49</td><td>&lt;42.65</td><td>-0.02 ± 0.51</td><td>0.9</td><td>MLAGN</td></tr><tr><td>CV8</td><td>10:00:13.9</td><td>+02:22:49.4</td><td>0.22</td><td>0.24</td><td>1.56</td><td>22.63</td><td>42.25</td><td>0.73 ± 0.46</td><td>0.347</td><td>HLAGN</td></tr><tr><td>CV9</td><td>10:00:47.6</td><td>+02:09:58.5</td><td>0.42</td><td>0.11</td><td>1.80</td><td>24.31</td><td>&lt;42.34</td><td>-1.19 ± 0.21</td><td>0.669</td><td>MLAGN</td></tr><tr><td>CV10</td><td>10:01:04.4</td><td>+02:04:37.1</td><td>0.32</td><td>0.18</td><td>2.04</td><td>23.51</td><td>&lt;42.34</td><td>-0.14 ± 0.34</td><td>0.668</td><td>MLAGN</td></tr><tr><td>CV11</td><td>10:01:32.4</td><td>+02:16:37.6</td><td>0.19</td><td>0.15</td><td>2.10</td><td>24.08</td><td>&lt;43.08</td><td>-0.34 ± 0.29</td><td>1.364</td><td>MLAGN</td></tr><tr><td>CV12</td><td>10:01:49.8</td><td>+02:28:32.0</td><td>0.17</td><td>0.21</td><td>2.61</td><td>23.76</td><td>&lt;42.99</td><td>0.24 ± 0.4</td><td>1.248</td><td>MLAGN</td></tr><tr><td>CV13</td><td>10:01:42.6</td><td>+02:07:52.9</td><td>0.71</td><td>0.11</td><td>6.75</td><td>24.39</td><td>&lt;42.64</td><td>-0.51 ± 0.21</td><td>0.89</td><td>MLAGN</td></tr><tr><td>CV14</td><td>10:01:39.8</td><td>+02:25:48.6</td><td>0.5</td><td>0.12</td><td>8.11</td><td>22.34</td><td>41.79</td><td>-0.398 ± 0.244</td><td>0.124</td><td>MLAGN</td></tr><tr><td>CV15</td><td>10:01:24.0</td><td>+02:20:04.7</td><td>0.56</td><td>0.13</td><td>14.97</td><td>23.92</td><td>&lt;42.34</td><td>-0.406 ± 0.248</td><td>0.666</td><td>MLAGN</td></tr><tr><td>CV16</td><td>10:01:53.5</td><td>+02:11:52.4</td><td>2.87</td><td>0.11</td><td>144.57</td><td>23.07</td><td>43.01</td><td>-0.28 ± 0.21</td><td>0.405</td><td>HLAGN</td></tr><tr><td>CV17</td><td>10:02:24.1</td><td>+02:16:21.3</td><td>6.18</td><td>0.11</td><td>343.04</td><td>23.55</td><td>41.13</td><td>-0.67 ± 0.21</td><td>0.122</td><td>...</td></tr><tr><td>CV18</td><td>10:01:20.1</td><td>+02:34:43.6</td><td>8.43</td><td>0.11</td><td>803.20</td><td>25.77</td><td>43.77</td><td>0.097 ± 0.218</td><td>1.555</td><td>HLAGN</td></tr></tbody></table>

Note. The column headers are as follows: “ID”: the CHILES VERDES ID assigned to the object; “R.A.,decl.”: the International Celestial Reference System R.A. and decl. found by TraP; “⟨F⟩”: the median 1.4 GHz flux density measured by CHILES VERDES for all 172 epochs; “V, η”: the coefficient of variability and the reduced χ² statistic defined in Section 4.2; “Log L<sup>S10</sup><sub>1.4</sub>”: the VLA-COSMOS 1.4 GHz luminosity from Smolčić et al. (2017a); “Log L<sup>M16</sup><sub>X</sub>”: the 0.5–8 keV X-ray luminosity from Marchesi et al. (2016); “α”: the 1.4–3 GHz spectral index; “z”: redshift; and “Classification”: the classification of the objects as HLAGNs, MLAGNs, or star-forming galaxies, as discussed in Section 5.

(This table is available in machine-readable form.)

we cannot rule out the variability of J0943–0819 itself, we note that this calibrator’s flux destiny scale is set by 3C286. Therefore, any very low-level variability exhibited by 3C286 could have been transferred to J0943–0819. A similar scatter is also observed in the polarization time series of J0943–0819, which was also attributed to variability in 3C286 (Hales & Stephenson 2019). Such low-level (<1%) variability in 3C286 is consistent with the archival measurements in Perley & Butler (2013a, 2013b, 2017).

The variability in the calibrators will likely make a relatively small contribution to the variability observed in our targets. This is particularly evident from the V-η values of J0943–0819 (the blue star in Figure 6), compared to the CHILES VERDES sources. J0943–0819 has a high η because it is bright (2.7 Jy at 1.45 GHz), observed at a high S/N, but it has V ≈ 0.007 calculated over the full CHILES VERDES observations—more than an order of magnitude lower than the V of the CHILES targets. We also compare the fractional deviation in flux density from the median flux density per epoch (Δ F/F) for J0943–0819, and for the summed flux density of all sources brighter than 1 mJy. We define the quantity as:

$$ \Delta F/F = \frac{F_i - \langle F \rangle}{\langle F \rangle}, \qquad (6) $$

where Fᵢ is the flux density in an epoch i (of the calibrator or the mJy-bright CHILES sources) and ⟨F⟩ is the median flux density during a particular CHILES semester. Figure 13 shows that the scatter in Δ F/F for the CHILES VERDES sources is larger than that for the gain calibrator (consistent with Figure 6), and mostly uncorrelated with the gain calibrator variability. An exception is the first semester, where some anomalously high flux densities of the calibrator (also seen in Figure 11) are correlated with the variability of the CHILES

sources. However, even if we exclude the 2013 semester, the impact on the V-η distribution in Figure 6 is minimal, with an average difference of 6% in the V-η measurements with and without the 2013 semester. We therefore conclude that, statistically speaking, issues with calibration are unlikely to be producing the observed variability in CHILES VERDES.

# 5. Multiwavelength Characteristics

## 5.1. AGN Classification

We can determine the physical characteristics of our radio variable sources using the plethora of multiwavelength data available in the COSMOS field.<fnref n="13" /> A large number of surveys have targeted the COSMOS field for characterizing galaxies, but in this paper we will closely follow the work of (Smolčić et al. 2017a, henceforth S17), who provide the largest catalog of z < 6 radio-selected galaxies classified using observations spanning radio to X-ray wavelengths. S17 uses the VLA-COSMOS 3 GHz survey, which detected 10,830 sources above 5σ with 384 hr of VLA data in the full 2.2 sq deg of the COSMOS field (Smolčić et al. 2017b). S17 obtained 8035 multiwavelength counterparts for the 3 GHz sources in the unmasked areas of the COSMOS field (the masked areas mainly demarcate regions with bright foreground stars, gaps between chips, etc. in data at nonradio wavelengths).

The radio sources in extragalactic surveys are mainly star-forming galaxies and AGNs (Condon 1992). The radio emissions from star-forming galaxies mainly trace the integrated cosmic ray population produced by supernovae and supernova remnants in massive star-forming regions (Murphy 2009; Delhaize et al. 2017). For AGNs, the radio

<footnote n="13">¹³ https://cosmos.astro.caltech.edu/page/datasets</footnote>
<table><caption>**Table 3**
The Multiband and Multiwavelength Properties of Our 40 CHILES VERDES Low-variability Sources</caption><thead><tr><th>ID</th><th>R.A.</th><th>Decl.</th><th>⟨F⟩<br>(mJy)</th><th>V</th><th>η</th><th>Log L<sup>S10</sup><sub>I,4</sub><br>(W Hz<sup>-1</sup>)</th><th>Log L<sup>M16</sup><sub>X<sup>-1</sup></sub><br>(erg s<sup>-1</sup>)</th><th>α</th><th>z</th><th>Classification</th></tr></thead><tbody><tr><td>CV19</td><td>10:02:27.1</td><td>+02:21:19.2</td><td>0.95</td><td>0.039</td><td>1.01</td><td>24.89</td><td>&lt;42.95</td><td>-0.83 ± 0.11</td><td>1.202</td><td>MLAGN</td></tr><tr><td>CV20</td><td>10:01:13.6</td><td>+02:06:53.6</td><td>0.41</td><td>0.074</td><td>1.03</td><td>23.2</td><td>&lt;41.61</td><td>-0.55 ± 0.15</td><td>0.323</td><td>MLAGN</td></tr><tr><td>CV21</td><td>10:01:01.3</td><td>+02:01:18.0</td><td>1.31</td><td>0.052</td><td>1.05</td><td>25.72</td><td>&lt;43.41</td><td>-1.28 ± 0.12</td><td>1.89</td><td>MLAGN</td></tr><tr><td>CV22</td><td>10:01:33.0</td><td>+02:21:09.8</td><td>1.05</td><td>0.019</td><td>1.05</td><td>24.36</td><td>&lt;42.39</td><td>-0.96 ± 0.08</td><td>0.7</td><td>MLAGN</td></tr><tr><td>CV23</td><td>10:01:28.0</td><td>+02:40:29.3</td><td>1.25</td><td>0.045</td><td>1.12</td><td>24.92</td><td>42.79</td><td>-0.75 ± 0.11</td><td>1.108</td><td>HLAGN</td></tr><tr><td>CV24</td><td>10:01:52.5</td><td>+02:19:54.2</td><td>0.30</td><td>0.081</td><td>1.36</td><td>22.87</td><td>&lt;41.34</td><td>-0.71 ± 0.17</td><td>0.246</td><td>HLAGN</td></tr><tr><td>CV25</td><td>10:01:04.5</td><td>+02:02:03.5</td><td>1.51</td><td>0.044</td><td>1.48</td><td>...</td><td>...</td><td>-1.09 ± 0.11</td><td>...</td><td>...</td></tr><tr><td>CV26</td><td>10:02:10.1</td><td>+02:16:38.0</td><td>0.60</td><td>0.061</td><td>1.48</td><td>24.28</td><td>&lt;42.63</td><td>-0.53 ± 0.13</td><td>0.879</td><td>MLAGN</td></tr><tr><td>CV27</td><td>10:02:32.7</td><td>+02:20:00.0</td><td>0.51</td><td>0.095</td><td>1.58</td><td>23.86</td><td>&lt;42.25</td><td>-0.34 ± 0.20</td><td>0.607</td><td>MLAGN</td></tr><tr><td>CV28</td><td>10:00:14.2</td><td>+02:13:12.1</td><td>1.34</td><td>0.051</td><td>1.66</td><td>25.35</td><td>43.63</td><td>-1.80 ± 0.12</td><td>1.14</td><td>HLAGN</td></tr><tr><td>CV29</td><td>10:00:46.9</td><td>+02:07:26.5</td><td>1.85</td><td>0.027</td><td>1.72</td><td>25.31</td><td>&lt;42.99</td><td>-1.02 ± 0.08</td><td>1.252</td><td>HLAGN</td></tr><tr><td>CV30</td><td>10:01:10.6</td><td>+02:24:58.4</td><td>0.63</td><td>0.045</td><td>1.83</td><td>25.71</td><td>&lt;43.30</td><td>-1.83 ± 0.11</td><td>1.698</td><td>MLAGN</td></tr><tr><td>CV31</td><td>10:01:49.6</td><td>+02:33:34.8</td><td>1.89</td><td>0.025</td><td>1.84</td><td>24.88</td><td>&lt;42.52</td><td>-1.29 ± 0.08</td><td>0.792</td><td>MLAGN</td></tr><tr><td>CV32</td><td>10:02:28.4</td><td>+02:32:30.2</td><td>1.31</td><td>0.063</td><td>1.86</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CV33</td><td>10:01:02.4</td><td>+02:05:27.9</td><td>1.58</td><td>0.036</td><td>2.23</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CV34</td><td>10:01:34.2</td><td>+02:09:17.5</td><td>0.71</td><td>0.052</td><td>2.46</td><td>25.32</td><td>&lt;43.57</td><td>-0.76 ± 0.12</td><td>2.222</td><td>MLAGN</td></tr><tr><td>CV35</td><td>10:02:30.7</td><td>+02:09:21.6</td><td>1.32</td><td>0.077</td><td>2.54</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CV36</td><td>10:01:12.1</td><td>+02:41:06.6</td><td>2.21</td><td>0.048</td><td>3.01</td><td>25.46</td><td>43.12</td><td>-1.30 ± 0.11</td><td>1.258</td><td>HLAGN</td></tr><tr><td>CV37</td><td>10:01:09.0</td><td>+02:28:15.8</td><td>0.61</td><td>0.064</td><td>3.2</td><td>24.32</td><td>&lt;42.66</td><td>-0.57 ± 0.15</td><td>0.911</td><td>MLAGN</td></tr><tr><td>CV38</td><td>10:01:44.8</td><td>+02:04:09.1</td><td>0.89</td><td>0.087</td><td>3.28</td><td>24.5</td><td>&lt;42.72</td><td>-0.43 ± 0.18</td><td>0.965</td><td>MLAGN</td></tr><tr><td>CV39</td><td>10:01:36.5</td><td>+02:26:41.6</td><td>0.69</td><td>0.083</td><td>3.52</td><td>22.97</td><td>41.06</td><td>-1.25 ± 0.15</td><td>0.123</td><td>MLAGN</td></tr><tr><td>CV40</td><td>10:01:59.8</td><td>+02:39:04.8</td><td>3.24</td><td>0.037</td><td>3.73</td><td>25.04</td><td>&lt;42.53</td><td>-0.79 ± 0.10</td><td>0.8</td><td>MLAGN</td></tr><tr><td>CV41</td><td>10:02:29.9</td><td>+02:32:25.1</td><td>2.47</td><td>0.047</td><td>3.93</td><td>24.55</td><td>43.37</td><td>-0.76 ± 0.09</td><td>0.432</td><td>HLAGN</td></tr><tr><td>CV42</td><td>10:02:29.8</td><td>+02:09:10.3</td><td>2.21</td><td>0.061</td><td>4.03</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CV43</td><td>10:02:14.5</td><td>+02:35:10.0</td><td>4.28</td><td>0.03</td><td>5.24</td><td>26.07</td><td>&lt;43.51</td><td>-0.67 ± 0.09</td><td>2.099</td><td>MLAGN</td></tr><tr><td>CV44</td><td>10:01:24.1</td><td>+02:17:06.3</td><td>1.67</td><td>0.027</td><td>5.69</td><td>25.4</td><td>&lt;43.29</td><td>-0.70 ± 0.09</td><td>1.677</td><td>MLAGN</td></tr><tr><td>CV45</td><td>10:00:05.4</td><td>+02:30:29.0</td><td>3.36</td><td>0.052</td><td>6.13</td><td>24.96</td><td>&lt;42.46</td><td>-0.77 ± 0.12</td><td>0.746</td><td>MLAGN</td></tr><tr><td>CV46</td><td>10:01:51.5</td><td>+02:25:32.2</td><td>0.70</td><td>0.078</td><td>6.20</td><td>24.22</td><td>42.55</td><td>-0.65 ± 0.16</td><td>0.827</td><td>HLAGN</td></tr><tr><td>CV47</td><td>10:00:00.6</td><td>+02:15:31.0</td><td>3.91</td><td>0.041</td><td>6.7</td><td>25.68</td><td>44.47</td><td>0.33 ± 0.10</td><td>2.45</td><td>HLAGN</td></tr><tr><td>CV48</td><td>10:01:09.3</td><td>+02:17:21.6</td><td>3.64</td><td>0.019</td><td>10.32</td><td>26.35</td><td>&lt;43.72</td><td>-0.98 ± 0.07</td><td>2.582</td><td>MLAGN</td></tr><tr><td>CV49</td><td>10:02:28.8</td><td>+02:17:21.8</td><td>3.41</td><td>0.041</td><td>14.03</td><td>26.07</td><td>43.68</td><td>-0.57 ± 0.10</td><td>2.625</td><td>HLAGN</td></tr><tr><td>CV50</td><td>10:00:49.8</td><td>+02:16:54.8</td><td>0.95</td><td>0.084</td><td>14.27</td><td>24.53</td><td>&lt;42.67</td><td>-0.53 ± 0.17</td><td>0.917</td><td>MLAGN</td></tr><tr><td>CV51</td><td>09:59:58.0</td><td>+02:18:09.3</td><td>4.20</td><td>0.064</td><td>14.43</td><td>26.18</td><td>&lt;43.30</td><td>-1.02 ± 0.14</td><td>1.698</td><td>MLAGN</td></tr><tr><td>CV52</td><td>10:01:47.3</td><td>+02:03:14.1</td><td>6.01</td><td>0.032</td><td>14.70</td><td>24.36</td><td>42.16</td><td>-0.74 ± 0.09</td><td>0.323</td><td>HLAGN</td></tr><tr><td>CV53</td><td>10:02:01.2</td><td>+02:13:27.0</td><td>3.35</td><td>0.043</td><td>21.40</td><td>...</td><td>...</td><td>-0.99 ± 0.10</td><td>...</td><td>...</td></tr><tr><td>CV54</td><td>10:01:22.4</td><td>+02:01:11.9</td><td>9.04</td><td>0.04</td><td>31.90</td><td>25</td><td>&lt;41.88</td><td>...</td><td>0.425</td><td>HLAGN</td></tr><tr><td>CV55</td><td>10:00:16.6</td><td>+02:26:38.3</td><td>4.15</td><td>0.065</td><td>44.22</td><td>26.3</td><td>&lt;43.66</td><td>-0.77 ± 0.14</td><td>2.436</td><td>HLAGN</td></tr><tr><td>CV56</td><td>10:01:31.1</td><td>+02:29:24.7</td><td>3.44</td><td>0.07</td><td>105.14</td><td>24.31</td><td>42.81</td><td>-0.78 ± 0.15</td><td>0.349</td><td>HLAGN</td></tr><tr><td>CV57</td><td>10:02:09.1</td><td>+02:16:02.4</td><td>5.58</td><td>0.062</td><td>162.70</td><td>25.34</td><td>&lt;42.68</td><td>-0.71 ± 0.13</td><td>0.928</td><td>...</td></tr><tr><td>CV58</td><td>10:01:31.4</td><td>+02:26:39.5</td><td>9.39</td><td>0.062</td><td>460.92</td><td>24.79</td><td>&lt;41.68</td><td>-0.76 ± 0.13</td><td>0.348</td><td>MLAGN</td></tr></tbody></table>

Note. The column headers have the same meanings as in Table 2.

(This table is available in machine-readable form.)

emissions mainly trace the synchrotron emissions produced by the collimated relativistic jets from the central supermassive black hole (SMBH), although outflows and accretion disk coronae can also contribute to the radio emissions in radio-quiet systems (Panessa et al. 2019). Such radio-selected AGNs have been differentiated into two main categories, based on the appearance of excitation lines in optical spectra: high-excitation and low-excitation (Smolčić 2009; Heckman & Best 2014). The two categories are believed to trace different phases of AGN–galaxy coevolution, with high-excitation sources emitting across a wide range of the electromagnetic spectrum, consistent with radiatively efficient SMBH accretion at high bolometric luminosity; these are often associated with galaxies in the “green valley” of the color–magnitude diagram. Meanwhile, low-excitation sources exhibit lower bolometric luminosities, as expected for radiatively inefficient SMBH accretion,

and tend to be associated with the red passive sequence of galaxies (Smolčić 2009; Heckman & Best 2014).

S17 further expands on this classification, using the abundant multiwavelength information available in the COSMOS field, including the latest (z<sup>++</sup>YJHK<sub>s</sub>) photometry from COSMOS2015 (Laigle et al. 2016); the i-band catalog from Capak et al. (2007); the 3.6 µm Spitzer(IRAC)-COSMOS catalog from Sanders et al. (2007); the FIR Herschel data (Griffin et al. 2010; Poglitsch et al. 2010; Oliver et al. 2012), the submillimeter photometry from various observational campaigns (Bertoldi et al. 2007; Scott et al. 2008; Aretxaga et al. 2011; Smolčić et al. 2012; Casey et al. 2013; Miettinen et al. 2015); and the X-ray catalog of point sources from the Chandra-COSMOS Legacy Survey (Marchesi et al. 2016). These data sets were used to classify the 3 GHz radio sources into three broad categories: high-to-moderate luminosity AGNs (HLAGNs), moderate-to-low luminosity AGNs (MLAGNs),
<figure><img src="image_12.png" /><figcaption>**Figure 11.** The flux density of the gain calibrator J0943–0819, measured using the CASA (version 4.7.2) task *fluxscale* at the end of the calibrations.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** An example steady or nonvariable source in our CHILES VERDES data set with $\eta < 1$.</figcaption></figure>

and star-forming galaxies. HLAGNs and MLAGNs are expected to be high-redshift counterparts to the high- and low-excitation galaxies discussed earlier, and the classification scheme is based on bolometric luminosity, which also serves as a proxy for the SMBH accretion rate (Figure 9 in S17). We refer the reader to S17 and Delvecchio et al. (2017) for details of the selection criteria and individual categories, and provide only a brief description of each category below.

VLA-COSMOS 3 GHz sources were classified as HLAGNs if they satisfied at least one of three criteria—(i) 0.5–8 keV X-ray luminosity >10<sup>42</sup> erg s<sup>-1</sup> (Szokoly et al. 2004); (ii) mid-IR colors consistent with the AGN selection criteria of Donley et al. (2012); or (iii) a significant AGN contribution to the optical-to-millimeter spectral energy distribution (SED; Delvecchio et al. 2017). The individual criteria are limited in their completeness, but together they are believed to effectively select AGNs with high radiative luminosities (Figure 9 of S17).

The remaining sources not classified as HLAGNs were grouped into MLAGNs and star-forming galaxies based on their restframe colors and the presence of “radio excess,” i.e., 1.4 GHz luminosity exceeding (by more than 3σ) the redshift-dependent L<sub>1.4</sub>-star formation rate relation (SFR), where the SFR was determined from the integrated IR (1–1000 µm) luminosity from the SED fitting (Delhaize et al. 2017; Delvecchio et al. 2017). A non-HLAGN object was classified as a star-forming galaxy depending on its dust-corrected restframe colors (i.e., either $M_{NUV} - M_{r^+} < 3.5$ mag or $M_{NUV} - M_{r^+} > 3.5$ mag but with detection in Herschel bands) and if it exhibited no radio excess. The remaining radio sources were classified as MLAGNs, and consist of two subcategories: (i) galaxies deemed to be red/quiescent (based on their dust-

corrected restframe colors being $M_{NUV} - M_{r^+} > 3.5$ mag and with nondetection in Herschel bands); and (ii) star-forming galaxies (also based on restframe colors, but with a >3σ radio excess).

In addition to the AGN classifications, we also discuss the redshifts and radio spectral indices from the VLA-COSMOS 1.4 and 3 GHz surveys, which are described in the subsequent sections.

## 5.2. Cross-matching with S17 Using TOPCAT

We cross-matched our sources with the S17 catalog using TOPCAT, a GUI analysis package for working with large tabular data sets (Taylor et al. 2005). The CHILES VERDES sources were matched to their nearest S17 sources within a radius of 1", which is roughly equal (but rounded-off) to the TraP-based positional errors discussed in Section 4.1. We note that S17 took into account possible astrometric uncertainties and false match probabilities (e.g., due to chance coincidences) when constructing the catalog of multiwavelength counterparts for the VLA-COSMOS radio sources. Since the CHILES VERDES and S17 catalogs are both based on VLA observations of a common patch of sky, we did not employ any additional cross-matching verification beyond visual checks and simply identifying the nearest neighbors for our sources.

## 5.3. Results

### 5.3.1. Host Galaxy Properties

Of the 370 pointlike CHILES VERDES sources, TOPCAT found cross-matches for 316 sources in the S17 catalog; their demographics are shown in Figure 14. The remaining 54
<figure><img src="image_14.png" /><figcaption>Figure 13. Comparison between the fractional deviation from median flux density (Δ F/F) for the gain calibrator J0943–0819 and the sum of all sources brighter than 1 mJy.</figcaption></figure>

sources did not cross-match for several reasons, such as the higher image resolution of the VLA-COSMOS 3 GHz survey (i.e., some unresolved CHILES VERDES sources were resolved into two or more point sources in 3 GHz that were spatially separated on scales larger than our search radius of 1"), the faintness of the 3 GHz counterpart, proximity to a bright, saturated object (thus appearing in a masked region), or a genuine lack of optical/near-IR counterparts at the source location. However, a majority of the CHILES VERDES pointlike sources have counterparts, and form a large enough sample to be sufficiently representative of the different source populations that we discuss below.

Of the variable sources, cross-matches were obtained for 16 of the 18 moderate-variability sources, and 33 of the 40 low-variability sources. The CHILES VERDES variables fall into both HLAGN and MLAGN categories, as seen in Figure 14, with MLAGN featuring more variables in both the moderate- and low-variability categories as compared to HLAGN. We note that two of our brightest variables in the sample—CV16 and CV18—are HLAGNs. This may suggest that radio variability is correlated with the SMBH accretion rate, although our relatively small sample size makes this a low-significance result.

A notable, if not surprising, result is that all of the sources that are “clean” (without radio excess) star-forming galaxies are also “steady,” i.e., not showing any radio variability. This result is consistent with the expectation that star formation occurs over extended temporal and spatial scales in galaxies, and therefore the integrated radio emission should not be variable on our timescales. This result is also confirmation that our thresholds on V and η (described in Section 3.2) are effective at separating variable sources from otherwise steady objects.

### 5.3.2. Spectral Indices

To calculate spectral indices, we collected flux density information from the VLA-COSMOS 1.4 GHz (Schinnerer et al. 2010) and VLA-COSMOS 3 GHz (Smolčić et al. 2017b) surveys. Although these data sets have different image resolutions, the majority of our sources are pointlike (as enforced in Section 4). The 1.4 GHz data set consists of sources detected above 5σ (rms of 12 µJy beam⁻¹) in the combined observations from the VLA-COSMOS Large and Deep Projects, while the 3 GHz data set consists of sources detected above 5σ (rms of 2.3 µJy). Both surveys provide the integrated flux densities (thus allowing for extended sources) and their uncertainties.

While we could have acquired spectral indices from the luminosities available in the S17 catalog, it did not list the flux

density uncertainties for the 1.4 and 3 GHz measurements. We therefore cross-matched our 370 CHILES VERDES sources with the 1.4 and 3 GHz data sets, using the same 1" radius circle criteria as described in Section 5.2. We only selected targets that cross-matched to within 1" for both the 1.4 and 3 GHz surveys, bringing us to a total of 339 sources. For a spectral index of the form $F_ν ∝ ν^α$, we calculate the uncertainty on α as:

$$ \sigma_{\alpha} = \frac{[(\sigma_{F_1}/F_1)^2 + (\sigma_{F_2}/F_2)^2]^{1/2}}{\ln(\nu_1/\nu_2)}, \qquad (7) $$

where $ν_1, ν_2$ are the observed frequencies (1.4 GHz and 3 GHz, respectively); $F_1, F_2$ are the observed flux densities at these frequencies; and $σ_{F_1}, σ_{F_2}$ are the uncertainties on $F_1, F_2$. We assume that the uncertainties in the flux density measurements are due to statistical uncertainties ($σ_{F,st}$), as well as the intrinsic variability of the source ($σ_{F,var}$), added in quadrature

$$ \sigma_F = \sqrt{\sigma_{F,st}^2 + \sigma_{F}^2}. \qquad (8) $$

The statistical uncertainties are listed in the catalogs. Since V is the fractional rms variation in the median flux density of our sources, we assume $σ_{F,var} = VF$. We also assume, for the sake of simplicity, that the variability V at 3 GHz is the same as that at 1.4 GHz for any source.

The relation between α and variability V of the 339 CHILES VERDES sources is shown in Figure 15. We designate as “flat spectrum” any source with $α > -0.5$, consistent with the traditional definition (Condon 1984). The steady sources show a range of spectral indices, with a mean and standard deviation of $α = -0.88 ± 0.45$, and with nearly 97% of the sources being consistent with or below $α = -0.5$ within their 1σ uncertainties. This is consistent with the majority of steady sources being optically thin synchrotron emitters.

The variable sources, however, show a correlation between spectral index and variability, with the moderate-variability sources showing higher α on average as compared to the low-variability sources (Figure 15(a)). No such correlation is apparent between α and η, however (Figure 15(b)). We note that the uncertainties in α are larger for sources with higher variability, because they are dominated by V. About 9 of the 18 moderate-variability sources have a flat spectrum at the 1σ level (i.e., the median α is at least 1σ above -0.5), and 4 have flat spectra at the 2σ level. In comparison, only 1 of the 34 low-variability sources (for which we obtained spectral indices) has a flat spectral index (although it is flat at the 2σ level). Such flat-spectrum sources could be core-dominated systems;
<figure><img src="image_15.png" /><figcaption>**Figure 14.** The distribution of CHILES VERDES sources in the V–η parameter space, separated by the multiwavelength classification of the host galaxies. Each panel has the same “steady–moderate-variability–low-variability” parameter spaces as defined in Figure 6, with the colors indicating the log 1.4 GHz median flux density. For details of the host galaxy classification, refer to Section 5.1.</figcaption></figure>

systems with aligned jets, such as blazars; or young jets, such as those associated with gigahertz-peaked spectrum sources (O'Dea 1998; Hovatta et al. 2014).

### 5.3.3. Redshifts

One of the unique advantages of deep-drilling surveys is the opportunity to probe the cosmological evolution of sources. This is nicely demonstrated in the S17 catalog, which provides radio-selected sources with reliable redshifts out to $z \approx 6$, sampling most of the cosmic history of star formation and black hole growth, and resolving the peak epochs at $z \approx 1.5–3$ (Shankar et al. 2009; Heckman & Best 2014; Madau & Dickinson 2014; Smolčić et al. 2017c).

Here we investigate the 1.4 GHz luminosity versus redshift distribution of the CHILES VERDES sources. Redshifts for all sources were compiled by S17 from an exhaustive list of photometric and spectroscopic observations in the COSMOS field, obtained from both the VLA-COSMOS 3 GHz team and archival studies, and we refer the reader to S17 and Delvecchio et al. (2017) for the references and details of the measurements. The 1.4 GHz luminosities were measured by S17 using the restframe 1.4 GHz flux densities and spectral indices discussed in Section 5.3.2.

We show the 1.4 GHz luminosity–redshift distribution of our CHILES VERDES sources in Figure 16. Our sources span nearly 6 orders of magnitude in luminosity (from $10^{21}–10^{27}$ W Hz$^{-1}$). The steady sources, being more numerous, span most of the observed range of redshifts from $z = 0.07$ to $z = 3.86$. The moderate-variability sources are observed between $z = 0.12–1.55$, while the low-variability sources extend to higher redshifts, with the maximum $z = 2.63$. The luminosity evolution shows the characteristic shape of flux-limited surveys, whereby less luminous sources are discovered at lower redshifts. The distribution of radio luminosities at a given redshift range is similar to that of the steady sources. However, the moderate-variability sources are generally less luminous as compared to the low-variability sources. The low-variability sources are among the most luminous sources in the sample across most of the observed redshift range. We note that, compared to the full S17 sample, the CHILES VERDES

sources are generally more luminous, mainly because of the shallower flux density limit of the individual CHILES VERDES epochs ($3\sigma$ limits in the range of the 21–81 µJy beam$^{-1}$), compared to 36 µJy for the VLA-COSMOS 1.4 GHz survey, and 6.9 µJy for the VLA-COSMOS 3 GHz survey. Similarly, while the full S17 sample extends to $z \approx 6$, our catalog of steady sources is limited to $z < 4$, since there are very few sources (only 95 of the 8035 S17 sources with redshifts) beyond $z = 4$.

## 6. Discussion

### 6.1. Variability Statistics

With 370 pointlike sources within our CHILES VERDES search region, a tally of 18 moderate-variability sources suggests a variability fraction of $18/370 \approx 4.9\%$ (we excluded low-variability sources to be consistent with previous radio variability surveys that generally include sources with variability amplitude above 10%). We estimate a Poisson statistical uncertainty $=\sqrt{18}/370 \approx 1.1\%$. Our selection criteria may also have added uncertainty to the variable fraction, but we expect this to be of the order of the Poisson uncertainty (e.g., if we relax our positional uncertainty criteria for pointlike sources in Section 4.1 to $\Delta\theta_{pos} \le 1.5''$, we get 616 sources and 21 moderate-variability sources, making the variability fraction $\approx 3.4\%$). For the mJy sources, we have 3 moderate-variability objects out of 32 pointlike sources, giving a variability fraction for sources brighter than 1 mJy as $(9.3 \pm 5.4)\%$. In the submJy regime, we have 15 moderate-variability objects out of 338 pointlike sources, giving a variability fraction for submJy sources of $(4.4 \pm 1.1)\%$.

These fractions are larger than those in most other 1.4 GHz surveys in the literature, where the typical variability fraction is constrained to less than 1% (e.g., Croft et al. 2010; Thyagarajan et al. 2011; Mooley et al. 2013; Hancock et al. 2016; Bhandari et al. 2018; similar statistics are echoed in other frequency bands as well).<fnref n="14" /> One caveat applying to comparisons with

<footnote n="14">¹⁴ We restricted our comparison to 1.4 GHz surveys for consistency, since variability can depend on the frequency band observed. The low variability fraction, however, has been reported in other frequencies as well, and a full list of studies is available here: http://www.tauceti.caltech.edu/kunal/radio-transient-surveys/.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 15.** (a) 1.4–3 GHz spectral index $\alpha$ vs. $V$ measured for our CHILES VERDES sources, color-coded by their variability categories from Figure 6. (b) The same as in panel (a), except with $\alpha$ plotted against $\eta$ for our CHILES VERDES sources. We note that error bars for the steady sources are not shown for clarity, and error bars for the variable sources include the effects of variability.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 16.** 1.4 GHz luminosity vs. redshift for the CHILES VERDES sources, obtained from cross-matching with the S17 catalog. The source categories are labeled as usual.</figcaption></figure>

previous variability surveys involves the different sensitivities, timescales, and variability thresholds used by the studies. Additionally, because of our single pointing and the effect of cosmic variance, the source counts may differ by about 10%–35% in another similar patch of the sky (Heywood et al. 2013).

Regardless, our result agrees with those of previous findings, in that only a small fraction of the extragalactic radio population exhibits any statistically significant variability. Our survey is likely picking up a variability fraction >1% because of the combination of depth (which helps to probe the rich extragalactic population between 0.1–1 mJy at >5σ), cadence (allowing us to resolve weekly–monthly timescale variability), and survey length (making our survey also

sensitive to the slow year-long variability that appears to be common in AGNs). Previous studies have hinted at the possible dependence of the variability fraction on factors like sensitivity (Carilli et al. 2003; Thyagarajan et al. 2011; Hancock et al. 2016; Radcliffe et al. 2019), cadence, and the timescales probed (Hodge et al. 2013; Radcliffe et al. 2019).

We also found that no variable source brightened by more than a factor of 2 during the course of our observations, which is consistent with the previously mentioned deep-field studies (e.g., Mooley et al. 2013; Hancock et al. 2016; Bhandari et al. 2018; Radcliffe et al. 2019), implying that dramatic brightening events, such as those expected from flares and renewed jet activity in AGNs discovered in wider-area surveys (Mooley et al. 2016), are relatively rare. Given the 370 variability candidates recovered by our pipeline, and an observing baseline of ~5.5 yr, we estimate that an object brightening by more than a factor of 2 happens at a rate of less than twice per square degree per 1000 yr.

## 6.2. Sources of Variability Observed in CHILES VERDES

From Figures 8, 9, and 10, it is clear to see that the radio variables in CHILES VERDES vary on a range of timescales, from a few days/week to years, and that individual sources themselves may be associated with variability on multiple timescales (e.g., CV17).

Variability in radio sources has been attributed to both extrinsic effects (such as scintillation by the Galactic ISM) and intrinsic effects (i.e., arising from shocks in the black hole jet). We can place a rough estimate on the contribution from scintillation using the models of Walker (1998) and the distribution of ISM electron density from Cordes & Lazio (2002). Walker (1998) parameterizes the amount of scattering induced in a radio wave front by the turbulent ISM (approximated as a thin phase-scattering screen) in terms of the scattering strength $\zeta = (\nu_0/\nu)^{17/10}$, where $\nu$ is the observing frequency and $\nu_0$ is the line of sight-dependent transition frequency between weak ($\zeta \ll 1$) and strong ($\zeta \gg 1$)
scattering regimes. Physically, broad band scintillations are greatest when $\nu = \nu_0$, and therefore $\zeta = 1$, implying a modulation index (corresponding to $V$) = 1. In the direction of the CHILES field, $\nu_0 = 9.4$ GHz,<fnref n="15" /> which means $\zeta \approx 24 \gg 1$ for $\nu = 1.45$ GHz (our observing frequency), putting us in the strong scattering regime. In this regime, scintillation can be refractive or diffractive, but only the former will be captured in broad band observations. For refractive scintillation, the modulation index $V = (\nu/\nu_0)^{17/30} \approx 0.35$ and the timescale of variation $t_r \sim (2 \text{ hr})(\nu_0/\nu)^{11/5} \approx 5$ days for $\nu = 1.45$ GHz (marked as an orange dashed line in the SFs plotted in Figures 8, 9, and 10). Hints of a slope break on timescales $\sim t_r$ are seen in the SFs of some of the variables, and are particularly prominent in the bright sources CV16 and CV17 (Figure 9), the fainter source CV4 (Figure 8), as well as the low-variability sources CV38 and CV50 (Figure 10). Several sources show variability on timescales (based on plateaus in their SFs) almost an order of magnitude above the scintillation timescale, such as CV3, CV5, CV6, CV10, CV13, CV17, and CV18, along with some of the low-variability sources. Such variability on longer timescales is therefore not likely entirely driven by scintillation.

Another constraint on the variability timescale comes from the fact that the brightness temperature for incoherent synchrotron sources is limited to $<3 \times 10^{11}$ K, due to inverse Compton losses (Kellermann & Pauliny-Toth 1969; Readhead 1994). The brightness temperature depends on the source angular extent ($\theta$) as $T_B \propto \theta^{-2}$. By a light travel-time argument, $\theta \le c\tau$, where $\tau$ is the variability timescale of the source, putting a lower limit on $T_B$ (Pietka et al. 2015). However, this lower limit needs to be below the inverse Compton limit, which puts a physical lower limit on $\tau$. Following Bhandari et al. (2018) Equation 4, this limit is defined as:

$$ \tau_{\min} \ge \sqrt{\frac{\Delta FD^2}{2k_B \nu^2 T_{B,max}}}, \qquad (9) $$

where $k_B$ is the Boltzmann constant, $\Delta F$ is the amplitude of flux density variation, $D$ is the luminosity distance to the source, $\nu$ is the observing frequency, and $T_{B,max} = 3 \times 10^{11}$ K is the maximum allowed brightness temperature. Note that this does not include any correction for Doppler-boosting effects. Assuming $\Delta F = V\bar{T}$, we show $\tau_{\min}$ in the SFs of our variables (the green dashed lines on the SFs in Figures 8, 9, and 10). We find $\tau_{\min}$ larger than the characteristic timescales in many of the sources, like CV6, CV10, CV11, CV13, and CV18, and in all the low-variability sources in Figure 10. If a source is varying on these timescales, it would exceed the brightness temperature limit.

This implies that the aforementioned monthly-yearly variability is either driven by long-timescale scintillation, or by Doppler-boosted intrinsic variability. The scintillation timescale derived earlier assumes that the extent of radio emission from sources ($\theta_s$) is smaller than the first Fresnel zone $\theta_r$ (the characteristic scale over which 1 rad of phase-delay is induced in the wave front). If components larger than the central engine (e.g., star-forming regions, lobes) also make a significant contribution to the radio emission, then $\theta_s$ could be larger than $\theta_r$, and the scintillation timescale becomes larger by a factor ($\theta_s/\theta_r$). However, the modulation index on these

<figure><img src="image_18.png" /><figcaption>Figure 17. The effect of sparse sampling demonstrated by the $V-\eta$ distribution for $10^4$ randomly drawn flux densities from the light curve of CV18 in sets of 3, 5, and 15 epochs (shown in colors). The top and side histograms show the distributions of $V$ and $\eta$, respectively, for the randomly drawn points (see Section 6.3 for details). The black circle with crosshairs shows the measured $V$ and $\eta$ for CV18 from all 172 epochs.</figcaption></figure>

timescales will also be decreased by a factor $(\theta_r/\theta_s)^{7/6}$, whereas from the SFs (e.g., CV18 or CV13) we see that longer timescales have more power than shorter ones. Scintillation on such long timescales at the latitude of the CHILES pointing (42°4) at 1.4 GHz is ruled out by more recent scintillation timescale measurements (see Hancock et al. 2019). The remaining possibility is that the variability is intrinsic, but affected by Doppler boosting, which increases $T_B$ for a given source extent by a multiplicative factor $\delta^3$ that depends on the jet kinematics (Cohen et al. 2003; Kellermann et al. 2007; Hovatta et al. 2009), and therefore would decrease the variability timescale according to Equation 9. Doppler boosting has been observed in sources for even modest inclinations (Zensus 1997). Further comparison with the very long baseline interferometry measurements of the CHILES field (e.g., Herrera Ruiz et al. 2017) could help to measure the Doppler factors and deduce the presence of boosted variability in our objects.

## 6.3. Reflections on Future Surveys

Our work underscores the importance of future deep-field surveys for studying transients and variables. While wide-field surveys discover more transient events, the ability of deep-field surveys to sample radio variability with more epochs over the same time allocation is crucial for characterizing the temporal behavior of the radio sky. This is particularly important for characterizing extragalactic variables that show flux density variations on a range of timescales; sampling a few epochs can lead to less precise determinations of their variability. As an example, we randomly sample $10^4$ sets of $N=3, 5, \text{ and } 15$ epochs from the light curve of CV18 (mimicking the epoch intervals of previous studies), and calculate $V$ and $\eta$ for each set (shown in Figure 17). In all cases, the mean $V$ and $\eta$ are similar to the measured values

<footnote n="15">¹⁵ Based on https://www.nrl.navy.mil/rsd/RORF/ne2001/.</footnote>
from 172 epochs, but the percentage uncertainty (the standard deviation divided by the mean) in the measured V and η is 51% and 93%, respectively, for 3 epochs, and drops to 18% and 33%, respectively, for 15 epochs. Statistically, therefore, there is a higher chance of over/underestimating the significance and variability of a source with only a few epochs. Although the number of epochs per source depends on the survey area and total time allocation, it would be advisable to have ≥10 epochs to accurately characterize the parameter space of variables, and, by extension, steady sources as well as bright transients.

## 7. Conclusion

We have conducted a unique deep-field radio variability survey in a 0.44 deg² area of the data-rich COSMOS field using 1–2 GHz VLA continuum observations. Our sources were observed from 2013 October to 2019 April, and sampled with 172 epochs, with rms noise in the range of 7–27 µJy and synthesized beam resolution of ~4.′′3. On average, the cadence was <2 days (except for the year-long gaps between each B-configuration semester), and the observing time totaled almost 960 hr of on-source time.

A total of 370 pointlike sources were detected using TraP, of which 18 were designated “moderate-variability” sources (>10% variability), 40 were designated “low-variability” sources (<10%, but still significant, variability), and the rest were designated “steady” sources (no significant variability). The classification of these sources was carried out using the coefficient of variation V (equivalent to the rms fractional flux density variation, or the modulation index) and η (the reduced χ² statistic). The fraction of >80 µJy radio sources that are variable at the level of 10% or more is 4.9%. We characterized our sources using their radio light curves, SFs, spectral indices obtained by matching the VLA-COSMOS 1.4 GHz and 3 GHz catalogs, and multiwavelength characteristics obtained by cross-matching with the S17 catalog. This catalog provided redshifts, intrinsic luminosities, and host galaxy classifications for VLA-COSMOS 3 GHz sources by cross-matching with archival photometric, spectroscopic, and multiwavelength catalogs. Radio sources were classified as HLAGNs, MLAGNs, and star-forming galaxies, based on the thresholds in X-ray luminosity, mid-IR colors, optical-to-millimeter SEDs, and radio excess over the expectations from IR-derived star formation rates.

The 18 moderate-variability sources show complex time series behavior, including slow month–year long variability (e.g., CV17, CV18) and faster day–week timescale variability (e.g., CV16). The flux density variations are observed at the level of ≥10% for the brightest (mJy) sources and up to ~30% in the submJy regime. The SFs of these sources show the expected plateaus at short timescales, due to measurement errors, and at long timescales in the white noise/flicker noise regime of variability (roughly equal to the variance of flux densities). The intermediate timescales show a variety of shapes and breaks, possibly indicating the superposition of processes occurring at multiple characteristic timescales. We note that such detailed SFs were possible in a blind survey, due to the large number of epochs sampling the light curves.

The moderate-variability sources generally have a maximum-to-median flux density ratio between 1.3–2, indicating that brighter flaring events in AGNs are quite rare. The host

galaxies for all the moderate-variability sources have AGN signatures, the majority of which are also associated with radio excess, which further corroborates the presence of a relativistic jet producing radio emission. The host galaxies are a mixture of HLAGNs and MLAGNs, spanning a redshift range of z = 0.22–1.56. About half of the moderate-variability sources are consistent with having flat spectra in the 1–4 GHz range ($F_ν ∝ ν^{\alpha}$, with $\alpha > -0.5$), and are candidates for core-dominated systems, blazars, and/or young jets associated with gigahertz-peaked spectra.

Most of the 40 low-variability sources show flux density variation at levels of 2%–8%. Physically there are similarities between the host galaxy properties of low- and moderate-variability sources, in that they are associated with HLAGNs and MLAGNs in similar proportions. However, low-variability sources are generally brighter and more luminous at radio wavelengths, extend to farther distances ($z \approx 2.8$), have smaller maximum-to-median flux density ratios (≤1.5), and have generally steeper 1–4 GHz spectra. A few interesting objects (e.g., CV55, CV50, and CV57) show variability patterns similar to our objects classified as moderate-variability sources, indicating that V–η thresholds can miss interesting objects, and it is important to manually inspect the light curves of sources near such thresholds.

Steady sources form the bulk of our time series sample. The largest group comprises star-forming galaxies (45%), and the rest of the sample contains roughly equal proportions of HLAGNs and MLAGNs. This indicates that not all AGNs are associated with radio variability (at least with the cadence, depths, and duration explored in this paper).

For the 58 CHILES VERDES sources that demonstrate significant variability, their light curves can be explained as a combination of scintillation at short timescales and intrinsic Doppler-boosted variability on longer timescales. Our observing frequency, bandwidth, and line of sight off the Galactic plane implies that any scintillation would be refractive, with variability at the level of <35% over a characteristic timescale ~5 days, based on the Walker (1998) model. Several variables show breaks in their SFs at these timescales, indicating the role of scintillation. However, breaks at longer timescales of months and years are also observed, indicating intrinsic AGN variability. However, these observed timescales are lower than the minimum timescale needed to remain below the inverse Compton brightness temperature limit of $3 × 10^{11}$ K. This implies that scintillation still dominates at longer timescales, in excess of model predictions, or that the variability is intrinsic but Doppler-boosted.

We are grateful to Jacqueline van Gorkom for her leadership of the CHILES survey, Ian Heywood for discussions about source counts, Vernesa Smolčić and Ivan Delvecchio for guidance with the VLA-COSMOS 3 GHz catalog, Dillon Dong and Kristina Nyland for discussions about radio variability, and Kirill Sokolovsky for discussions on Doppler boosting. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This paper includes observations obtained at the Southern Astrophysical Research (SOAR) telescope, which is a joint project of the Ministério da Ciência, Tecnologia, Inovações e Comunicações do Brasil (MCTIC/LNA), the US National Science Foundation’s (NSF’s) NOIRLab (NOIRLab), the University of North Carolina at
Chapel Hill (UNC), and Michigan State University (MSU). This research has made use of the NASA/IPAC Infrared Science Archive, which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology. This research has made use of NASA's Astrophysics Data System Bibliographic Services.

S.S., E.T., and L.C. are grateful for support from NSF grants AST-1412549 and AST-1907790. J.S. and S.S. acknowledge support from the Packard Foundation. C.P. and E.W. are grateful for support from NSF grant AST-1413099.

Facilities: VLA, SOAR, HST, IRSA.

Software: numpy (van der Walt et al. 2011), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), astropy (Astropy Collaboration et al. 2013, 2018), CASA (McMullin et al. 2007).

## Appendix

### Redshift Determination for CV17

CV17 is one of the brightest variable objects, but its redshift was not listed in the COSMOS catalogs. This is likely due to its location in a masked area of optical data, because of a bright foreground source (Figure 18). We obtained an optical spectrum of CV17 with the Goodman Spectrograph (Clemens et al. 2004) on the SOAR telescope on 2019 November 4. Two 20-minute exposures were obtained, using the 4001 mm⁻¹ grating and covering a wavelength range of ~3850–7850 Å at a resolution of about 5.9 Å. The spectra were optimally extracted in the usual way and combined, with a first-order flux calibration applied. The resulting spectrum is shown in Figure 18. It shows the characteristic absorption lines of an intermediate-to-old stellar population, as well as narrow emission lines. By cross-correlating the emission lines with an emission-line galaxy template from XCSAO,<fnref n="16" /> we find z = 0.122.

<figure><img src="image_19.png" /><figcaption>Figure 18. w? > SOAR spectrum of CV17 showing the location of redshifted Hα + [N II]. The inset shows an HST/ACS image of CV17 and the bright foreground object.</figcaption></figure>

## ORCID iDs

Sumit K. Sarbadhicary https://orcid.org/0000-0002-4781-7291

Evangelia Tremou https://orcid.org/0000-0002-4039-6703

Adam J. Stewart https://orcid.org/0000-0001-8026-5903

Laura Chomiuk https://orcid.org/0000-0002-8400-3705

Charee Peters https://orcid.org/0000-0001-5826-6803

Chris Hales https://orcid.org/0000-0002-3733-2565

Jay Strader https://orcid.org/0000-0002-1468-9668

Emmanuel Momjian https://orcid.org/0000-0003-3168-5922

## References

Andrews, S. K., Driver, S. P., Davies, L. J. M., et al. 2017, MNRAS, 464, 1569

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

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Bannister, K. W., Murphy, T., Gaensler, B. M., Hunstead, R. W., & Chatterjee, S. 2011, MNRAS, 412, 634

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

Bell, M. E., Fender, R. P., Swinbank, J., et al. 2011, MNRAS, 415, 2

Bell, M. E., Murphy, T., Hancock, P. J., et al. 2019, MNRAS, 482, 2484

Bell, M. E., Murphy, T., Kaplan, D. L., et al. 2014, MNRAS, 438, 352

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

Bhandari, S., Bannister, K. W., Murphy, T., et al. 2018, MNRAS, 478, 1784

Blue Bird, J., Davis, J., Luber, N., et al. 2020, MNRAS, 492, 153

Booth, R. S., & Jonas, J. L. 2012, AfrSk, 16, 101

Bower, G. C., Saul, D., Bloom, J. S., et al. 2007, ApJ, 666, 346

Bower, G. C., Whysong, D., Blair, S., et al. 2011, ApJ, 739, 76

Bower, R. G., Vernon, I., Goldstein, M., et al. 2010, MNRAS, 407, 2017

Brunthaler, A., Menten, K. M., Reid, M. J., et al. 2009, A&A, 499, L17

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

Carilli, C. L., Ivison, R. J., & Frail, D. A. 2003, ApJ, 590, 192

Casey, C. M., Chen, C.-C., Cowie, L. L., et al. 2013, MNRAS, 436, 1919

Chevalier, R. A., & Fransson, C. 2017, in Handbook of Supernovae, ed. A. W. Alsabti & P. Murdin (Cham: Springer), 875

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

Clemens, J. C., Crain, J. A., & Anderson, R. 2004, Proc. SPIE, 5492, 331

Cohen, M. H., Russo, M. A., Homan, D. C., et al. 2003, in ASP Conf. Ser. 300, Variability and Velocity of Superluminal Sources, ed. J. A. Zensus, M. H. Cohen, & E. Ros (San Francisco, CA: ASP), 177

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

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

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

Cordes, J. M., & Lazio, T. J. W. 2002, arXiv:astro-ph/0207156

Croft, S., Bower, G. C., Ackermann, R., et al. 2010, ApJ, 719, 45

Croft, S., Bower, G. C., Keating, G., et al. 2011, ApJ, 731, 34

Croft, S., Bower, G. C., & Whysong, D. 2013, ApJ, 762, 93

Davies, L. J. M., Driver, S. P., Robotham, A. S. G., et al. 2015, MNRAS, 447, 1014

de Ruiter, H. R., Willis, A. G., & Arp, H. C. 1977, A&AS, 28, 211

de Vries, W. H., Becker, R. H., White, R. L., & Helfand, D. J. 2004, AJ, 127, 2565

Delhaize, J., Smolčić, V., Delvecchio, I., et al. 2017, A&A, 602, A4

Delvecchio, I., Smolčić, V., Zamorani, G., et al. 2017, A&A, 602, A3

Dodson, R., Vinsen, K., Wu, C., et al. 2016, A&C, 14, 8

Donley, J. L., Koekemoer, A. M., Brusa, M., et al. 2012, ApJ, 748, 142

Driessen, L. N., McDonald, I., Buckley, D. A. H., et al. 2020, MNRAS, 491, 560

Emmanoulopoulos, D., McHardy, I. M., & Uttley, P. 2010, MNRAS, 404, 931

Fender, R., Stewart, A., Macquart, J. P., et al. 2015, in Proc. of Advancing Astrophysics with the Square Kilometre Array (AASKA14) (Giardini Naxos, Italy), 51

Fender, R., Woudt, P. A., Armstrong, R., et al. 2017, in Proc. of MeerKAT Science: On the Pathway to the SKA (Stellenbosch, South Africa), 13

Fernández, X., Gim, H. B., van Gorkom, J. H., et al. 2016, ApJL, 824, L1

Fernández, X., van Gorkom, J. H., Hess, K. M., et al. 2013, ApJL, 770, L29

Frail, D. A., Kulkarni, S. R., Berger, E., & Wieringa, M. H. 2003, AJ, 125, 2299

Frail, D. A., Kulkarni, S. R., Ofek, E. O., Bower, G. C., & Nakar, E. 2012, ApJ, 747, 70

<footnote n="16">¹⁶ http://tdc-www.harvard.edu/iraf/rvsao/Templates</footnote>
Gal-Yam, A., Ofek, E. O., Poznanski, D., et al. 2006, ApJ, 639, 331

Gregory, P. C., & Taylor, A. R. 1986, AJ, 92, 371

Griffin, M. J., Abergel, A., Abreu, A., et al. 2010, A&A, 518, L3

Hales, C. A. & (Chiles Con Pol Collaboration) 2014, Exascale Radio Astronomy, 2, 50301

Hales, C. A. 2016a, antintflag v1.0, Zenodo, doi:10.5281/zenodo.163546

Hales, C. A. 2016b, interpgain v1.1, Zenodo, doi:10.5281/zenodo.163002

Hales, C. A. 2016c, plot3d v1.5, Zenodo doi:10.5281/zenodo.162986

Hales, C. A., & Middelberg, E. 2014, pieflag: CASA task to efficiently flag bad data. Astrophysics Source Code Library, ascl:1408.014

Hales, C. A., & Stephenson, N. 2019, EVLA Memo 208: JVLA calibration stability at L-band over 5.5 years, https://library.nrao.edu/public/memos/evla/EVLAM_208.pdf

Hancock, P. J., Charlton, E. G., Macquart, J.-P., & Hurley-Walker, N. 2019, arXiv:1907.08395

Hancock, P. J., Drury, J. A., Bell, M. E., Murphy, T., & Gaensler, B. M. 2016, MNRAS, 461, 3314

Heckman, T. M., & Best, P. N. 2014, ARA&A, 52, 589

Herrera Ruiz, N., Middelberg, E., Deller, A., et al. 2017, A&A, 607, A132

Hess, K. M., Luber, N. M., Fernández, X., et al. 2019, MNRAS, 484, 2234

Heywood, I., Jarvis, M. J., & Condon, J. J. 2013, MNRAS, 432, 2625

Hobbs, G., Heywood, I., Bell, M. E., et al. 2016, MNRAS, 456, 3948

Hodge, J. A., Becker, R. H., White, R. L., & Richards, G. T. 2013, ApJ, 769, 125

Hovatta, T., Aller, M. F., Aller, H. D., et al. 2014, AJ, 147, 143

Hovatta, T., Valtaoja, E., Tornikoski, M., & Lähteenmäki, A. 2009, A&A, 494, 527

Hughes, P. A., Aller, H. D., & Aller, M. F. 1989, ApJ, 341, 68

Hughes, P. A., Aller, H. D., & Aller, M. F. 1992, ApJ, 396, 469

Hunter, J. D. 2007, CSE, 9, 90

Jaeger, T. R., Hyman, S. D., Kassim, N. E., & Lazio, T. J. W. 2012, AJ, 143, 96

Jenet, F. A., & Gil, J. 2003, ApJL, 596, L215

Kellermann, K. I., Kovalev, Y. Y., Lister, M. L., et al. 2007, Ap&SS, 311, 231

Kellermann, K. I., & Pauliny-Toth, I. I. K. 1969, ApJL, 155, L71

Koay, J. Y., Macquart, J. P., Jauncey, D. L., et al. 2018, MNRAS, 474, 4396

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

Lazio, T. J. W., Clarke, T. E., Lane, W. M., et al. 2010, AJ, 140, 1995

Levinson, A., Ofek, E. O., Waxman, E., & Gal-Yam, A. 2002, ApJ, 576, 923

Lovell, J. E. J., Jauncey, D. L., Bignall, H. E., et al. 2003, AJ, 126, 1699

Luber, N., van Gorkom, J. H., Hess, K. M., et al. 2019, AJ, 157, 254

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

Marchesi, S., Civano, F., Elvis, M., et al. 2016, ApJ, 817, 34

Marscher, A. P., & Gear, W. K. 1985, ApJ, 298, 114

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007a, in ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

Metzger, B. D., Williams, P. K. G., & Berger, E. 2015, ApJ, 806, 224

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

Mooley, K. P., Frail, D. A., Ofek, E. O., et al. 2013, ApJ, 768, 165

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

Murphy, E. J. 2009, ApJ, 706, 482

Murphy, E. J., Bolatto, A., Chatterjee, S., et al. 2018, in ASP Conf. Ser. 517, The ngVLA Science Case and Associated Science Requirements, ed. E. Murphy (San Francisco, CA: ASP), 3

Murphy, T., Chatterjee, S., Kaplan, D. L., et al. 2013, PASA, 30, e006

O'Dea, C. P. 1998, PASP, 110, 493

Ofek, E. O., Frail, D. A., Breslauer, B., et al. 2011, ApJ, 740, 65

Oliver, S. J., Bock, J., Altieri, B., et al. 2012, MNRAS, 424, 1614

Panessa, F., Baldi, R. D., Laor, A., et al. 2019, NatAs, 3, 387

Perley, R. 2016, EVLA Memo 195: Jansky Very Large Array Primary Beam Characteristics, https://library.nrao.edu/public/memos/evla/EVLAM_195.pdf

Perley, R. A., & Butler, B. J. 2013a, ApJS, 204, 19

Perley, R. A., & Butler, B. J. 2013b, ApJS, 206, 16

Perley, R. A., & Butler, B. J. 2017, ApJS, 230, 7

Perley, R. A., Chandler, C. J., Butler, B. J., & Wrobel, J. M. 2011, ApJL, 739, L1

Pietka, M., Fender, R. P., & Keane, E. F. 2015, MNRAS, 446, 3687

Plotkin, R. M., Miller-Jones, J. C. A., Chomiuk, L., et al. 2019, ApJ, 874, 13

Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2

Radcliffe, J. F., Beswick, R. J., Thomson, A. P., et al. 2019, MNRAS, 490, 4024

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

Readhead, A. C. S. 1994, ApJ, 426, 51

Rickett, B. J. 1990, ARA&A, 28, 561

Rowlinson, A., Bell, M. E., Murphy, T., et al. 2016, MNRAS, 458, 3506

Rowlinson, A., Stewart, A. J., Broderick, J. W., et al. 2019, A&C, 27, 111

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

Scargle, J. D., Norris, J. P., Jackson, B., & Chiang, J. 2013, ApJ, 764, 167

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

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

Shankar, F., Weinberg, D. H., & Miralda-Escudé, J. 2009, ApJ, 690, 20

Simonetti, J. H., Cordes, J. M., & Heeschen, D. S. 1985, ApJ, 296, 46

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

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

Smolčić, V., Delvecchio, I., Zamorani, G., et al. 2017a, A&A, 602, A2

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

Smolčić, V., Novak, M., Delvecchio, I., et al. 2017c, A&A, 602, A6

Stewart, A. J., Fender, R. P., Broderick, J. W., et al. 2016, MNRAS, 456, 2321

Swinbank, J. D., Staley, T. D., Molenaar, G. J., et al. 2015, A&C, 11, 25

Szokoly, G. P., Bergeron, J., Hasinger, G., et al. 2004, ApJS, 155, 271

Taylor, G. B., Momjian, E., Pihlström, Y., Ghosh, T., & Salter, C. 2005, ApJ, 622, 986

Thyagarajan, N., Helfand, D. J., White, R. L., & Becker, R. H. 2011, ApJ, 742, 49

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, CSE, 13, 22

van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Walker, M. A. 1998, MNRAS, 294, 307

Williams, P. K. G., & Berger, E. 2016, ApJ, 821, L22

Zensus, J. A. 1997, ARA&A, 35, 607