Quasars That Have Transitioned from Radio-quiet to Radio-loud on Decadal Timescales Revealed by VLASS and FIRST

2.1.1. FIRST

2.1.2. VLASS

VLASS is a wide-area (33,885 deg²), high-resolution (25), full-polarization, broadband radio continuum survey currently being carried out by NRAO at S band (2–4 GHz). The relative sky footprints of FIRST and VLASS are illustrated in Figure 1. In order to enable transient science, a key science goal of the survey, VLASS observations are being taken over a series of three epochs with a cadence of 32 months. Each VLASS epoch reaches a 1σ sensitivity of 0.12 mJy beam⁻¹, comparable to the depth of FIRST.

Zoom In Zoom Out Reset image size

VLASS Epoch 1 observations have recently been completed (2017–2019), and preliminary “QuickLook” images have been made publicly available on the NRAO archive.²⁶ The VLASS QuickLook images have flux uncertainties up to ∼20% and are not final survey data products (see Lacy et al. 2020 for details). Despite the limitations of the QuickLook images, they are valuable for identifying transient and variable sources that have experienced substantial changes in flux (larger than a factor of 2) in the time period between FIRST and VLASS.

Currently, no “official” source catalog for VLASS has been made publicly available. We therefore produced our own source catalog over a subset of VLASS Epoch 1.2 covering an area of ∼3440 deg². We used the source finding algorithm PyBDSF to compile a raw VLASS catalog based on the QuickLook images of sources brighter than 1 mJy. Although our goal is to identify sources that have recently brightened in the radio, we note that similar studies of reverse transients (i.e., sources that were previously detected in FIRST but discovered to have declined drastically in flux when reobserved in VLASS one to two decades later) are also currently in progress (e.g., Law et al. 2018; Cendes 2020).

We performed a positional cross-match (within a radius of 1″) with the FIRST catalog. The high prevalence of artifacts (uncleaned sidelobes) in the VLASS QuickLook images poses a challenge for reliable source extraction. We manually vetted our preliminary catalog of VLASS-FIRST transients by eye to remove spurious VLASS sources associated with image artifacts or diffuse emission (radio lobes) with inherently ill-defined positions. Such spurious sources account for the vast majority (∼80%–90%) of our preliminary transient catalog. Our manually vetted catalog covering the 3440 deg² of VLASS Epoch 1 analyzed thus far contains ∼2000 candidate transients that are compact and brighter than 1 mJy in VLASS but below the formal detection threshold (1 mJy) of FIRST. Additional details will be presented in D. Z. Dong et al. (2020, in preparation).

To select AGNs with variable radio emission, we considered both optical and infrared AGN selection techniques to capture both the unobscured and obscured populations. In the optical, we used the compilation of spectroscopically verified quasars from the Sloan Digital Sky Survey (SDSS; York et al. 2000) DR14 (Pâris et al. 2018) and found 52 candidate radio-variable quasars within <1″ of the VLASS position. Obscured AGNs that have recently brightened in the radio were identified on the basis of their infrared colors using data from the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010). We used the “R90” diagnostic criteria from Assef et al. (2018), which were designed to identify AGNs with 90% reliability. A total of 144 obscured AGNs with variable radio emission were identified using the R90 catalog, of which 29 sources also met our optical selection criteria, making a total of 167 AGNs associated with our manually vetted FIRST/VLASS sample.

With the AGN associations confirmed, we selected a subset of sources with peak flux densities brighter than 3 mJy beam⁻¹ in VLASS. This selects a more reliable sample of radio-variable AGNs, since the implied spectral index between FIRST and VLASS would be steeper than ν^2.5, thereby ruling out steady sources that are optically thick²⁷ owing to synchrotron self-absorption (SSA), though free–free absorption (FFA) is still a possibility. Following the 3 mJy cut, the final sample of radio-variable AGNs contains 26 objects, and it is these that we consider in this paper. We provide an illustration of our selection criteria in Figure 2.

Zoom In Zoom Out Reset image size

We summarize the basic properties of our sample in Tables 1 and 2. Of these 26 sources, 13 are selected from SDSS spectroscopy and the remaining 13 are obscured AGNs selected based on the WISE criteria described in Section 2.2. The 13 optically selected AGNs are classified as broad-line quasars with redshifts in the range 0.2 < z < 3.2. Shen et al. (2011) report bolometric luminosities of $\mathrm{log}({L}_{\mathrm{bol}}/\mathrm{erg}\,{{\rm{s}}}^{-1})\approx 45.2\mbox{--}46.8$ from SDSS spectroscopy, virial SMBH masses in the range $\mathrm{log}(M/{M}_{\odot })\approx 8\mbox{--}9.7$, and Eddington ratios of 6%–20% (consistent with radiatively efficient SMBH accretion; e.g., Heckman & Best 2014). All of our sources would have previously been classified as radio-quiet based on their upper limits in FIRST. The observed increase in radio flux between FIRST and VLASS suggests that they are no longer in a radio-quiet state. Given the high VLASS radio luminosities (${L}_{3\mathrm{GHz}}={10}^{40\mbox{--}42}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$), our sources are now consistent with radio-loud quasars (Kellermann et al. 2016).

Table 1.  Sample

Source	R.A.	Decl.	z	LBol	log(MSMBH)	${L}_{\mathrm{Bol}}/{L}_{\mathrm{Edd}}$	Type
	(J2000)	(J2000)		(erg s−1)	(${M}_{\odot }$)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
J0742+2704	115.701796	27.070113	0.6264	45.57	8.67	−1.20	SDSS, WISE
J0751+3154	117.877557	31.904037	1.8640	46.57	9.47	−1.00	SDSS, WISE
J0800+3124	120.044683	31.415871	1.9368	46.93	9.56	−0.73	SDSS, WISE
J0807+2102	121.767970	21.035786	1.5588	46.31	9.34	−1.13	SDSS, WISE
J0832+2302	128.198513	23.042816	0.9430	⋯	⋯	⋯	SDSS, WISE
J0950+5128	147.653173	51.477212	0.2142	45.18	7.98	−0.90	SDSS
J1023+2219	155.842720	22.321296	⋯	⋯	⋯	⋯	WISE
J1037−0736	159.491631	−7.607362	⋯	⋯	⋯	⋯	WISE
J1112+2809	168.055854	28.165137	⋯	⋯	⋯	⋯	WISE
J1157+3142	179.388528	31.707514	0.8910	⋯	⋯	⋯	SDSS, WISE
J1204+1918	181.218737	19.306199	2.3440	⋯	⋯	⋯	SDSS, WISE
J1208+4741	182.241352	47.698944	1.0915	⋯	⋯	⋯	SDSS, WISE
J1246+1805	191.518270	18.089036	⋯	⋯	⋯	⋯	WISE
J1254−0606	193.521438	−6.109209	⋯	⋯	⋯	⋯	WISE
J1333−0349	203.334269	−3.832307	⋯	⋯	⋯	⋯	WISE
J1347+4505	206.836689	45.098706	⋯	⋯	⋯	⋯	WISE
J1357−0329	209.443830	−3.484817	⋯	⋯	⋯	⋯	WISE
J1413+0257	213.454848	2.953648	⋯	⋯	⋯	⋯	WISE
J1447+0512	221.817580	5.207433	1.7475	46.36	9.02	−0.76	SDSS, WISE
J1507−0549	226.962335	−5.819718	⋯	⋯	⋯	⋯	WISE
J1512−0533	228.096364	−5.550100	⋯	⋯	⋯	⋯	WISE
J1514+4000	228.520730	40.013724	2.1226	46.80	9.74	−1.04	SDSS
J1546+1500	236.641541	15.010193	⋯	⋯	⋯	⋯	WISE
J1603+1809	240.828155	18.151432	3.2460	⋯	⋯	⋯	SDSS
J1609+4306	242.441175	43.106462	⋯	⋯	⋯	⋯	WISE
J2109−0644	317.320264	−6.743461	1.0812	46.38	9.05	−0.77	SDSS

