High-resolution VLA Imaging of Obscured Quasars: Young Radio Jets Caught in a Dense ISM

We obtained spectroscopic redshifts for 71 out of 80 attempted sources using several telescopes (see Lonsdale et al. 2015, for details). The remaining 9 sources were too faint to provide a reliable redshift. Figure 1 shows the redshift distribution, which is seen to be approximately flat from 0.5 < z < 2 with a possible decline from 2 < z < 2.8. The median value is z_med ∼ 1.53. While the subset of sources targeted for redshift is likely biased to the optically brighter sources, it is unclear whether or not this translates to a bias in redshift—while optically brighter galaxies might be at lower redshift, optically brighter quasars might be at higher redshift. Taken at face value, our redshift distribution indicates that many of our sources lie close to the epoch of peak star formation and black hole fueling; some are nearer (z ≲ 1) and may be suitable for detailed follow-up observations.

Zoom In Zoom Out Reset image size

The 870 μm Atacama Large Millimeter/Submillimeter Array (ALMA) imaging of 49 sources (Lonsdale et al. 2015) and 850 μm James Clerk Maxwell Telescope (JCMT) Submillimetre Common-User Bolometer Array imaging of 30 sources (Jones et al. 2015) yielded 26/49 ALMA and 4/30 JCMT detections. Overall the MIR–submillimeter spectral energy distributions (SEDs) of our sample are likely to be dominated in the MIR by AGN-heated thermal dust emission. The selection of extremely red optical WISE colors and bright 22 μm emission revealed that these sources have high IR and bolometric luminosities (L_bol ∼ 10^11.7–10^14.2 L_⊙), with a few reaching the hyperluminous infrared galaxy regime. AGN populations identified using ultrared WISE color diagnostics are now known to belong to a class of IR-luminous obscured quasars such as Hot DOGs (e.g., Eisenhardt et al. 2012; Wu et al. 2012a; Assef et al. 2015). The MIR signatures and high-ionization lines in the spectra of our sample (Kim et al. 2013; E. Ferris et al. 2020, in preparation) are consistent with a population of radiative-mode obscured quasars. We refer our readers to Lonsdale et al. (2015) for more details on the MIR and submillimeter properties of our sample.

We observed 167 sources from Lonsdale et al. (2015) at X band (8–12 GHz) with the VLA in the A- and B-arrays through projects 12B-127 and 12A-064, respectively. Due to the complexity of dynamic scheduling for such a large sample, 12 sources were not observed in any array, and 32 were observed in only one array. Therefore, the sample discussed in this paper consists of 155 sources, 26 of which lack imaging with the A-array and 6 of which lack imaging with the B-array. The A-array observations were divided into 13 separate scheduling blocks (SBs), and a total of 129 sources were observed between 2012 October and December. The B-array observations were divided into seven different SBs, and 149 sources were observed from 2012 June to August.

Sources closer to each other on the sky were scheduled in groups, with phase calibrators interleaved. However, to maximize observing efficiency, the same calibrator was not always reobserved after each target. This strategy was worth the inherent risk of failing to obtain phase closure for a few targets, because most of the sources were expected to be bright enough for self-calibration.

The observations took place during the Open Shared Risk Observing period when maximum bandwidths were limited to ∼2 GHz. Our WIDAR correlator setup consisted of two basebands with central frequencies 8.6 and 11.4 GHz, respectively. The bandwidth of each baseband was 1024 MHz divided among eight 128 MHz wide spectral windows. The total bandwidth of our observations was 2 GHz. The correlator setup was kept identical for both of the arrays. Our observing strategy aimed to obtain snapshot imaging of the full sample with about 5 minutes of integration time per source with a theoretical rms noise level of about ∼13 μJy beam⁻¹.

We used the Common Astronomy Software Applications package (CASA; McMullin et al. 2007) version 4.7.0 for data editing, calibration, and imaging. The initial step was to remove bad data with the help of the VLA operators log,¹⁸ followed by visual inspection of the data in the uv-plane using the task PLOTMS. Hanning smoothing was performed prior to calibration to remove the rigging effect from the Gibbs phenomenon caused by strong radio frequency interference. The data were calibrated using the CASA VLA calibration pipeline¹⁹ (version 1.3.9).

We then used the pipeline weblog and test images of the targets and phase calibrators to examine the quality of the calibration. If necessary, additional flagging was done, followed by a rerun of the calibration pipeline. We then used the CASA task SPLIT to separate the uv-data for each target into individual data sets for self-calibration and final imaging.

