STORM IN A "TEACUP": A RADIO-QUIET QUASAR WITH ≈10 kpc RADIO-EMITTING BUBBLES AND EXTREME GAS KINEMATICS

We observed the "Teacup" AGN between 2013 December and 2014 May using the VLA at L- and C-band (≈1–2 GHz and ≈4–8 GHz respectively), in both the B- and A-configurations. These observations were taken under ID 13B-127, and were scheduled dynamically to ensure appropriate weather conditions. We bracketed our scans of the "Teacup" with two-minute scans of the bright, nearby BL Lacertae (BL Lac) object J1415+1320, which we use for phase referencing. We calibrated both the bandpass and the absolute flux scale using scans of J1331+3030 (3C 286) taken at the beginning and the end of each session. At L-Band, we recorded 16 contiguous spectral windows in full polarization of 64 × 1 MHz channels each, that yields a total instantaneous bandwidth of 1024 MHz. Accounting for losses due to flagging of radio frequency interference, the band center is 1.52 GHz. At C-band we recorded data in two frequency bands of eight contiguous spectral windows in full polarization, centered on ≈5 GHz and ≈7 GHz, with each spectral window consisting of 64 × 2 MHz channels, yielding 1024 × 2 MHz channels in total. The raw visibilities of the VLA data were calibrated locally using version 1.2.0 of the casa VLA Calibration pipeline.⁸

Below, we describe a number of maps that we created from our data at different resolutions and with different frequency centroids. The details of the synthesized beams and the noise in each map are summarized in Table 1. We also provide details of how we calculated flux densities, their uncertainties, and spectral indices for the radio structures that we identified. These values are summarized in Table 2.

Table 2. Radio Properties of Selected Structures

Structure	D	PA	S1.5 GHz	S5 GHz	S7 GHz	$\alpha ^{5}_{1.5}$	$\alpha ^{7}_{5}$	Variation in α
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
Matched resolution (Figure 1)
Core	1.5[2]	0	11.49[4]	3.462[16]	2.410[15]	−0.988[5]	−1.04[2]	±0.08
Eastern bubble	7.9[5]	66	9.85[4]	3.350[16]	2.219[15]	−0.888[5]	−1.18[2]	±0.3
Western bubble	6.9[5]	−90	2.45[4]	0.961[16]	0.685[15]	−0.77[2]	−0.98[8]	±0.6
High resolution (Figures 3 and 4)
HR-A	≲0.4	⋅⋅⋅	⋅⋅⋅	2.36[3]	1.59[3]	⋅⋅⋅	−1.13[6]	⋅⋅⋅
HR-B	≲0.4	⋅⋅⋅	⋅⋅⋅	0.59[3]	0.37[3]	⋅⋅⋅	−1.3[2]	⋅⋅⋅

Notes. The sizes, flux densities, and spectral indices of the structures identified in the VLA data of the Teacup AGN (see Figures 1 and 3). (1) The projected size of each structure (in arcseconds) based on the 5 GHz images, measured by taking a vector from the central core peak to the edge of each region (in the case of HR-A and HR-B we quote the HPBW); (2) the PA (in degrees) of the axis that the distances, D, were measured (for the eastern bubble this is in the direction to the brightest knot); (3)–(5) the observed-frame flux densities in mJy; (6) and (7) the spectral indices calculated from the quoted flux densities; and (8) the scatter of the α values within the defined regions (see Figure 1). In each column the number in square brackets gives the uncertainty on the last decimal place. We estimate an additional ≈5% (≈15%) uncertainty on the flux densities of the core and eastern bubble (western bubble) due to the uncertain definitions of the chosen regions (Figure 1).

Download table as:  ASCIITypeset image

In the top panel of Figure 1 we show the L-band, A-configuration VLA data (band center is 1.52 GHz), mapped with a robust = 0.5 weighting scheme and no additional tapering in the uv plane. We achieve a synthesized beam of ≈1 arcsec. In the second and third panels of Figure 1, we show maps from two subsets of the C-band data, centered at 5.12 GHz and 7.26 GHz, that are mapped by weighting the uv data with Gaussian tapers of fwhm = 120 kλ and 100 kλ, respectively, to achieve a close approximation with the resolution of the L-band imaging (see Table 1). In the bottom panel of Figure 1, we show a radio spectral index map (between 7 GHz and 5 GHz; $\alpha ^{7}_{5}$) from the full C-band spectral data set that was mapped using the multi-scale, multi-frequency synthesis (msmfs) imaging mode in casa. The resolution of this spectral index map is 2.05 × 1.78 arcsec (PA = 74°).