Note. Column (1): source name. Columns (2) and (3): source R.A. and decl. Column (4): spectroscopic redshift from SDSS DR14 from Pâris et al. (2018). Column (5): bolometric quasar luminosity. Column (6): virial SMBH mass estimate from SDSS spectroscopy from Shen et al. (2011). Column (7): Eddington ratio. Column (8): AGN type. Sources are classified as AGNs on the basis of optical spectroscopy from SDSS (Pâris et al. 2018) and/or infrared colors based on data from WISE (Assef et al. 2018).

Download table as:  ASCIITypeset image

Table 2.  FIRST and VLASS Properties

Source	DateFIRST	σFIRST	FFIRST	DateVLASS	σVLASS	FVLASS	$\mathrm{log}({L}_{\mathrm{VLASS}})$
		(mJy beam−1)	(mJy beam−1)		(mJy beam−1)	(mJy beam−1)	(erg s−1)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
J0742+2704	1995 Nov 10	0.166	0.53	2019 Apr 10	0.146	9.19	41.68
J0751+3154	1995 Oct 23	0.148	<0.44	2019 Apr 13	0.148	3.00	42.36
J0800+3124	1994 Jun 4	0.123	0.38	2019 Apr 13	0.125	4.54	42.58
J0807+2102	1998 Sep	0.157	<0.47	2019 Apr 14	0.137	3.69	42.26
J0832+2302	1995 Dec 28	0.157	0.48	2019 Apr 13	0.130	3.08	41.64
J0950+5128	1997 Apr 29	0.148	<0.44	2019 Apr 18	0.126	8.77	40.57
J1023+2219	1996 Jan	0.205	0.65	2019 Apr 18	0.145	3.50	⋯
J1037−0736	2002 Jun	0.136	<0.41	2019 Apr 16	0.151	12.95	⋯
J1112+2809	1995 Nov 4	0.156	0.49	2019 Apr 13	0.120	3.19	⋯
J1157+3142	1994 Jun 6	0.601	2.02	2019 Apr 14	0.135	8.58	42.03
J1204+1918	1999 Nov	0.149	0.48	2019 Mar 21	0.145	5.09	42.84
J1208+4741	1997 Apr 5	0.157	<0.47	2019 Apr 15	0.117	3.20	41.82
J1246+1805	1999 Nov 24	0.147	<0.44	2019 Apr 12	0.122	7.09	⋯
J1254−0606	2001 Apr	0.144	0.46	2019 Mar 7	0.247	8.76	⋯
J1333−0349	1998 Sep	0.14	<0.42	2019 Mar 24	0.175	3.09	⋯
J1347+4505	1997 Mar 3	0.151	<0.45	2019 Mar 19	0.129	3.58	⋯
J1357−0329	1998 Sep	0.156	0.49	2019 Mar 24	0.170	3.11	⋯
J1413+0257	1998 Jul	0.142	0.47	2019 Apr 2	0.229	6.91	⋯
J1447+0512	2000 Feb	0.16	0.49	2019 Mar 12	0.152	4.30	42.45
J1507−0549	2001 Apr	0.139	<0.42	2019 Mar 7	0.159	4.73	⋯
J1512−0533	2001 Apr	0.233	<0.70	2019 Mar 7	0.200	3.25	⋯
J1514+4000	1994 Aug 21	0.141	<0.42	2019 Mar 28	0.164	3.49	42.57
J1546+1500	1999 Dec	0.144	<0.43	2019 Apr 11	0.122	4.69	⋯
J1603+1809	1999 Nov	0.148	<0.44	2019 Apr 12	0.116	3.61	43.03
J1609+4306	1994 Jul 24	0.124	0.45	2019 May 4	0.126	4.53	⋯
J2109−0644	1997 Feb 28	0.153	<0.46	2019 Apr 5	0.170	18.30	42.57

Note. Column (1): source name. Column (2): date of FIRST observation. For some sources, only the year and month are provided in the header of the FIRST image. Column (3): local 1σ rms noise in FIRST. Column (4): peak flux density at 1.4 GHz from FIRST. Newly identified sources in VLASS that have faint FIRST counterparts with peak flux densities >3σ_FIRST are reported as detections. All other sources are reported as FIRST upper limits. We note that none of the sources in this table meet the selection criteria of the FIRST catalog (Helfand et al. 2015). Column (5): date of VLASS observation. Column (6): local 1σ rms noise in VLASS. Column (7): peak flux density at 3 GHz from VLASS. Column (8): radio luminosity at 3 GHz from VLASS.

Download table as:  ASCIITypeset image

In Figure 3, we show image cut-outs from FIRST and VLASS. The 1σ rms sensitivities of the FIRST and VLASS images are typically 0.15 and 0.12 mJy beam⁻¹, respectively. The local rms noise level measured in each image is given in Table 2. The quality of the FIRST images is generally quite good. However, we note that one source in particular (J1157+3142) appears to have unusually poor image quality in the FIRST survey, making its classification as a variable AGN uncertain. We emphasize that while none of our sources met the formal detection threshold criteria (F_peak > 1 mJy) of the final release of the FIRST source catalog (Helfand et al. 2015), 12 sources have faint FIRST emission at the location of the VLASS source at the 3σ–3.6σ level (J1204+1918). We list sources with peak FIRST flux densities >3σ as FIRST “detections” in Table 2.

Zoom In Zoom Out Reset image size