We ran a few rounds of phase-only self-calibration and one round of amplitude and phase calibration to correct artifacts due to residual calibration errors. We used the CASA task CLEAN to produce the final continuum image. Because of the wide bandwidths made available by the new correlator, we formed images using the multifrequency synthesis mode with two Taylor coefficients (by setting the CLEAN parameter nterms = 2) to more accurately model the spectral dependence of the sky. Also, to mitigate the effects of non-coplanar baselines during imaging, we used the W-projection algorithm with 128 w-planes. The FWHMs of the synthesized beam of the final images in the A- and B-arrays are typically θ_b ∼ 02 and ∼06, respectively.

Despite our careful calibration and imaging strategy, a total of 13 targets (11 in A-array and 2 in B-array) suffered from severe phase closure issues. As a result, 110 sources have imaging in both arrays, 8 sources have only A-array imaging, and 37 sources have only B-array imaging. Thus, the analysis presented in the remainder of this paper is based on 155 sources.

To determine source parameters such as peak flux density, integrated flux, deconvolved shape parameters, and all corresponding uncertainties, we used the JMFIT task available in the 31DEC18 version of the Astronomical Image Processing Software (AIPS). In most cases, the radio sources have either single- or multicomponent Gaussian-like morphologies, and their flux and shape parameters may be estimated by fitting one or more two-dimensional elliptical Gaussian models. For sources with extended, complex structures, we manually estimated the source parameters using the CASA Viewer.²⁰ The flux measurement uncertainties were calculated by adding the error provided by JMFIT and the 3% VLA calibration error (Perley & Butler 2013) in quadrature. We provide the clean beam dimensions, peak flux, and total flux from our A- and B-array observations in Table C1.

The total flux distributions in A- and B-array observations span the ranges 0.18–45 mJy and 0.13–60 mJy, respectively, with similar medians of ∼3.3 mJy. Figure 2 compares the integrated fluxes of the 110 sources with high-quality flux measurements from both A- and B-arrays. The designation "high-quality" here simply indicates no hint of image artifacts.

Zoom In Zoom Out Reset image size

We find that, for most of our sample, the total flux measurements from each array are in good agreement. There are four sources that lie below the unity line in Figure 2 and have less flux recovered in the longer-baseline A-array observations. These sources may have a diffuse emission component that has not been recovered in the A-array data.²¹ There is also one outlier in Figure 2 with significantly higher flux in the B-array data compared to the A-array, possibly as a result of intrinsic source variability or calibration error.

We used the JMFIT task in AIPS to measure the angular sizes of our sources. For resolved sources, JMFIT²² requires that (1) the integrated flux be larger than the peak flux density and (2) the deconvolved major axis is greater than zero (within the relevant uncertainties). If neither of these criteria was satisfied, the source was classified as unresolved. The source-fitting algorithm gives a cautionary message when only one of two criteria is satisfied. We discuss our morphological classification in the next section, including our approach to sources with ambiguous JMFIT results.

Deconvolved source sizes were taken directly from JMFIT. The uncertainties were calculated based on the formalism given by Murphy et al. (2017):

$$\begin{eqnarray}&&\displaystyle \frac{{\sigma }_{\theta }}{{\sigma }_{\phi }}={\left[1-{\left(\displaystyle \frac{{\theta }_{b}}{\phi }\right)}^{2}\right]}^{-1/2},\end{eqnarray} \tag{ 1 }$$

where σ_θ and σ_ϕ are the rms errors on the deconvolved (θ) and measured (ϕ) source sizes, respectively. The parameter θ_b is the FWHM of the synthesized beam. For the unresolved sources, we consider the maximum deconvolved angular size provided by JMFIT to be an upper limit on the source size. For extended sources with non-Gaussian morphologies, we measured the angular sizes using the CASA Viewer. Table C2 provides the deconvolved source sizes and morphological classification. For sources with more than one component, separate measurements are given for each component.

As described in the previous section, the JMFIT task in AIPS uses two basic criteria to determine whether a source is formally resolved: the peak/total flux ratio and the deconvolved source size compared to the clean beam size. We use these criteria but modify the first to be more conservative by including a 3% uncertainty in the flux calibration (see Section 4.1).

We classify as "unresolved, U" sources that satisfy both criteria, deconvolved sizes consistent with zero in both axes and peak/total flux ratio of unity within the uncertainties. We classify as "slightly resolved" sources that show finite size along one of the two axes and a peak/total flux ratio consistent with 1. We classify as "resolved, R" sources that show finite size along both axes and a peak/total flux ratio less than 1 (within 1σ, following Owen 2018). Sources with more than a single distinct component are classified as double-, triple-, or multicomponent morphologies. Figure 3 shows the distribution of morphologies in our sample. We note that the entire analysis is performed separately for the A- and B-array data, and when possible, A-array results are preferred for the morphological classification and further analysis. In summary, we categorize our sample sources into the following morphological classes:

• 1.  
  Unresolved (UR): The source is unresolved along both the major and minor axes, and the peak/total flux ratio is unity within the 1σ uncertainty.