Zoom In Zoom Out Reset image size

From the three matched-resolution maps described above (Figure 1), we identify three spatial structures: a bright eastern bubble, a fainter western bubble, and a bright core. We measured the flux densities at 1.52 GHz, 5.12 GHz and 7.26 GHz of these features, from the maps shown in Figure 1, by integrating in polygons around each feature (with the edges of the polygons motivated by the 2σ contours in the 5.12 GHz map; see Figure 1). The measured flux densities are reported in Table 2, where the uncertainties correspond to the noise in the respective maps. Using our three-band radio photometry, we fit piecewise power-law SEDs between 1.52–5.12 GHz and between 5.12–7.26 GHz (Figure 2) and calculate the spectral indices ($\alpha ^{5}_{1.5}$ and $\alpha ^{7}_{5}$; Table 2). To reflect the variation of the spectral indices within each region, in Table 2, we also quote the 1σ scatter, within these regions, derived from the spectral index map shown in Figure 1. We note that the total 1.52 GHz flux density across all regions, seen in this ≈1 arcsec resolution L-band map, is S_1.5 GHz = 23.79 ± 0.07 mJy; however, by re-imaging our L-band data to a resolution 4.59 × 4.18 arcsec we recovered a total of S_1.5 GHz = 25.9 ± 0.1 mJy. We thus conclude that ≈2 mJy of the L-band radio flux in the "Teacup" AGN is produced in a diffuse, low surface brightness halo surrounding the core and bubbles, whose emission is resolved out in our higher-resolution continuum maps.

Zoom In Zoom Out Reset image size

To search for spatial structures on scales smaller than ≈1 arcsec, we make use of the C-band data at full resolution. In Figure 3 we show a robust = 0.5 weighted map, created from the B- and A-configuration, across the entire C-band, centered at 6.22 GHz. In this ≈0.5 arcsec resolution map the fainter western bubble is resolved out, while the core appears to be slightly elongated along a position angle of PA ≈ 60°, i.e., in line with the position angle connecting the core to the brightest region in the eastern bubble (see Table 2). Inset (a) of Figure 3 shows our highest, ≈0.3 arcsec, resolution C-band map of the core, comprising only of the A-configuration uv data that is mapped with super-uniform weighting in casa. In this map, the core is clearly resolved into two structures, which we label "High Resolution" (HR)-A and HR-B. HR-B, with peak position: 14^h30^m29902, +13°39'1215, is located 0.5 arcsec (≈0.8 kpc) to the northeast, at PA ≈ 60°, of HR-A which has peak position: 14^h30^m29873, +13°39'1190. To measure the 5 GHz and 7 GHz flux densities of these two structures we fit super-uniform weighted 5 GHz and 7 GHz maps with two beams (i.e., two unresolved components; Figure 4). The majority of the flux in the core is indeed inside two unresolved structures (corresponding to HR-A and HR-B) at both 5 GHz and 7 GHz. Therefore, these structures are unresolved at a half-power beam width of HPBW ≈ 0.4 arcsec (i.e., ≈0.6 kpc), which can be taken as very conservative upper limits on their intrinsic sizes. Figure 4 reveals low surface brightness emission (at ≈3σ) between the structures in the residual image. Furthermore, the combined 5 GHz flux density of HR-A and HR-B (S_5 GHz = 2.95 mJy; Table 2) is less than that of the total core structure from the lower resolution map (S_5 GHz = 3.46 mJy), implying a further ≈0.5 mJy of diffuse material the may exist inside the core that is not part of these two dominant structures.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We observed the "Teacup" AGN with the VIMOS (Le Fèvre et al. 2003) instrument installed on the ESO/VLT telescope between 2014 March 9 and 10 (Program ID: 092.B-0062). The seeing was <0.9 arcsec throughout. The observations were carried out in IFU mode using the 0.67 arcsec fiber, which provides a field-of-view of 27 × 27 arcsec, and the HR-Orange grism, which provides a wavelength range of 5250–7400 Å at a spectral resolution of ≈2650. By using the widths of sky-lines from our observations we found that observed resolution was FWHM = 105 ±8 km s⁻¹ around λ ≈ 5400 Å (i.e., around the observed wavelength of [O iii] λ5007). We corrected our measurements for this instrumental resolution. During the observations the target was dithered around the four quadrants of the VIMOS IFU. There were nine on-source exposures performed, each of 540 s, resulting in a total on-source exposure time of 4860 s. The on-source exposures were interspersed with three on-sky exposures, which were used for sky-subtraction. Flat field frames and arc-lamp frames were observed periodically during the observations. A standard star, taken under similar conditions to the science observations, was also observed.