The VLASS QuickLook image cut-outs shown in Figure 3 highlight the radio variability of our sources on timescales of roughly one to two decades, as well as their compactness on roughly arcsecond scales. Some of the VLASS QuickLook images suffer from strong imaging artifacts, such as prominent sidelobes due to insufficient cleaning and/or the presence of residual phase errors. This is due to the rough nature of the QuickLook images,²⁸ which we emphasize are not the final VLASS survey products. Despite the known limitations of the QuickLook images, we find the peak VLASS flux densities of our compact sources to typically be consistent with our independent VLA follow-up imaging at S band (see Section 4.2 for a discussion of source variability at S band). Thus, our study supports the viability of the use of VLASS QuickLook images for science as long as users take the known limitations into account.

We calibrated our data using the scripted VLA calibration pipeline²⁹ for the Common Astronomy Software Applications (CASA) package (McMullin et al. 2007) version 5.3.0. In order to ensure the use of an optimal calibration solution interval for each band, the pipeline was run separately for each band of each SB.

Imaging and self-calibration were performed in CASA 5.6.0 following standard heuristics for wide-field, broadband VLA data. We used the TCLEAN task to produce a full field-of-view image at each band and employed the w-projection algorithm (Cornwell et al. 2005) to correct for non-coplanar baselines. To avoid bandwidth smearing over the large fractional bandwidths of our observations, deconvolution was performed using the MTMFS (multiterm, multifrequency synthesis) algorithm (Rau & Cornwell 2011). We adopted a Briggs weighting scheme (Briggs 1995) with a robust parameter between −0.5 and 0.5 to achieve the best compromise among sidelobe levels, resolution, and sensitivity given the quality of the uv-coverage and resulting point-spread function (PSF) of each data set.³⁰

Commensal low-frequency (<1 GHz) data were also recorded during our VLA follow-up campaign. This commensal system, known as the VLA Low-band Ionosphere and Transient Experiment (VLITE; Clarke et al. 2016; Polisensky et al. 2016), records data at 340 MHz simultaneously with regular VLA observing programs.³¹ During our VLA observations, commensal VLITE data with a maximum of 18 antennas were recorded.

VLITE data are automatically calibrated and imaged through an analysis pipeline based on the Obit (Cotton 2008) data reduction package. The VLITE pipeline images from our snapshot observations typically achieve a 1σ depth of 3 mJy beam ⁻¹ and a maximum angular resolution of ≈5″. We visually inspected all VLITE imaging pipeline products and found that none of our sources were detected, as expected given their curved or flat radio SEDs. We discuss constraints on the radio SEDs and underlying emission mechanisms using the upper limits from VLITE in Section 4.3.

Source flux and shape measurements for the VLA data above 1 GHz were made using the CASA IMFIT task and are reported in Table A1. For the majority of sources, the flux errors include a standard 3% uncertainty in the absolute flux scale (Perley & Butler 2017), added in quadrature with the flux error reported by IMFIT. Due to scheduling limitations, one of our sources (J2109−0644) was observed using 3C 48 as the flux calibrator. Because 3C 48 is currently flaring, we conservatively report a flux uncertainty of 20% for each measurement of J2109−0644 following current NRAO guidelines.³²

Our sources are characterized by compact (≲01) radio continuum emission over the full frequency range of our VLA study, which provides a maximum angular resolution of θ_max ≲ 01 (for the Ku-band/A-configuration data). This corresponds to an upper limit on the intrinsic linear source sizes that is <1 kpc.

We illustrate the flux variations of our sample at L and S band over two epochs probing two different timescales in Figures 4 and 5. We also show the variability at X band for the three sources (J0807+2102, J0832+2302, and J0950+5128) in our sample with a gap of several weeks between observations from 1–12 GHz (L, S, C, and X band) and 12–18 GHz (Ku band). Tables 2 and 3 provide additional details on the observing dates of the data from FIRST, VLASS, and our new multiband JVLA follow-up observations.

Zoom In Zoom Out Reset image size

At L band, the two epochs shown in the top panel of Figure 4 are from the FIRST survey (with observing dates ranging from 1995 to 2002 for our sample) and our 2019 JVLA follow-up study. The L-band source flux densities from our new JVLA observations range from ∼1 to 11 mJy. Thus, compared to the FIRST flux measurements made 17–24 yr earlier, our sources have L-band flux densities that have increased dramatically by factors of ∼2–25. The left panel of Figure 5 shows the variability amplitude at L band between FIRST and our 2019 observations for each source. The weighted median value is ∼200%. However, we emphasize that since half of our sources remain undetected below the 3σ level in FIRST, this is likely to be a lower limit. Thus, the true L-band variability amplitude on timescales of roughly one to two decades could be even higher.

Zoom In Zoom Out Reset image size

At S band, the two observing epochs shown in Figure 4 span a much narrower range of time (∼months) compared to the variability timescale constraints at L band (∼decades). The S-band flux densities between our JVLA follow-up observations and VLASS are approximately constant over timescales of 3–7 months. This is clearly illustrated in the right panel of Figure 5, which shows that the S-band variability amplitude has a median value of ∼15%. Thus, our follow-up JVLA S-band data are typically in good agreement within the current ∼20% flux uncertainties of VLASS. We note that the most substantial outlier in the right panel of Figure 5 is the source J2109−0644, which, as discussed in Section 3, may have unusually large absolute flux errors due to the variability of the flux calibrator.

At X band, the two observing epochs shown in Figure 4 for J0807+2102, J0832+2302, and J0950+5128 span an even narrower range of time (∼weeks) compared to the variability timescale constraints at L and S band. The X-band flux densities are roughly constant, although all three sources exhibit a slight increase in flux of ∼10% over a period of 58 days. Although an increase in flux density of this magnitude could be due to the evolution of a nascent jet (see Section 5), the presence of residual errors in the absolute flux scale of similar magnitude cannot be entirely ruled out. This is because the first epoch of X-band observations was taken during the move from the hybrid BnA configuration to the standard A array, which led to poor PSF quality that ultimately posed challenges for self-calibration and deconvolution (see Section 3). Nevertheless, we conclude that the X-band flux densities are in close enough agreement to justify the inclusion of the observations above 12 GHz in our radio SED analysis.

Although our time domain assessment is currently crude and would benefit from higher-cadence monitoring in the future, we conclude that our sources are characterized by large-amplitude (2–25 times) variability at L band on timescales of one to two decades but maintain roughly constant flux densities over timescales of a few months at S band. Thus, the typical variability timescale, τ_var, likely falls in the range of 3–7 months < τ_var < 17–24 yr. We discuss implications of these constraints on the origin of the radio variability in Section 5.

In Figure 6, we show the radio SEDs of the 14 sources from our sample with quasi-simultaneous VLA observations from 1 to 18 GHz. The radio SEDs have peaked spectral shapes with varying degrees of curvature easily discernible by eye. This indicates the presence of an underlying absorption mechanism, such as SSA or FFA, as is commonly associated with compact radio sources (Callingham et al. 2017; Collier et al. 2018).

Zoom In Zoom Out Reset image size