• 2.  
  Slightly resolved (SR): The source is unresolved along one of the axes, and the peak/total flux ratio is unity within the 1σ uncertainty.
• 3.  
  Fully resolved (R): The radio source is resolved along both the axes, and the peak/total flux ratio is <1.
• 4.  
  Double (D): The source consists of two distinct components, each of which may be unresolved, slightly resolved, or fully resolved.
• 5.  
  Triple (T): The source consists of three distinct components, resembling the core-jet or core-lobe emission seen in large-scale radio galaxies.
• 6.  
  Multiple (M): The source consists of more than three distinct components.

Zoom In Zoom Out Reset image size

Figure 3 shows the morphological classifications of the 155 sources in our final sample. Expressed as percentages, 55.5% ± 9.3% are unresolved, 13.5% ± 4.6% are slightly resolved, 7.7% ± 3.5% are fully resolved single sources, 14.8% ± 4.8% are double sources, 6.4% ± 3.1% are triple sources, and 1.3% ± 1.4% are multicomponent sources. Figure 4 shows the distribution of angular sizes for each morphological class. There is a wide range of upper limit sizes for the unresolved sources owing to the large span of source declinations and the use of both A- and B-array data. Deconvolved sizes are plotted for the slightly resolved sources, and outermost peak separation sizes are given for double, triple, and multiple sources.

Zoom In Zoom Out Reset image size

Our VLA X-band observations capture a wide range of frequencies, 8–12 GHz, offering the possibility of measuring "in-band" spectral indices, α (defined as f_ν ∼ ν^α). Although CASA generates a spectral index map with errors, we chose not to use it since its errors are calculated only as uncertainties to a polynomial fit and are less reliable at lower S/N (Cornwell et al. 2005; Rau & Cornwell 2011). Instead, we have chosen a more classical approach to estimate the in-band spectral index and its uncertainty. By dividing our bandwidth into two halves (centered at ν₁ = 8.6 and ν₂ = 11.4 GHz), we imaged each half separately using identical CLEAN parameters. We smoothed each 11.4 GHz image to match the resolution of the 8.6 GHz image using the task IMSMOOTH. We then regridded the smoothed 11.4 GHz image using the corresponding 8.6 GHz image as a template (using the CASA task IMREGRID) to ensure matched coordinate systems in the two images. Finally, we ran JMFIT to obtain source flux and shape measurements of all images.

The in-band spectral index was estimated using the following equation:

$$\begin{eqnarray}&&{\alpha }_{\mathrm{IB}}=\displaystyle \frac{{\mathrm{log}}_{10}({S}_{{\nu }_{1}}/{S}_{{\nu }_{2}})}{{\mathrm{log}}_{10}({\nu }_{1}/{\nu }_{2})},\end{eqnarray} \tag{ 2 }$$

where ν₁ and ν₂ are 11.4 and 8.6 GHz, respectively. Using standard propagation of errors, the uncertainty in the in-band spectral index is

$$\begin{eqnarray}&&{\sigma }_{{\alpha }_{\mathrm{IB}}}=\displaystyle \frac{{\left[{\left({\sigma }_{{S}_{1}}/{S}_{{\nu }_{1}}\right)}^{2}+{\left({\sigma }_{{S}_{2}}/{S}_{{\nu }_{2}}\right)}^{2}\right]}^{1/2}}{{\mathrm{log}}_{10}({\nu }_{1}/{\nu }_{2})}.\end{eqnarray} \tag{ 3 }$$

The left panel in Figure 5 shows the resulting uncertainty, ${\sigma }_{{\alpha }_{\mathrm{IB}}}$, plotted against the average S/N of the 8.6 and 11.4 GHz images. As expected, lower S/N yields larger uncertainties in α_IB with a threshold of S/N ≳ 70 for ${\sigma }_{{\alpha }_{\mathrm{IB}}}\lesssim 0.1$, which we take as a threshold of reliability for the calculated values of α_IB.

Zoom In Zoom Out Reset image size

Condon (2015) gives a theoretical analysis of in-band spectral indices and their uncertainties that broadly confirms our simple approach above. Combining Equations (48) and (49) from Condon (2015) for an in-band spectral index α_IB over a bandwidth of 8–12 GHz, we find

$$\begin{eqnarray}&&{\sigma }_{{\alpha }_{\mathrm{IB}}}\times {\rm{S}}/{\rm{N}}=\displaystyle \frac{\sqrt{12}}{\mathrm{ln}({\nu }_{\max }/{\nu }_{\min })}\,\sim 8,\end{eqnarray} \tag{ 4 }$$

where S/N is the signal-to-noise ratio of the source and ${\nu }_{\max }$ and ${\nu }_{\min }$ are the upper and lower ends of the observing bandwidth, respectively. The middle panel of Figure 5 shows the product ${\sigma }_{{\alpha }_{\mathrm{IB}}}\times {\rm{S}}/{\rm{N}}$ for our data and broadly confirms this result, with values near 7–8 for a range of in-band spectral indices.