The data reduction on the VIMOS data was performed using the standard esorex pipeline, which includes bias subtraction, flat-fielding, and wavelength calibration. The standard star was reduced in an identical manner to the science frames and esorex was used to apply the flux calibration. Data cubes were constructed from the individually sky-subtracted, reduced science frames and consequently the final data cube was created by median combining these cubes using a three-sigma clipping threshold.

Our IFU data for the "Teacup" AGN covers ≈40 × 40 kpc, which means that we can trace the kinematic structure over the full extent of the EELR. The analysis performed on the IFU data in this work closely follows the methods described in detail in Section 3 of Harrison et al. (2014); therefore, we only give brief details here.

To account for the complexity in the emission-line profiles, we opted to use non-parametric definitions to characterize the emission-line profiles within each pixel of the IFU data. In this work, we use the following definitions:

• 1.  
  The peak signal-to-noise ratio (S/N), (S/N)_peak, which is the S/N of the emission-line profile at the peak flux density. This allows us to identify the spatial distribution of the emission-line gas, including low surface-brightness features.
• 2.  
  The "peak velocity" (v_p), which is the velocity of the emission-line profile at the peak flux density. For the "Teacup" AGN, this traces the velocity structure of the dominant narrow component of the emission-line profile (see Harrison et al. 2014).
• 3.  
  The line width, W₈₀, which is the velocity that contains 80% of the emission-line flux. This is equivalent to 1.088×FWHM for a single Gaussian profile.
• 4.  
  The asymmetry value (A; see Liu et al. 2013), which is defined as

  $$\begin{equation} A \equiv \frac{(v_{90}-v_{50})-(v_{50}-v_{10})}{W_{80}}, \end{equation} \tag{ 1 }$$

  where v₁₀, v₅₀, v₉₀, are the velocities of at 10%, 50% and 90% of the cumulative flux, respectively. A very negative value of A means that the emission-line profile has a strong blue wing, while A = 0 indicates that it is symmetric.

To illustrate the kinematic structure of the ionized gas, using our IFU data, we use the [O iii] λ5007 emission line as it is the brightest line in the spectra and it is not blended with other emission lines. In Figure 5 we show four maps, produced using this [O iii] emission line, corresponding to the four non-parametric values described above. The values were calculated at each pixel of the data cube following the methods of Harrison et al. (2014). Briefly, we first fit the emission-line profiles (weighting against wavelengths with bad sky lines) with multiple Gaussian components before using the emission-line fits to define each of these four values. By using the emission-line fits, we are able to define our non-parametric values even in regions of low S/N (see Harrison et al. 2014). Values are assigned only to pixels where the emission-line was detected with (S/N)_peak ⩾ 5. In Figure 5, we re-sampled the pixel scale of the maps up by a factor of three, solely for the purposes of a clearer comparison to the morphology of the radio data.

Zoom In Zoom Out Reset image size

We extracted spectra from different spatial regions of the IFU data cube, where the regions were motivated by the morphology of the radio emission. These spectra were obtained by summing the spectra from the pixels in the vicinity of the three regions shown in Figure 6: (1) the eastern edge of the eastern bubble; (2) the western edge of the western bubble; and (3) the region around the HR-B radio structure. We fit the [O iii] emission-line profiles of these spectra in an identical manner to that outlined above for the individual spatial pixels. The fits to the spectra are shown in Figure 6, and the non-parametric values, derived from these fits are given in Table 3. We calculated uncertainties by fitting to 10⁴ mock spectra that were created by taking the best fit model and adding random Gaussian noise. We found that these formal uncertainties were very small (i.e., ≲10%) for the fits to emission-line profiles of the eastern bubble region and HR-B region and do not reflect the true uncertainties, which are dominated by the range in values within each region. Therefore, in Table 3 we provide the range of each value across the individual pixels within the regions.