4.3.1. Radio SED Modeling

To quantify the location of the spectral peak and degree of curvature, we performed least-squares fits to the quasi-simultaneous VLA data above 1 GHz using two basic synchrotron emission models:

• 1.  
  A standard nonthermal power-law model defined by S_ν = a ν^α, where S_ν is the flux at frequency ν, a represents the amplitude, and α is the spectral index.
• 2.  
  A curved power-law model defined by ${S}_{\nu }=a{\nu }^{\alpha }{e}^{q{(\mathrm{ln}\nu )}^{2}}$, where q represents the degree of spectral curvature (the breadth of the peak of the radio SED). The q parameter is defined by ${\nu }_{\mathrm{peak}}={e}^{-\alpha /2q}$, where ν_peak is the turnover frequency. Significant spectral curvature is typically defined as $| q| \geqslant 0.2$ (Duffy & Blundell 2012; Callingham et al. 2017).

We note that our use of the curved power-law model is phenomenological in nature. More physically motivated spectral curvature models, such as SSA or FFA (or models with contributions from multiple electron populations; Tingay et al. 2015), require more than three free parameters, and thus have only a few degrees of freedom given the number of bands (typically 5; see Table 3) available for each source.

We summarize our spectral modeling analysis in Table 4. Based on the reduced chi-square values provided in this table, a curved power-law model provides a better fit compared to a simple power-law model for all sources. The observed spectral turnover frequencies for our sources range from 2.5 to 22.7 GHz, with a median value of ν_peak = 6.6 GHz. For sources with measured redshifts, the corresponding rest-frame turnover frequencies lie in the range of ${\nu }_{\mathrm{peak},\mathrm{rest}}=6.8\mbox{--}58.1\,\mathrm{GHz}$. Likewise, the values of q from the curved power-law fits span the range $0.18\lt | q| \lt 1.09$, confirming that essentially all are strongly curved. The exception is J2109−0644, which is only marginally below the formal limit of $| q| \gt 0.2$, and the VLITE upper limit at 340 MHz supports a curved, not flat, radio SED.

Table 4.  Radio Spectral Modeling Parameters

Source	θmax	Linear Size	νpeak	${\nu }_{\mathrm{peak},\mathrm{rest}}$	${\chi }_{\mathrm{red},\mathrm{PL}}^{2}$	${\chi }_{\mathrm{red},\mathrm{CPL}}^{2}$	q	${\alpha }_{\mathrm{thick}}$	αthin
	(arcsec)	(pc)	(GHz)	(GHz)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)
J0742+2704	<0.17	<1193	8.6	14.0	107.75	9.23	−0.65 ± 0.09	1.59 ± 0.10	−0.62 ± 0.17
J0807+2102	<0.11	<956	21.9	56.2	67.53	35.97	−0.36 ± 0.16	1.25 ± 0.11	⋯
J0832+2302	<0.17	<1379	11.4	22.1	39.49	3.96	−0.40 ± 0.06	1.36 ± 0.11	⋯
J0950+5128	<0.16	<575	13.2	16.0	34.53	5.78	−0.34 ± 0.07	1.43 ± 0.10	⋯
J1037−0736	<0.17	⋯	6.2	⋯	74.88	5.74	−0.54 ± 0.07	1.35 ± 0.11	−0.44 ± 0.17
J1204+1918	<0.16	<1341	4.6	15.5	132.82	3.99	−0.75 ± 0.06	1.27 ± 0.10	−1.16 ± 0.17
J1208+4741	<0.12	<1006	5.6	11.8	89.17	6.44	−0.58 ± 0.08	1.32 ± 0.11	−0.61 ± 0.16
J1246+1805	<0.13	⋯	9.0	⋯	46.30	2.52	−0.40 ± 0.05	1.28 ± 0.10	−0.09 ± 0.17
J1254−0606	<0.15	⋯	6.6	⋯	236.24	2.25	−1.09 ± 0.04	2.57 ± 0.10	−1.16 ± 0.17
J1413+0257	<0.21	⋯	6.4	⋯	215.22	1.96	−1.02 ± 0.04	2.20 ± 0.10	−1.22 ± 0.17
J1447+0512	<0.18	<1562	2.6	7.3	43.75	29.84	−0.34 ± 0.17	⋯	−0.45 ± 0.17
J1546+1500	<0.16	⋯	4.8	⋯	43.58	11.91	−0.38 ± 0.10	0.84 ± 0.10	−0.29 ± 0.17
J1603+1809	<0.17	<1308	8.1	34.6	109.40	0.23	−0.66 ± 0.01	1.76 ± 0.10	−0.47 ± 0.17
J2109−0644	<0.17	<1423	6.5	13.5	13.74	9.94	−0.18 ± 0.09	0.15 ± 0.11	−0.41 ± 0.16

Note. Column (1): source name. Column (2): angular size upper limit from the highest-resolution VLA data available. Column (3): linear size upper limit for sources with known redshifts. Column (4): spectral peak (or “turnover”) frequency (ν_peak) from a fit to a curved power-law model. Column (5): same as Column (4), but in the rest frame for sources with known redshifts. Column (6): reduced chi squared value for the fit to a simple power-law model. Column (7): reduced chi squared value for the fit to a curved power-law model. Column (8): spectral curvature parameter, q, from a fit to a curved power-law model. Column (9): optically thick spectral index, α_thick, estimated by a power-law fit to the two lowest-frequency VLA bands (L and S band). The uncertainties were calculated using standard propagation of errors. We required ν_peak > 4 GHz to estimate α_thick. Column (10): estimates of the optically thin spectral index, α_thin (and the corresponding uncertainties), from a power-law fit to the two highest-frequency VLA bands (either X and Ku bands or Ku and K bands) above ν_peak. We required quasi-simultaneous VLA measurements in at least two bands above ν_peak to estimate α_thin.

Download table as:  ASCIITypeset image

In Table 4, we also provide the optically thick and optically thin spectral index values (α_thick and α_thin) below and above the turnover frequency, respectively. We estimated α_thick by fitting a power-law model to the quasi-simultaneous flux densities at the two lowest-frequency VLA bands (L and S band) below the turnover frequency. Estimates for α_thick are only provided for sources with ν_peak > 4 GHz (spectral peaks above the VLA S band). For α_thin, we performed a power-law fit to the quasi-simultaneous flux densities at the two highest-frequency VLA bands (either X and Ku bands or Ku and K bands)³³ above the turnover frequency. We required the quasi-simultaneous VLA data to include at least two independent bands above the turnover frequency to estimate α_thin. The uncertainties in α_thick and α_thin provided in Table 4 were calculated using standard propagation of errors.