The right panel of Figure 5 shows the distribution of in-band spectral indices with values above/below our S/N threshold color-coded as orange/blue. The distribution is strongly peaked near the median value of α_IB = −1.0, with 80% of the high-quality values within the range −1.7 to −0.5. We will discuss these spectral indices, together with the overall radio SEDs in a companion paper (P. Patil et al. 2020, in preparation). Briefly, the median spectral index is broadly consistent with optically thin synchrotron emission (α ∼ −0.7 near 1 GHz; e.g., Condon & Ransom 2016), perhaps steepened somewhat via radiative losses and inverse Compton scattering from either the cosmic microwave background or local infrared radiation fields. About 5% of our sources might plausibly have a flat spectrum (i.e., α > −0.5), consistent with an unresolved synchrotron core. This is also consistent with the absence of evidence for short-timescale variability typical of beamed sources, indicated by the good overall agreement between the fluxes measured in our A- and B-configuration observations. We will address the spectral characteristics and the role of beamed core emission more thoroughly in the SED paper.

Our sample was selected to have compact emission in the NVSS and FIRST catalogs. As discussed in Section 4, the majority of our sources have compact morphologies in our new high-resolution X-band observations. However, the presence of diffuse, extended emission on scales of a few arcseconds (which could be associated with earlier episodes of AGN activity) cannot be definitively ruled-out on the basis of the X-band data alone owing to surface brightness sensitivity limitations.

5.1.1. Constraints from Radio Surveys

5.1.2. NVSS and FIRST Flux Ratios

As a further test for missed emission in our X-band observations, we compare in Figure 6 the 1.4 GHz NVSS and FIRST fluxes of our sources. Excluding six sources that are resolved in FIRST but not in NVSS (J1025+61, J1138+20, J1428+11, J1651+34, J2145−06, J2328−2), the fluxes are in good agreement above 30 mJy with slight (∼5%) scatter to lower FIRST fluxes for weaker sources, with two outlier sources, J2322−00 and J1717+53, with flux ratios of 0.54 and 0.65, respectively. Neither of these sources shows any extended emission in TGSS, VLASS, or FIRST, and since both NVSS and FIRST were corrected for "CLEAN Bias," it cannot explain the offsets. We note that other sources of bias exist for measurements at low S/N (e.g., Hopkins et al. 2015). Variability might explain some of the outliers (Mooley et al. 2016), though we emphasize that Figures 2 and 6 indicate that the majority of our sources are not likely to be variable on the timescales sampled by our data.

Zoom In Zoom Out Reset image size

Figure 7 shows the distribution of physical source sizes, with the sample divided into resolved (including both slightly and fully resolved) and unresolved source morphologies. With the exception of 12 double or triple sources larger than 10 kpc, the rest are smaller than 5 kpc. Roughly 55% of the sources are unresolved with a median upper limit near 0.6 kpc. Given that our radio selection only requires sources to be compact on 40'' scales (NVSS, 100% of the sample) or 5'' scales (FIRST, 30% of the sample), we find that essentially all our sources are significantly more compact than these size limits, suggesting that our joint selection with luminous and red WISE MIR emission is preferentially associated with compact radio sources. A further check of whether the MIR selection is associated with compact radio emission is to ask whether an MIR blind radio survey with similar flux threshold and redshift range yields many compact sources.

Zoom In Zoom Out Reset image size

Such a survey exists. The CENSORS sample of Best et al. (2003) used NVSS to select sources brighter than 7.8 mJy and cross-matched these with the ESO Imaging Survey (EIS). The resulting sample of 150 has similar median redshift and radio luminosity to our sample. However, the median radio source size for the CENSORS sample is 6'', which is significantly larger than our own median source size of 01–02. Since the redshift distribution and flux cut for the two samples are similar, we conclude that the smaller source size of our sample is tied to the additional selection criteria of extreme MIR colors and luminosities.