Zoom In Zoom Out Reset image size

Table 3. [O iii] Emission-line Profile Properties from Selected Regions

Region	(S/N)peak	vp	W80	A
		(km s−1)	(km s−1)
HR-B region	640	147$_{-40}^{+11}$	920$_{-120}^{+50}$	−0.19$_{-0.15}^{+0.04}$
East bubble edge	805	150$_{-40}^{+70}$	290$_{-50}^{+20}$	0.14$_{-0.15}^{+0.02}$
West bubble edge	5	−110 ± 20	320 ± 70	0(⋆)

Notes. The peak S/N (S/N_peak), peak velocity (v_p), line width (W₈₀), and asymmetry (A) values for the [O iii] emission-line profiles extracted from labeled regions in Figure 6. The upper and lower limits on the values of HR-B and the eastern bubble reflect the full range in values of the individual pixels inside the defined regions. As the western bubble is not detected in individual pixels (Figure 5), we provide the 2σ uncertainties on the values for this region (see Section 3.2). ^⋆ As only one Gaussian component is fit to the emission-line profile of the western bubble, A = 0 by default (see Section 3.2).

Download table as:  ASCIITypeset image

Star formation is a possible candidate to produce both the radio emission in the core and also drive outflows via stellar winds and supernovae. The radio luminosity in the core of the "Teacup" AGN is L_1.4 GHz = 2 × 10²³ W Hz⁻¹, which is close to the proposed luminosity boundary between star-formation dominated and AGN-dominated radio emission in low redshift quasars (i.e., L_1.4 GHz ≈ 3 × 10²³ W Hz⁻¹ for 0.2 < z < 0.3; Kimball et al. 2011; Condon et al. 2013). The "Teacup" AGN has an FIR luminosity due to star formation of L_{IR, SF} = 3 × 10⁴⁴ erg s⁻¹ (Harrison et al. 2014). Assuming that the majority of this FIR emission is found in the central few kiloparsecs, the radio flux in the core region is a factor of ≈6 above the radio–FIR correlation of star-forming galaxies (e.g., Helou et al. 1985; Bell 2003; Ivison et al. 2010), which would indicate that star formation does not dominate the radio emission in the core.

Another possibility to explore is whether or not the HR-B radio structure (Figure 3) is associated with a compact star-forming region that is responsible for driving the high-velocity gas we see co-located with this structure (with v = −740 km s⁻¹ and v_max ≈ 1100 km s⁻¹; Figure 6). Indeed, extremely compact star-forming regions, with SFRs of ≈200 M_☉ yr⁻¹ and sizes ≈100 pc, have been attributed to driving ≈1000 km s⁻¹ outflows (Heckman et al. 2011; Diamond-Stanic et al. 2012; Sell et al. 2014). However, maximum velocities as high as that observed around the HR-B radio structure (i.e., v_max ≈ 1100 km s⁻¹) are not seen in non-AGN host galaxies with SFRs comparable to that of the "Teacup" AGN (i.e., SFRs ≈ 10 M_☉ yr⁻¹; e.g., Genzel et al. 2011; Arribas et al. 2014) and compact starbursts are often associated with flatter radio spectral indices than that observed here (Murphy et al. 2013). We further disfavor the star formation scenario because: (1) there are not two continuum peaks observed in the HST imaging (i.e., to correspond to both the HR-A and HR-B radio structures); (2) HR-B is unlikely to be an obscured star-forming region because the dust in this source is not located at this position (see Figure 6; also see Keel et al. 2014); and (3) the high-velocity component has an [O iii]/Hβ emission-line flux ratio that is consistent with being dominated by AGN photoionization rather than a H ii region (see Harrison et al. 2014).