Based on the radio spectral modeling analysis described here, we conclude that the radio SEDs of our sources are best represented by curved power-law models (as opposed to flat power-law models lacking curvature), thus supporting the presence of significant spectral curvature. In a future study, we will incorporate new data from forthcoming VLA and Giant Metre-wave Radio Telescope (GMRT) observations spanning a broader frequency range (the full contiguous frequency range of the VLA from 1 to 50 GHz and deep measurements below 1 GHz using the VLA and GMRT). This will enable a more rigorous radio SED modeling analysis that will allow us to test whether SSA or FFA is responsible for the observed spectral curvature (Mhaskey et al. 2019a, 2019b).

4.3.2. Radio Color–Color Diagram

In Figure 7, we show our sources on a radio color–color diagram (α_thick vs. α_thin). All of our sources reside in the second quadrant of this figure, which characterizes peaked radio spectral shapes indicative of an underlying absorption mechanism commonly associated with source compactness owing to youth and/or confinement by an external medium (O’Dea 1998; Orienti 2016; O’Dea & Saikia 2020). Most of our sources also meet the more strict selection criteria for peaked-spectrum radio sources from Callingham et al. (2017) of α_thick > 0.1 and α_thin < −0.5.

Zoom In Zoom Out Reset image size

Our sources have optically thick spectral index constraints that are consistent with SSA, though improved spectral coverage at low frequencies (<1 GHz) will ultimately be needed (Callingham et al. 2015; Collier et al. 2018; Keim et al. 2019; Mhaskey et al. 2019a, 2019b). The identification of FFA associated with any of the newly radio-loud AGNs in our sample would be of particular interest in the context of galaxy evolution since simulations have shown that it may arise from jet−ISM interactions (Bicknell et al. 1997, 2018). Such jet−ISM interactions may subsequently influence the ambient star formation rate and/or efficiency, as argued by recent observational studies in low-redshift galaxies (Morganti et al. 2013; Nyland et al. 2013; Mukherjee et al. 2018; Husemann et al. 2019; Jarvis et al. 2019; Ruffa et al. 2019).

4.3.3. Infrared Color–Color Diagram

In Figure 8, we show the distribution of the 14 sources from our multiband VLA follow-up campaign in WISE infrared color space (W2 − W3 vs. W1 − W2). A summary of the infrared properties of our sources is provided in Table 5. All but one these sources meet the infrared color selection criteria for obscured quasars defined by Mateos et al. (2012). The single exception is J1603+1809, which is identified as a type 1 quasar in SDSS but does not meet the WISE AGN selection criteria of the Assef et al. (2018) R90 catalog.

Zoom In Zoom Out Reset image size

Table 5.  Infrared Source Properties

Source	W1	W2	W3	W4	$\mathrm{log}({L}_{6\mu {\rm{m}},\mathrm{rest}})$
	(mag)	(mag)	(mag)	(mag)	(erg s−1)
(1)	(2)	(3)	(4)	(5)	(6)
J0742+2704	14.42 ± 0.03	13.68 ± 0.04	11.25 ± 0.15	⋯	45.05
J0751+3154	15.98 ± 0.06	14.86 ± 0.07	⋯	⋯	46.24
J0800+3124	15.91 ± 0.05	14.70 ± 0.07	12.26 ± 0.41	⋯	46.30
J0807+2102	15.91 ± 0.06	14.72 ± 0.07	11.64 ± 0.25	⋯	46.14
J0832+2302	13.84 ± 0.03	12.46 ± 0.03	9.45 ± 0.04	7.21 ± 0.11	46.07
J0950+5128	13.43 ± 0.02	12.80 ± 0.03	10.40 ± 0.08	8.22 ± 0.25	44.13
J1023+2219	15.85 ± 0.05	14.77 ± 0.07	11.32 ± 0.20	8.47 ± 0.34	⋯
J1037−0736	14.82 ± 0.03	13.83 ± 0.04	11.22 ± 0.16	⋯	⋯
J1112+2809	16.33 ± 0.07	15.40 ± 0.10	⋯	⋯	⋯
J1157+3142	14.98 ± 0.03	13.82 ± 0.04	11.26 ± 0.16	8.62 ± 0.41	45.40
J1204+1918	16.07 ± 0.06	15.17 ± 0.09	12.08 ± 0.33	8.85 ± 0.41	46.54
J1208+4741	16.16 ± 0.06	15.18 ± 0.07	12.70 ± 0.43	⋯	45.26
J1246+1805	16.53 ± 0.07	15.43 ± 0.10	12.72 ± 0.47	⋯	⋯
J1254−0606	15.92 ± 0.05	14.65 ± 0.06	10.35 ± 0.06	7.83 ± 0.16	⋯
J1333−0349	15.34 ± 0.04	14.35 ± 0.05	11.99 ± 0.25	⋯	⋯
J1347+4505	15.22 ± 0.04	13.87 ± 0.03	10.82 ± 0.09	⋯	⋯
J1357−0329	13.87 ± 0.03	12.77 ± 0.03	9.87 ± 0.05	7.39 ± 0.10	⋯
J1413+0257	18.14 ± 0.25	16.11 ± 0.17	11.65 ± 0.17	8.94 ± 0.32	⋯
J1447+0512	16.20 ± 0.06	15.14 ± 0.07	⋯	⋯	46.03
J1507−0549	15.97 ± 0.06	15.12 ± 0.10	12.04 ± 0.29	⋯	⋯
J1512−0533	16.66 ± 0.08	15.19 ± 0.10	12.43 ± 0.39	9.10 ± 0.47	⋯
J1514+4000	16.22 ± 0.05	15.54 ± 0.08	12.13 ± 0.19	9.53 ± 0.52	46.21
J1546+1500	15.97 ± 0.05	14.90 ± 0.07	⋯	⋯	⋯
J1603+1809	16.10 ± 0.05	15.47 ± 0.09	12.59 ± 0.48	⋯	47.05
J1609+4306	16.38 ± 0.05	15.41 ± 0.08	13.12 ± 0.43	9.62 ± 0.48	⋯
J2109−0644	15.14 ± 0.04	14.07 ± 0.04	11.21 ± 0.17	8.22 ± 0.31	45.71

Note. Column (1): source name. Columns (2)–(5): WISE W1 (3.4 μm), W2 (4.6 μm), W3 (12 μm), and W4 (22 μm) band magnitudes and errors from the AllWISE source catalog (Cutri et al. 2013). Column (6): estimated rest-frame 6 μm luminosity extrapolated from a power-law fit to the AllWISE photometry.

Download table as:  ASCIITypeset image

We also highlight the region occupied by extremely red and powerful AGNs defined by (W1 − W2) + 1.25(W2 − W3) > 7 (Lonsdale et al. 2015) in Figure 8. Sources with such extreme mid-infrared colors are believed to be heavily obscured, ultraluminous quasars. Of our 14 sources, one source, J1413+0257, meets the Lonsdale et al. (2015) selection criteria. Although this source currently lacks redshift information, Patil et al. (2020) reported redshifts in the range 0.4 < z < 3.0 (with a median value of z ∼ 1.5) and luminosities of $\mathrm{log}({L}_{\mathrm{bol}}/{L}_{\odot })\sim 11.7\mbox{--}14.2$ for a sample of quasars in the Lonsdale et al. (2015) color space also harboring compact radio sources.