Having established that our radio sources are compact, are there any previously established classes of radio sources that closely resemble our sources? Clearly they are different from the classical Fanaroff–Riley (FR) type I and II (Fanaroff & Riley 1974) radio sources, which are much more extended. Similarly, our sources, with their steep spectral index (Section 4.4), are also different from the compact flat-spectrum sources. There are four known classes of steep-spectrum radio sources that approximately match the angular and physical scales of our sample. These are the GPS (gigahertz-peaked spectrum; e.g., Fanti et al. 1990; O'Dea et al. 1991; Snellen et al. 1998; Fanti 2009; Collier et al. 2018), CSS (compact steep-spectrum; e.g., Peacock & Wall 1982; Spencer et al. 1989; Fanti et al. 1990, 2001; Sanghera et al. 1995), HFP (high-frequency peaker; e.g., Dallacasa et al. 2000; Stanghellini et al. 2009; Orienti & Dallacasa 2014), and FR0 classes (e.g., Sadler et al. 2014; Baldi et al. 2015). Of these, the FR0 class is significantly less luminous (<10²⁴ W Hz⁻¹; Baldi et al. 2018), and while the available GPS/CSS samples are somewhat more luminous than our sample (see next section), an SED analysis (P. Patil et al. 2020, in preparation) confirms that a significant fraction of our sources have curved or peaked spectra in the GHz range, similar to the GPS/CSS sources. Thus, since our sample seems to share a number of properties with the GPS/CSS sources, we will use these as a point of comparison in the following discussion.

Figure 8 presents the 1.4 GHz radio luminosity of our sample, which spans the range $25\lesssim \mathrm{log}({L}_{1.4\mathrm{GHz}}/{\rm{W}}\,{\mathrm{Hz}}^{-1})\,\lesssim 27.5$, with a median of $\mathrm{log}({L}_{1.4\mathrm{GHz}}/{\rm{W}}\,{\mathrm{Hz}}^{-1})\approx 26.3$. We also use Figure 8 to compare with other well-known samples of radio AGNs to help place our own sample within a wider "zoo" of radio sources.

Zoom In Zoom Out Reset image size

A representative sample of local (z < 0.3) radio AGNs was presented by Best & Heckman (2012), who cross-matched NVSS and FIRST sources with SDSS (radio luminosities calculated assuming a spectral index of −0.7). Clearly, our sample is roughly 2 dex more luminous than the local sample, confirming that our sample is much more luminous than the typical local radio AGN.

Next, we compare with the well-known low-frequency 3CRR survey, which is complete above S_{178 MHz} = 10.9 Jy (Laing et al. 1983). These span a wide range of redshift and luminosity and broadly divide into large-scale FR I and FR II radio sources (Fanaroff & Riley 1974). Our sample is, on average, 1.4 dex less luminous than the FR II galaxies and 1.3 dex more luminous than the FR I galaxies, though there is considerable overlap with both these samples.

Turning to radio sources that are, perhaps, better matched to the redshifts and physical scales of our own sources, the right side of Figure 8 includes samples of CSS and GPS sources (Spencer et al. 1989; Sanghera et al. 1995; O'Dea 1998; Fanti et al. 2001) and HFP sources (Dallacasa et al. 2000; Stanghellini et al. 2009). These samples show considerable overlap, though the median luminosities of the CSS, GPS, and HFP samples are larger by ∼1, 1.8, and 0.5 dex, respectively.

Overall, then, while our sample is significantly more radio luminous than typical radio AGNs, it has intermediate luminosity when compared to samples of powerful radio-loud AGNs.

An important property of a radio source that affects how it develops is its internal pressure. To first order, the measured pressure likely reflects the pressure of the surrounding medium into which the radio source is expanding. If the radio source is overpressured relative to the surrounding medium, perhaps being fed by a nuclear jet, then the radio source will expand.

To estimate the internal lobe pressures in our sample sources, we use relations derived from synchrotron theory given in Moffet (1975) and Miley (1980):

$$\begin{eqnarray}&&{P}_{l}\approx (7/9)({B}_{\min }^{2}/8\pi ),\end{eqnarray} \tag{ 5 }$$

where P_l is the pressure in the lobe of a radio source in dyne cm⁻² and B_min is the magnetic field in the magnetoionic plasma in G (gauss), derived using the common "minimum energy" or "equipartition" assumption that energy is shared approximately equally between the particles and the magnetic field. The equation for this magnetic field strength in G can be written as

$$\begin{eqnarray}&&{B}_{\min }\approx 2.93\times {10}^{-4}{\left[\displaystyle \frac{a}{{f}_{\mathrm{rl}}}\displaystyle \frac{{\left(1+z\right)}^{4-\alpha }}{{\theta }_{{rx}}{\theta }_{{ry}}}\displaystyle \frac{{S}_{\nu }}{{\nu }^{\alpha }}\displaystyle \frac{{X}_{0.5}(\alpha )}{{\theta }_{{ry}}{r}_{\mathrm{co}}}\right]}^{2/7},\end{eqnarray} \tag{ 6 }$$

where the radio source has flux S_ν in Jy with spectral form S_ν ∝ ν^α and angular size θ_rx × θ_ry arcsec, z is the redshift of the source, and r_co is the comoving distance in Mpc. We choose the filling factor for the relativistic plasma, f_rl, and the relative contribution of the ions to the energy, a, to be 1 and 2, respectively. The function X_0.5(α) handles integration over the frequency range from ν_l to ν_h, where ν_l = 0.01 GHz and ν_h = 100 GHz, and is defined as

$$\begin{eqnarray}&&{X}_{q}(\alpha )=({\nu }_{2}^{q+\alpha }-{\nu }_{1}^{q+\alpha })/(q+\alpha ),\end{eqnarray} \tag{ 7 }$$

where q is 0.5 in this case and represents the spectral shape function of the synchrotron emission.

Knowledge of the source size is required, since it feeds directly into the estimate of source pressure. For resolved single, double, or triple sources we take the measured region sizes directly from JMFIT. For slightly resolved or unresolved sources we take a conservative approach and use the beam major axis as an upper limit to source size. This yields a conservative lower limit for the source pressure. Higher-resolution Very Long Baseline Array (VLBA) images for a number of the unresolved sources (P. Patil et al. 2020, in preparation) usually reveal even smaller-scale double lobes with yet higher pressures. Thus, our current treatment of the VLA images yields useful, though conservative, lower limits to the radio source pressures in the unresolved sources.

Figure 9 shows the distribution of pressures for our sample, with lower limits for the unresolved sources. For the resolved sources, the median pressure is $\mathrm{log}({P}_{l}/(\mathrm{dyne}\,{\mathrm{cm}}^{-2}))$ = −7.2 or $\mathrm{log}[({P}_{l}/{k}_{B})/({\mathrm{cm}}^{-3}\,{\rm{K}})]$ = +8.7. For the lower limits, these values are $\mathrm{log}({P}_{l}/(\mathrm{dyne}\,{\mathrm{cm}}^{-2}))$ = −6.3 or $\mathrm{log}[({P}_{l}/{k}_{B})/({\mathrm{cm}}^{-3}\,{\rm{K}})]$ = +9.5.

Zoom In Zoom Out Reset image size

To help put our sample in context, we include the typical range of equipartition lobe pressures for a number of other classes of radio AGNs. On larger scales, the lobe pressures in FR I (e.g., Worrall & Birkinshaw 2000; Croston et al. 2008; Croston & Hardcastle 2014) and FR II (e.g., Croston et al. 2005; Harwood et al. 2016; Ineson et al. 2017; Vaddi et al. 2019) radio galaxies are roughly 3 dex lower than in our sample, almost certainly reflecting the much lower ambient pressures found on larger scales in the circumgalactic environment.

Figure 9 also shows the range of equipartition lobe pressures for CSS, GPS, and HFP sources taken directly from various studies (Mutel et al. 1985; Readhead et al. 1996; Orienti & Dallacasa 2014). There is a considerable overlap between our source pressures and those of the CSS, GPS, and HFP samples, possibly indicating a similarity in their properties and stage of development. However, a detailed comparison with these young radio AGNs is not straightforward because most measurements for the CSS, GPS, and HFP sources come from Very Long Baseline Interferometry observations with ∼milliarcsecond-scale angular resolution capable of identifying much more compact radio structures. Indeed, preliminary analysis of our own VLBA follow-up survey shows that many of our unresolved sources also have more compact source components with significantly higher pressures (∼1–3 dex; C. Lonsdale et al. 2020, in preparation). In all these comparisons, we have verified that our approach to measuring source pressures reproduces the source pressures given in these other papers.

The compact nature of the radio sources, together with their implied high pressures, seems to be a characteristic of the sample, and it is important to understand the origin of these high pressures. Unlike the lobes of extended FR I/FR II radio galaxies, the location of our radio sources deep within the host galaxy means that they are embedded within the relatively high-pressure environment of the central ∼1 kpc region. If the radio sources are in fact overpressured relative to the ambient ISM, then that overpressure may generate an expansion that, when coupled to the small size, may indicate a young source. We will present a more quantitative analysis of the source pressures and ages in Section 7.2 when we use a simple model of jet-driven lobe expansion to fit the observed source sizes and pressures.

Several models have been proposed to describe the temporal evolution of the observed properties of radio sources, such as luminosity and spectral turnover frequency (e.g., Falle 1991; Fanti et al. 1995; Readhead et al. 1996; Kaiser & Alexander 1997; O'Dea & Baum 1997; Snellen et al. 2000; Kaiser & Best 2007; Kunert-Bajraszewska et al. 2010; An & Baan 2012; Maciel & Alexander 2014). Many of the early models assumed self-similar expansion of radio jets as they move first through dense ISM during their initial growth until they emerge into the IGM and ICM to become large-scale, old sources (Kaiser & Alexander 1997). Early semianalytic models found that the radio source luminosity increases as ram-pressure-confined lobes expand within the galaxy. The luminosity reaches a maximum when the jets pass the boundary of the ISM, and then it decreases as the lobes expand into the ICM to become FR I/FR II sources.

We now explore the evolutionary stage of our sample and its connection to the FR I/FR II population by plotting our sources on the radio power versus linear size (PD) diagram in Figure 11. The range of linear extents of our sources covers multiple classes of medium- and compact-scaled radio sources, including CSS and GPS populations. It is clear from Figure 11 and Section 5.3 that the radio luminosities of our sample sources lie between those of the classical FR I and FR II populations. We also show the two tracks given by An & Baan (2012) that follow the high radio power (dotted) and low radio power (dashed–dotted) sources. Our sources, being intermediate in luminosity, lie between these two tracks in Figure 11. The dashed line shows the boundary between stable and unstable jets in the model of An & Baan (2012). The fact that all except one (J2318−25) of our sources lie above this line is consistent with them having stable jets that yield small-scale edge-brightened double or triple morphologies, as indeed we find in the majority of the resolved sources.

Zoom In Zoom Out Reset image size

Based on the evolutionary models given in An & Baan (2012) and following similar recent analyses (e.g., Jarvis et al. 2019), it seems that the position of our sources on the PD diagram relative to the jet instability criterion supports the possibility that they might eventually evolve into classical, FR I/II radio sources. This possibility is reinforced by the fact that our sources are heavily obscured, which points to a long-term fuel supply that could sustain the SMBH accretion for the ∼100 Myr time span necessary to create larger radio sources. However, a more careful discussion of possible evolutionary links between the WISE-NVSS sources and classical radio galaxies must consider the source ages. This we now attempt using a simple jet lobe expansion model.

There has been considerable work on models of radio source evolution in a variety of contexts, both analytic (e.g., Turner & Shabala 2015; Hardcastle 2018) and numerical (e.g., Mukherjee et al. 2016; Perucho 2019). Our sources may allow a relatively simple approach because the jets enter a dense, near-nuclear environment and are caught early in their development. While this may seem a potentially complex process, detailed simulations of just this situation (e.g., Mukherjee et al. 2016, 2018) suggest that the radio source develops in a quasi-spherical expansion, and in this case the analytic model of self-similar expansion is approximately correct.

A simple approach assumes purely adiabatic expansion, in which case the dynamics of the early phase of jet evolution can be approximated by the presence of a forward shock, a contact discontinuity, and an inner reverse shock. Following the mathematical treatment given in Weaver et al. (1977), a self-similar expansion of a spherical lobe can be expressed in terms of our observed parameters²⁴

$$\begin{eqnarray}&&{p}_{l}=7.76\times {10}^{-10}{F}_{43}{t}_{\mathrm{Myr}}{R}_{l}^{-3}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{p}_{l}=1.17\times {10}^{-9}\,{F}_{43}^{2/3}{n}_{a}^{1/3}{R}_{l}^{-4/3}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{p}_{l}=1.50\times {10}^{-12}{F}_{43}{\left({V}_{l}/c\right)}^{-1}{R}_{l}^{-2},\end{eqnarray} \tag{ 12 }$$

where p_l is the pressure inside the lobe expressed in dyne cm⁻², R_l is the radius of a lobe in kpc, F₄₃ is the mechanical jet power in units of 10⁴³ erg s⁻¹, n_a is the ambient number density in cm⁻³, ${t}_{\mathrm{Myr}}$ is a dynamical age in Myr, V_l is the lobe velocity, and c is the velocity of light.

While Sections 5.2 and 5.4 describe our estimates of radio source size, R_l, and pressure, p_l, estimating the jet power, F₄₃, is more uncertain. One approach is to assume that the jet power is related to the radio luminosity. While a number of studies have tried to establish such a link (e.g., Willott et al. 1999; Cavagnolo et al. 2010), others have argued that the relation is intrinsically quite scattered and has been amplified by selection bias (Godfrey & Shabala 2016). Bearing these caveats in mind, we cautiously adopt the relation given by Ineson et al. (2017):

$$\begin{eqnarray}&&{F}_{43}=5\times {10}^{3}{L}_{151}^{0.89\pm 0.09},\end{eqnarray} \tag{ 13 }$$

where L₁₅₁ is the rest-frame 151 MHz radio luminosity in units of 10²⁸ W Hz⁻¹ sr⁻¹. For our sample, we estimate L₁₅₁ using the 1.4 GHz luminosity from NVSS and a spectral index ${\alpha }_{1.4}^{10}$ derived from the NVSS flux and our X-band flux. We exclude sources with flat/inverted indices (${\alpha }_{1.4}^{10}\gt -0.3$) and low-S/N sources with very steep indices (${\alpha }_{1.4}^{10}\lt -2.0$) since the uncertainty in the extrapolation to rest-frame L₁₅₁ is large.

The left panel in Figure 12 shows the relation given in Equation (10) between source pressure, source size, and jet power, by plotting $3\mathrm{log}{R}_{l}-\mathrm{log}{F}_{43}$ against $\mathrm{log}{p}_{l}$ so that the source dynamical age, ${t}_{\mathrm{Myr}}$, appears as diagonal contours.

Zoom In Zoom Out Reset image size

Using estimates for R_l and p_l from Sections 5.2 and 5.4 and L₁₅₁ as described above, we find that the majority of our sources have dynamical ages in the range of 10⁴–10⁷ yr, with a median around 0.7 Myr. This is consistent with the overall picture that our sample contains young radio sources.

The middle panel of Figure 12 plots contours of external density, n_a, and the distribution of points reveals relatively high densities, spanning 1–10⁴ cm⁻³, comparable to the higher-density phases in spiral disks or near-nuclear ISM. Again, this is consistent with our overall picture of a young radio source emerging into a dense medium. Most radio sources are within a kiloparsec of a galaxy center, where we expect high average gas density, especially given the steep optical–MIR SED colors pointing to high columns. Indeed, we can combine the inferred ambient gas density with our measured source size to estimate a column density. The majority span $\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{2})\sim 22\mbox{--}25$, corresponding to A_V ∼ 5–5000, which is consistent with the red optical–MIR SEDs and the identification of Compton-thick columns in the related Hot DOG population (e.g., Stern et al. 2014; Ricci et al. 2017).

The right panel in Figure 12 plots contours of lobe expansion speeds. It seems that the sources expand with modest, subluminal speeds V_l ∼ 30–10,000 km s⁻¹ with a median near 450 km s⁻¹. We note that our velocities are also similar to those found in much more detailed simulations of a similarly powered jet interacting with a dense clumpy medium (Mukherjee et al. 2016, 2018).

The discussion of the growth of compact radio sources is often framed as two contrasting possibilities: the small sizes result from youth or from "frustrated" jets that cannot expand owing to a dense surrounding medium (e.g., van Breugel et al. 1984; Bicknell et al. 2018). Our analysis suggests that both perspectives might be relevant for our sources—the sources are indeed young, but the ISM is also dense, and this slows the source expansion.

Perhaps the most straightforward indication of youth would be to find a direct association with a short-lived phase in the host galaxy, such as a merger. Unfortunately, by selecting optically faint hosts (to avoid low-redshift sources), a simple inspection of the optical morphology is difficult. In the absence of direct observations, what might we expect? Despite early numerical simulations suggesting that luminous AGNs are associated with gas-rich mergers (e.g., Hopkins et al. 2008), the observational evidence has been mixed. For example, Cisternas et al. (2011) and Villforth et al. (2017) fail to find the AGN–merger connection. However, when the AGNs are selected to be dusty and obscured, such as WISE AGNs, the association with mergers is much clearer, particularly at high luminosity (e.g., Satyapal et al. 2014; Weston et al. 2017; Goulding et al. 2018).

Recent numerical simulations of galaxy mergers also support this association. Blecha et al. (2018) have tracked the evolution of WISE colors and luminosities for gas-rich mergers, finding the closest match to our sample's very red WISE colors during the brief final stage of coalescence.

Thus, our own sample of WISE-selected AGNs is very likely to contain a significant fraction of recent gas-rich mergers. Such a merger would be consistent with a newly triggered AGN with a radio jet.

Another approach that places the WISE-NVSS sources in a wider context is to use the RLF and dynamical age estimates to help establish a link to the other classes of radio source. First, the RLF analysis in Section 6 suggests that the WISE-NVSS sources have ∼300 times lower space density than classical radio galaxies. Second, comparing the median dynamical age of ∼10^5.8 yr to a typical age for a classical radio galaxy of ∼10⁷ yr suggests a source age ratio of ∼7%. Combining this age ratio with the ratio in space density of ∼0.3% indicates that ∼20% of the classical radio galaxies might have been born directly from a WISE-NVSS source.

Finding alternate compact progenitors that might evolve into classical radio galaxies is not hard. O'Dea (1998) performs a similar demographic analysis with GPS and CSS sources and shows that they actually overproduce the classical radio galaxies by a factor of ∼10. O'Dea (1998) interprets this apparent overproduction as evidence for recurrent activity in the GPS and CSS populations—meaning that there might be multiple phases of compact emission before the source finally evolves into a large-scale, classical, radio galaxy.

The relation between the WISE-NVSS sources and the GPS and CSS sources is not yet clear. There seems to be a systematic difference in the MIR properties (P. Patil et al. 2020, in preparation), suggesting that although all these sources may be dynamically young, the WISE-NVSS sources might be truly "newborn"—meaning that the radio source has emerged for the first time, into a dense near-nuclear ISM. In this case, the WISE-NVSS sources may either evolve directly into the classical radio galaxies or perhaps join the more common GPS and CSS classes, and from there ultimately evolve into a classical radio phase.