A possible origin of the radio emission in radio-quiet quasars is accretion disk winds (see, e.g., Shlosman et al. 1985; Murray et al. 1995; Proga et al. 2000; Crenshaw et al. 2003), which propagate through the surrounding interstellar medium (ISM) to produce non-thermal radio emission (Zakamska & Greene 2014; Nims et al. 2014). Based on the energy conserving outflow model of Faucher-Giguère & Quataert (2012), Nims et al. (2014) explored the observational signatures of these winds in the radio and X-rays. In their model, the wind causes a shock that propagates outward into the ambient medium and can reach distances of ≳1 kpc with velocities ∼1000 km s⁻¹ (at this point it interacts with the ISM and decelerates). They predict a steep-spectral index of α ≈ −1, which is consistent with our observations (see Table 2). It is less clear what this model predicts for the morphology of the radio emission; however, the propagation of the wind is predicted to be influenced by the properties and distribution of the ISM, and it will preferentially escape along the path of "least resistance" (Faucher-Giguère & Quataert 2012), which could lead to a moderately collimated outflow. For the "Teacup" AGN, we find that the position of the ionized outflow, and the associated radio emission, is directed away from the dusty regions (Figure 6), which may be consistent with these predictions.

In the fiducial model of Nims et al. (2014), the interaction of the winds with the ISM is predicted to result in a radio luminosity of

$$\begin{equation} \nu L_{\nu } \approx 10^{-5} \xi _{-2} L_{{\rm AGN}} \left(\frac{L_{{\rm kin}}}{0.05L_{{\rm AGN}}}\right), \end{equation} \tag{ 2 }$$

where ξ₋₂ is the scaled efficiency (i.e., 1% of the shock energy goes into relativistic electrons), L_AGN is the bolometric luminosity of the AGN and L_kin is the kinetic luminosity of the wind. For the Teacup AGN, Nims et al. (2014) would predict a radio luminosity of νL_ν ≈ 20 × 10³⁹ erg s⁻¹ (i.e., L_1.5 GHz ≈ 13 × 10²³ W Hz⁻¹), using their fiducial values (i.e., ξ₋₂ = 1; L_kin = 0.05L_AGN). This is a factor of ≈8 above the core radio luminosity that we observe of νL_ν(1.5 GHz) = 2.6 × 10³⁹ erg s⁻¹ (i.e., L_1.5 GHz = 1.7 × 10²³ W Hz⁻¹). Given the nature of the simple model assumptions (e.g., the assumed properties of the ISM; see Nims et al. 2014), we cannot use this discrepancy as strong evidence that the radio emission in the core of the "Teacup" AGN, and the associated outflow, are not driven by a quasar wind.

Compact radio emission (i.e., on ≲5 kpc scales) in radio-quiet quasars may be associated with nuclear coronal activity (i.e., the region around the AGN accretion disk) or compact radio jets (e.g., Kukula et al. 1998; Ulvestad et al. 2005; Leipski & Bennert 2006; Laor & Behar 2008). Nuclear coronal activity is predicted to produce flat spectral indices (i.e., α ≈ 0; Laor & Behar 2008), which is not what is observed in the "Teacup" AGN (Figure 2).

In contrast, the core emission in the "Teacup" AGN shares many characteristics with radio-loud (with L_5 GHz ≳ 10²⁶ W Hz⁻¹) compact, steep spectrum quasars, that show steep radio spectra, have jets and lobes on small scales and exhibit broad [O iii] emission-line profiles (e.g., Gelderman & Whittle 1994; see review in O'Dea 1998). In particular, it is plausible that the HR-B radio structure is a hot spot of a radio jet, where it interacts with the ionized gas in the local ISM and drives a high-velocity outflow (see Figure 6). This is also analogous to previous observations of a small number of AGNs showing compact jets interacting with the ionized gas (e.g., Capetti et al. 1999; Leipski & Bennert 2006; Barbosa et al. 2009; Stockton et al. 2007). We note that several radio-loud systems have also been observed to have high-velocity kinematic components of ionized gas, associated with radio jets, that are offset by several hundred km s⁻¹; furthermore, in the most extreme cases the line widths for these radio-loud systems can be even greater than observed in the "Teacup" AGN (i.e., FWHM ≳ 1000 km s⁻¹; e.g., Emonts et al. 2005; Holt et al. 2008; Nesvadba et al. 2008; Shih et al. 2013).

Overall the radio emission and outflow in the core of the "Teacup" AGN are consistent with being produced by radio jets or quasar winds that are interacting with the ISM in the host galaxy. This interaction could result in the high-velocity ionized gas that we observe (see Figures 6 and 7). The compact sizes and the narrow angular range of the radio emitting region in the core, in particular for the HR-B structure, may favor the jet scenario.