Finally, we use Figure 8 to investigate possible relationships between the WISE colors and two key parameters from our radio spectral modeling: the spectral curvature parameter, q, and the optically thick spectral index, α_thick. Qualitatively, Figure 8 suggests an association between both a higher degree of spectral curvature (a narrower spectral peak) and α_thick values for redder sources.

Previous studies have reported evidence for a connection between quasar reddening and radio jet properties possibly arising from hierarchical SMBH−galaxy coevolution (Georgakakis et al. 2012; Glikman et al. 2012; Klindt et al. 2019; Fawcett et al. 2020; Patil et al. 2020; Rosario et al. 2020). In such a scenario, a quasar may transition to a radio-loud phase following a gas-rich merger and subsequent change in SMBH accretion state and/or spin conducive to jet formation. In a future study, we will investigate this possibility in more detail by constraining the origin of the reddening in our sample through optical and infrared SED modeling. Ultimately, observations of a larger sample will be needed to more firmly establish the relationship between the infrared colors and radio SED properties in young radio quasars.

Radio waves traveling through an inhomogeneous plasma may be modulated by scattering or lensing phenomena. This may lead to variations in the observed properties of compact sources, such as flat-spectrum radio AGNs. Physical mechanisms responsible for such propagation effects include interplanetary scintillation (IPS; Morgan et al. 2018, and references therein), interstellar scintillation (ISS; Jauncey et al. 2016, and references therein), and extreme scattering events (ESEs; Fiedler et al. 1987).

IPS arises from turbulence in the solar wind plasma and leads to rapid flux variability on timescales of roughly seconds in low-frequency (∼100 MHz) radio observations. Although IPS has proven to be a powerful tool for identifying compact radio AGNs in low-frequency surveys (e.g., Chhetri et al. 2018), it is not compatible with the variability timescales and amplitudes observed in our sample above 1 GHz.

ISS is caused by the refraction or diffraction of radio waves emitted by a microarcsecond-scale source as they pass through fluctuations in the plasma density and/or magnetic field of the ionized ISM in our Galaxy (Stanimirović & Zweibel 2018). This leads to variable radio emission on timescales of hours to days with a strong dependence on galactic latitude.

Observational manifestations of ISS (i.e., variability amplitude and timescale) depend on the Galactic free electron density along a particular line of sight, as well as the observing frequency. The distribution of our sources in Galactic latitude is approximately flat (Figure 1), arguing against ISS as the main driver of the observed radio variability. As a more quantitative test of an ISS variability origin for each source, we computed the critical frequency,³⁴ ν₀, below which the modulation due to scintillation is in the strong regime.³⁵ At the positions of our sources we find ν₀ = 7.1–11 GHz. For observations at ν = 1.5 GHz, we are therefore well within the strong scattering regime and the equations from Section 3.2 of Walker (1998) apply. For point sources in this limit that are experiencing refractive scintillation, the expected flux modulation at 1.5 GHz is m = (ν₀/ν)^17/30 ∼ 30%–40%. This modulation is expected to occur on a timescale of ${t}_{r}\sim 2{({\nu }_{0}/\nu )}^{11/5}$ hr ∼ 2–7 days. This stands in sharp contrast to the observed modulations of ≳100% to ≳2500% occurring on observed timescales of months to decades for our sources, thus ruling out refractive ISS as the origin of the radio variability in our sample.

Diffractive ISS, which is associated with interference, leads to variability over an even narrower bandwidth and shorter timescale, making it considerably less plausible than refractive ISS. The fluctuations (of order unity) last for only ∼0.1–0.3 hr. Furthermore, at 1.5 GHz, any diffractive ISS would be decorrelated over a bandwidth of $\delta \nu \sim 1\mbox{--}5\,\mathrm{MHz}$. This modulation would thus be washed out by averaging over a single VLA frequency channel.

Another plasma lensing phenomenon known to produce radio variability is that of ESEs (Fiedler et al. 1987). In this case, refractive defocusing by the transverse passage of a discrete, high-density plasma lens in front of a compact radio source leads to a reduction of flux with a characteristic U-shaped light curve (Clegg et al. 1998; Bannister et al. 2016). The variability amplitude (a reduction in emission of ≲50%) and timescale (weeks to months) for an ESE are typically larger in magnitude compared to ISS. At least one case of an ESE with a threefold modulation in flux in the centimeter-wave radio regime is known (Bannister et al. 2016). This overlaps with the lower range in variability amplitude of our sample. However, we note that the timescale for the high-amplitude radio variability reported in Bannister et al. (2016) is a few months. This stands in sharp contrast to our sources, which vary on decades-long timescales but are steady over periods of a few months.

While we cannot totally rule out contributions to the observed radio variability by serendipitous propagation effects, we find such scenarios to be unlikely given inconsistencies with the high variability amplitudes and long timescales observed in our sample.

Variable radio emission may be associated with supernovae and gamma-ray bursts (GRBs), arising from collimated outflows of high-mass supernovae and compact object mergers (neutron stars and stellar-mass black holes; Weiler et al. 2002; Woosley & Bloom 2006; Berger 2014). However, the high luminosities of our sources (Table 2) rule out the possibility of a radio supernova afterglow (Palliyaguru et al. 2019) as a progenitor scenario. Regarding the possibility of radio variability associated with a GRB origin, we note that our sources are on par with, or more luminous than, on-axis GRBs, which have peak luminosities of ∼10³¹ erg s⁻¹ Hz⁻¹ (Chandra & Frail 2012). However, the typical variability timescale of radio emission associated with GRBs is around 1–2 weeks (Pietka et al. 2015). This is considerably shorter than the variability timescale constraints for our sources (greater than a few months), thus ruling out the possibility of a GRB progenitor.

Tidal disruption events (TDEs; e.g., Komossa 2015, and references therein) occur when a star passes within the tidal radius of an SMBH and is ripped apart and partially accreted onto the SMBH. This leads to a multiwavelength flare, which in some cases includes the production of radio continuum emission associated with nonthermal jets or thermal outflows (van Velzen et al. 2011, 2016; Anderson et al. 2020). The most radio luminous known TDE, Swift J1644+57, peaked at ∼10³²  erg s⁻¹ Hz⁻¹ (Eftekhari et al. 2018), which is roughly in line with the radio luminosities of our sources. In addition, recent studies have suggested that TDEs may be responsible for triggering a significant fraction of the changing-look AGN population, particularly at z > 2 (Padmanabhan & Loeb 2020). Because more massive black holes have weaker tidal fields near their event horizons, TDEs for main-sequence stars become increasingly rare above a critical mass of about 10⁸ M_⊙ (Hills 1975). Main-sequence stars that approach SMBHs above this limiting mass are “swallowed whole” (MacLeod et al. 2012) and thus do not lead to the production of an electromagnetic flare.³⁶ As shown in Table 1, the available SMBH mass estimates of our sources typically exceed 10⁸ M_⊙, suggesting that TDEs are unlikely progenitors for the majority of our sources.

On the other hand, evolved stars with large diffuse envelopes, such as red giants, are in principle susceptible to tidal disruption by an SMBH with M_SMBH > 10⁸ M_⊙. However, such events are expected to be rare owing to the relatively short duration of the red giant phase compared to main-sequence lifetimes. Recent studies have also argued that TDEs of giant stars by SMBHs with masses above 10⁸ M_⊙ may be less luminous than TDEs associated with lower-mass SMBHs (Bonnerot et al. 2016).

Although TDEs may provide a plausible mechanism for driving the extreme radio variability in some of our sources, the current literature consensus is that TDEs with powerful relativistic jets are rare (e.g., van Velzen et al. 2018), composing only a few percent of the known TDE population³⁷ (Alexander et al. 2020, and references therein). Thus, based on the high radio luminosities and SMBH masses of our sample, as well as the low space density of radio-loud TDEs, we conclude that the radio variability in our sources is most likely not associated with TDEs.

While the large-scale lobes of radio galaxies are believed to remain steady over megayear to gigayear timescales (Blandford et al. 2019), intrinsic radio variability on human timescales (days to years) is common among AGNs with compact (<1 kpc) jets such as blazars (Lister 2001), unbeamed radio quasars/galaxies with young jets (Torniainen et al. 2005; Kunert-Bajraszewska et al. 2020; Orienti & Dallacasa 2020), the cores of FR I/FR II (Fanaroff & Riley 1974) radio galaxies (Chatterjee et al. 2011; Maccagni et al. 2020), and low-luminosity AGNs (Brunthaler et al. 2005; Mundell et al. 2009). Although the physics of the intrinsic variability of compact AGN jets remains an unsolved problem, plausible mechanisms include the propagation of shocks along the jet (Marscher & Gear 1985) and changes in the SMBH accretion properties such as accretion disk instabilities (Czerny et al. 2009; Janiuk & Czerny 2011) or accretion state changes (Koay et al. 2016; Wołowska et al. 2017). In addition to variability mechanisms directly related to accretion, radio variability may also arise from the jet itself owing to jet reorientation and the evolution of a young, expanding jet (e.g., An et al. 2020; Kunert-Bajraszewska et al. 2020). We focus on accretion-driven radio variability in this section and discuss jet reorientation and youth in Sections 5.5 and 5.6.

Each radio variability scenario described above is characterized by differences in variability amplitude level, as well as temporal and spectral evolution. Thus, distinguishing among them is best accomplished through multiepoch, broadband radio studies. For instance, low-amplitude (∼tens of percent) radio flares occurring on timescales of days to months are popularly attributed to shock propagation. The low variability amplitudes and short timescales typically associated with this flaring mode contrast with the properties of our sources, which were selected to exhibit large (100% to >2500%) flux increases at ∼GHz frequencies over decadal timescales and were later found to be steady over timescales of a few months.

Furthermore, the radio luminosities of our sources are typical of flares from AGNs with bright, persistent counterparts (such as the ∼2 Jy, 6 × 10³³ erg s⁻¹ Hz⁻¹ blazar J0851+202; Pietka et al. 2015). Blazar flares typically represent less than a twofold change in quiescent flux, which we emphasize is considerably lower than the typical variability amplitude observed in our sample (more than an order of magnitude in the most extreme cases).

For jets aligned at small angles to our line of sight, special relativistic effects cause the source to be beamed, which in the time domain leads to an apparent increase in the variability amplitude and a decrease in the variability timescale (Lister 2001). A rapid reorientation of a compact jet toward our line of sight during the ∼20 yr between FIRST and VLASS Epoch 1 would lead to an apparent brightening of an intrinsically low-luminosity radio source due to an increase in the Doppler factor.

Potential underlying causes for jet reorientation that may lead to variability on human timescales include helical magnetic fields, flaring blazars, jet−ISM interactions, and SMBH orbital motion. Confirming or refuting the possibility of helical magnetic fields requires multiepoch observations with milliarcsecond-scale resolution to identify key signatures such as periodic changes in the position angle of the jet (Britzen et al. 2017). We discuss the remaining possibilities in this section.

5.5.1. Blazar Contamination

5.5.2. Interaction with a Dense ISM

5.5.3. SMBH Binary Orbital Motion

Young radio AGNs, such as gigahertz-peaked spectrum (GPS) sources, are characterized by compact morphologies and inverted radio SEDs below their turnover (peak) frequencies, which are typically in the GHz regime (O’Dea 1998), consistent with the morphologies and radio SEDs of the majority of our sources. After a jet is launched, models predict a rapid increase in luminosity (P_radio ∼ t^2/5) as the dominant energy-loss mechanism transitions from adiabatic to synchrotron losses (An & Baan 2012), making the identification of young radio AGNs in VLASS that have emerged in the time since the FIRST survey (17–24 yr for our sample) plausible. For a nascent radio jet that has been triggered within the past 20 yr, the model of An & Baan (2012) suggests an increase in radio luminosity of >300%. The identification of such young radio AGNs is not unprecedented; the youngest known sources have kinematic ages as low as 20 yr (Gugliucci et al. 2005).

Based on the currently available radio continuum constraints, which include large variability amplitudes (with flux increases from 100% to >2500%) at 1.5 GHz compared to FIRST over timescales of decades, steady flux densities on timescales less than a few months at 3 GHz compared to VLASS, source size constraints <01 (≲1 kpc), and curved radio SEDs peaking at ∼5–10 GHz, we find the radio properties of our sources to be consistent with young and compact radio AGNs.

If the jet youth scenario for our sources is supported by higher-cadence, multiband follow-up data, there are exciting prospects for forthcoming radio surveys. Wide-area, broadband, synoptic radio surveys above a few GHz conducted over timescales of years to decades would be particularly well tuned for identifying large statistical samples of AGNs with recently launched jets.

In Figure 9, we plot our sources on the turnover–size relation, along with a sample of young radio AGNs from the literature. This relationship is believed to arise from the evolving radio spectral shapes of expanding young radio sources. As a young and compact radio source expands, the opacity due to SSA or FFA decreases, which causes the spectral peak (or turnover) to shift to lower frequencies (e.g., Bicknell et al. 1997; O’Dea 1998; Tingay & de Kool 2003). Since we only have upper limits on the linear sizes of our sources, we cannot yet directly confirm whether they follow the turnover–size relation. However, if we assume that our sources do follow the relation, we can obtain a rough estimate of their sizes.

Zoom In Zoom Out Reset image size

As an example, we consider J0807+2102 at z = 1.5588. The observed spectral turnover of this source is 21.9 GHz, which corresponds to a rest-frame turnover frequency of 56.2 GHz. Based on Figure 9, this implies a source size of 1–10 pc. Assuming a conservative value for the jet advance speed of 0.1c, we obtain a rough estimate of the source age of ∼30–300 yr. Thus, the young jet model passes this important consistency check.

Very long baseline interferometric (VLBI) studies with submilliarcsecond resolution will ultimately be needed to test whether our sources are indeed compact on parsec scales, as expected for young jets. Nevertheless, the radio SED shapes and source size limits from the VLA data presented in this study are consistent with radio variability arising from the evolving radio SEDs of jets that have recently been triggered.

The large observed changes in radio flux may be associated with AGNs transitioning between a radio-quiet state to one in which synchrotron-emitting radio jets are present. Extreme radio variability is typical of Galactic radio sources such as X-ray binaries (e.g., Mirabel & Rodríguez 1999). In these sources, the short timescales associated with black holes with masses ∼1–10 M_⊙ mean that the change in state from a radio-quiet mode associated with soft X-ray emission to a mode with a hard X-ray spectrum and radio jets can happen on timescales of minutes. Given the ∼10⁷–10⁹ M_⊙ SMBH masses that are typical of AGNs, similar transitions may occur (Maccarone et al. 2003; Falcke et al. 2004; Nipoti et al. 2005; Körding et al. 2006), but the corresponding timescales may be longer.

Previous attempts have been made to identify radio sources whose jets may have recently switched off (e.g., Marecki & Swoboda 2011), but such objects still have radio emission from their lobes (though such sources are interesting in their own right in the context of restarted jets). An alternative approach is to identify AGNs with a new onset of radio activity potentially associated with jets that have recently switched on. By their very nature these objects are expected to be relatively rare, but by surveying the radio emission from a large number of AGNs in two or more well-separated epochs, it may be possible to find objects that are candidates for AGNs undergoing this transition.

Although rare, previous studies have identified quasars with high-amplitude radio variability over timescales of years to decades (de Vries et al. 2004; Barvainis et al. 2005; Prandoni et al. 2010; Bell et al. 2015). In the Caltech NRAO Stripe 82 Survey (CNSS) pilot survey, Mooley et al. (2016) reported the detection of a single radio-loud type 2 quasar at z = 1.65 (VTC233002-002736) that lacked any detectable emission in the FIRST survey. Follow-up observations revealed a nearly 10-fold flux increase at 1.5 GHz over a 15 yr period and a curved radio SED, very similar in nature to our sources. Another example of such behavior in the CNSS, found in the z = 0.94 quasar CNSS 013815+00, has been reported recently by Kunert-Bajraszewska et al. (2020). A newborn expanding radio jet is responsible for GPS-like characteristics of the radio spectrum of 013815+00. In addition, the transition from the radio-quiet to the radio-loud phase in 013815+00 coincides with changes in its accretion disk luminosity. Thus, the burst of radio activity in 013815+00 is interpreted as a result of an enhancement in the SMBH accretion rate. We note that only 013815+00 is identified in the Pâris et al. (2018) SDSS quasar catalog, but both of the CNSS quasars are classified as AGNs in the mid-infrared by Assef et al. (2018) and would therefore meet the AGN and radio selection criteria of our sample (see Section 2).

We note that only one source matching our selection criteria was identified in the 50 deg² pilot CNSS footprint.³⁸ This is loosely consistent (within a factor of a few) with our observed detection rate so far over 3440 deg² of the VLASS-FIRST footprint of one confirmed optical or infrared AGN with newly radio-loud activity per ∼20 ± 13 deg². A more formal assessment of the areal density of AGNs that have transitioned from radio-quiet to radio-loud on decades-long timescales that incorporates the entirety of VLASS Epoch 1 will be presented in a future study.

The large-scale radio jets and lobes launched by some AGNs are believed to have long lifetimes and duty cycles on the order of millions of years. The properties of such “classical” radio AGNs contrast sharply with those of radio jets featuring compact (subgalactic) extents, younger ages (decades to thousands of years), and shorter lifetimes that have previously been identified (Lähteenmäki et al. 2018; Kunert-Bajraszewska et al. 2020).

In the case of reorienting jets, the change in jet direction may facilitate the transfer of energy over a large volume of the ISM of the host galaxy, thus potentially influencing ambient ISM conditions (Gaibler et al. 2011; Mukherjee et al. 2016). However, the prevalence and basic physical properties of reorienting jets, regardless of origin (see Section 5.5), have not been well constrained owing to inherent observational challenges (the combined requirements of high angular resolution imaging and high-cadence monitoring). If reorienting compact jets are common, particularly at z > 1, they may contribute substantially to SMBH−galaxy coevolution via the regulation of star formation from jet−ISM feedback or a disruption in SMBH growth in response to the launching of a jet.

Intermittent, albeit short-lived, episodes of radio-loud activity associated with powerful AGNs may also be important in the context of galaxy evolution. Numerous studies over the past few decades have used a variety of arguments (such as the overrepresentation of compact jetted AGNs in flux-limited surveys) in favor of the existence of a large population of short-lived and/or rapidly retriggered compact radio AGNs (Czerny et al. 2009; Mooley et al. 2016; Jarvis et al. 2019; Patil et al. 2020; Shabala et al. 2020). As a rough assessment of whether the rate of newly radio-loud AGNs identified in VLASS so far is consistent with this possibility, we compare the areal densities of the overall radio-loud AGN population based on constraints from FIRST with our sample. Ivezić et al. (2002) found 1154 SDSS quasars detected by FIRST over 774 deg², corresponding to an areal density of radio-loud quasars of ∼1 deg⁻². If these quasars have a lifetime of ∼10⁶ yr, we would only expect a very small fraction, ∼1 in 10⁵, to be within their first decade of being identified as radio-loud. This translates to an areal density of ∼10⁻⁵ deg⁻² compared to our detection rate of 13 newly radio-loud type 1 quasars in 3440 deg², i.e., ∼4 × 10⁻³ deg⁻², consistent with a typical period of occurrence of ∼10⁵ yr.

We speculate that frequent episodes of short-lived AGN jets that do not necessarily grow to large scales could be associated with higher-efficiency jet-driven feedback into the hosts of high-z galaxies. A similar conclusion was reached by van Velzen et al. (2015), who found that the fraction of powerful jets that grow beyond ∼100 kpc scales decreases with redshift. Future studies investigating the ISM content and conditions of the hosts of newly radio-loud AGNs will be important for placing quantitative constraints on the energetic impact of feedback from compact jets at cosmic noon.