Double-peaked Emission Lines Due to a Radio Outflow in KISSR 1219

The VLA observations of KISSR 1219 (Director's Discretionary Project 15A−468) were carried out on June 4, 2015 with the transition $\mathrm{BnA}\to {\rm{A}}$ array-configuration at the L band (1.326 GHz). The full bandwidth was split into 64 channels with a bandwidth of 2 MHz each, and 10 SPWs (spectral windows) ranging from 0.942 to 1.966 GHz. The data integration time was 2 s and the total on-source time was ∼15 minutes. J1219+4829 and J1331+3030 were used as the phase and amplitude calibrators, respectively. We calibrated the data using the Common Astronomy Software Applications (CASA) package version 4.6.0. We followed the CASA tutorial on VLA wide-band wide-field imaging.¹⁰ Data from shadowed antennas and zero-amplitude data were flagged using FLAGDATA. Before flagging the data affected by RFI (radio frequency interference), we Hanning-smoothed the data using CASA task HANNINGSMOOTH. We subsequently carried out automatic RFI excision using RFLAG in FLAGDATA, where we examined all the SPWs and scans individually. We eventually removed all data from SPW 2 (centered at ν = 1.646 GHz) and most of the data from SPW 6 (centered at ν = 1.006 GHz), since it was heavily affected by RFI. Standard calibration tasks in CASA were used to calibrate the data.

We made several images using multi-scale clean, multi-scale-wide-field clean with W-projection, multi-scale-multi-frequency-synthesis, multi-scale-multi-frequency-wide-field clean, using imagermode = ''csclean'' (Cotton-Schwab CLEAN) as well as ''mosaic.'' The image made using multi-scale, multi-frequency-synthesis and imagermode = ''mosaic'' had the lowest rms noise and highest source flux density. This is presented in the top panel of Figure 1. Self-calibration did not work for this weak source.

Zoom In Zoom Out Reset image size

The VLBI observations were carried out in a phase-referencing mode with seven antennas¹¹ of the Very Long Baseline Array (VLBA) at 1.5 GHz on 2015 February 06 (Project ID: BK191A), and eight antennas¹² at 4.9 GHz on 2015 February 08 (Project ID: BK191B). The data were acquired with an aggregate bit rate of 2048 Mbits s⁻¹ (2 polarizations, 16 baseband channels, bandwidth 32 MHz, and a 2-bit sampling rate). 1218+444, which is 222 away from the source and has an x, y positional uncertainty of 0.16, 0.32 mas, was used as the phase reference calibrator. 4C 39.25 and 1206+416 were used as the fringe-finder and phase-check calibrator, respectively. A nodding cycle of 5 minutes (2 minutes on the phase calibrator and 3 minutes on the target) was used for the observations, interspersed by four 5-minute scans on 4C 39.25 and one 3-minute scan on the phase-check calibrator. The experiment lasted a total of ≈4.0 hr at each frequency.

The data reduction was carried out using the Astronomical Image Processing System (AIPS) following standard VLBA-specific tasks and procedures outlined in the AIPS cookbook.¹³ Los Alamos (LA), which was a stable antenna in the middle of the configuration, was used as the reference antenna for the calibration. The amplitude calibration was carried out using the procedure VLBACALA, while the delay, rate, and phase calibration were carried out using the procedures VLBAPANG, VLBAMPCL, and VLBAFRGP. The phase calibrator 1218+444 was iteratively imaged and self-calibrated on phase and phase+amplitude using AIPS tasks IMAGR and CALIB. The images were then used as models to determine the amplitude and phase gains for the antennas. These gains were applied to the target and the final images were made using IMAGR. A round of data-flagging was carried out using the task CLIP on the source SPLIT file, prior to making the images. The radio component detected at 1.5 GHz was offset from the center of the image by 00, −020 in R.A. and decl. We ran the task UVFIX on the source SPLIT file to shift the source to the center before producing the final image. We used the elliptical Gaussian-fitting task JMFIT to obtain the component flux densities.

As the source was weak at 1.5 GHz, we specified three boxes (using parameter CLBOX) in the task IMAGR to run CLEAN. The choice of these three components was made to produce an image with low and uniform residual noise. Self-calibration could not be carried out for this weak source. The final naturally weighted image made using ROBUST = 0 and a circular beam of 17 mas, is presented in the bottom panel of Figure 1.

The observed SDSS optical spectrum was corrected for reddening using the $E(B-V)$ value from Schlegel et al. (1998). The pPXF (Penalized Pixel-Fitting stellar kinematics extraction) code by Cappellari & Emsellem (2004) was used to obtain the best-fit model of the underlying stellar population, in order to isolate the AGN emission lines. Emission lines in the dereddened spectrum were masked and the underlying absorption spectrum was modeled as a combination of single stellar population templates with MILES (Vazdekis et al. 2010); these templates are available for a large range of metallicity (M/H ∼ −2.32 to +0.22) and age (63 Myr to 17 Gyr).

pPXF provides the best-fit template and the stellar velocity dispersion of the underlying stellar population (${\sigma }_{\star }$), which was 142 ± 6.9 km s⁻¹ for KISSR 1219. Although the major contribution to the observed spectrum is the underlying stellar population (∼80%), there could be other contributors like a power-law component from AGN and Fe lines (∼20%). While estimating the best fit to the underlying population, pPXF fits a polynomial along with the optimal template to account for the contribution from the power-law continuum. The Fe lines are too weak to affect further analysis. The output of the analysis is shown in the lower panel of Figure 2. The reddening-corrected spectra is shown in black and the best-fit model is over-plotted in blue. The best-fit model is subtracted from the dereddened spectrum to obtain the pure emission line spectrum, shown in the lower panel in green. The prominent emission lines in the pure emission line spectrum such as, [S ii] λλ6717, 6731 doublet, [N ii] λλ6548, 6584 doublet, Hα, [O i] λλ6300, 6364 doublet, [O iii] λλ5007, 4959 doublet, Hβ and Hγ are shown in the middle and upper panels of Figure 2. Signatures of double peaks are seen in all of these emission lines.

Zoom In Zoom Out Reset image size

In order to measure the line parameters, the pure emission line spectrum is analyzed by modeling the profiles as Gaussians. In general, the [S ii] doublet lines, which are well separated, are considered as good representations of the shape of the [N ii] and Hα narrow lines (Filippenko & Sargent 1988; Greene & Ho 2004). Hence, we first modeled the [S ii] lines to obtain a satisfactory fit to the line-shape and used that model as a template for other narrow lines. The two [S ii] lines are assumed to have equal width (in velocity space) and are separated by their laboratory wavelengths. We also tried to fit each [S ii] line with double and triple component Gaussian models. When multiple Gaussian components are included, the corresponding components of each line are assumed to have equal width (in velocity space) and are separated by their laboratory wavelengths. Also, the relative intensities of different components of each line are held fixed. The best-fit to each [S ii] line required three Gaussian components, as shown in the lower right panel of Figure 3. The two narrow components in each [S ii] line represent the two lines corresponding to the double peak feature. The third component represents the extended wing, to which we do not attribute any physical significance. This component was essential to fully describe the [S ii] line-shape and its inclusion is statistically justified by an improvement in the reduced χ² value by ∼20%.

Zoom In Zoom Out Reset image size

The [S ii] model was finally used as a template to fit the narrow Hα and [N ii] doublet lines. The widths of the different components of the [N ii] doublet and Hα narrow lines were assumed to be the same as that of the corresponding components of the [S ii] lines. The separation between the centroids of the [N ii] narrow components were held fixed and the flux of [N ii] λ6583 to [N ii] λ6548 was fixed at the theoretical value of 2.96. Because there are three components for each [S ii] line, the profiles of Hα and [N ii] doublet narrow lines were strictly scaled from [S ii]. The best-fit profile and individual components to the Hα and [N ii] narrow lines are shown in the lower-middle panel of Figure 3.

The Hβ and Hγ line profiles were also modeled using the [S ii] template, but now the width of the third Gaussian component corresponding to the extended wing, was kept as a free parameter. As the fit improved by 18% in reduced χ² compared to the model where we kept the width fixed, and as we do not attribute any physical significance to the third Gaussian component, the choice to keep the width as a free parameter was justified. The best-fit profile and individual components to the Hβ and Hγ lines are shown in upper-middle and upper-left panels of Figure 3, respectively.

Because the [O iii] profile does not typically match that of [S ii] (Greene & Ho 2005), we fit the [O iii] and [O i] independently. We initially used a two Gaussian component model and a third additional Gaussian component was included if it improved the corresponding fit by >20%. As shown in the lower-left and upper-right panels of Figure 3, each line of the [O iii] doublet is modeled using three Gaussian components and the [O i] doublet lines using two Gaussian components. In the case of the [O iii] doublet lines, the third component was significantly blueshifted and broad. Hence, these lines may be representing an outflow feature.

IDL programs, which use the MPFIT function for non-linear least-square optimization were used to fit the emission line profiles with Gaussian functions and to obtain the best-fit parameters. The fit parameters corresponding to each line are presented in Table 1. The errors associated with each parameter are those provided by MPFIT.

The 1.5 GHz VLBA radio core position of KISSR 1219 is at R.A. 12^h09^m0880956, decl. 44° 00' 113057. The peak intensity of the 1.5 GHz radio core is ∼0.26 mJy beam⁻¹; the total flux density of the entire core-jet structure is ∼1.58 mJy. While KISSR 1219 was not detected at 5 GHz with the VLBA, the phase calibrator 1218+444 was detected. The final rms noise at 1.5 and 5 GHz is ∼40 μJy beam⁻¹ and ∼30 μJy beam⁻¹, respectively.

Using the total flux densities of the core and two jet components J1 and J2, the brightness temperature (T_B) at 1.5 GHz (see the relation in Ulvestad et al. 2005) turn out to be ∼6 × 10⁶ K, ∼3 × 10⁶ K, and ∼4 × 10⁶ K, respectively, for the unresolved components. These relatively high brightness temperatures support non-thermal AGN-related emission. While VLBI has indeed detected individual radio supernovae and supernova remnants in starburst galaxies with similarly high brightness temperatures (≥10⁶ K, Lonsdale et al. 2006; Pérez-Torres et al. 2009), the close alignment of the parsec- and kiloparsec-scale radio jets in KISSR 1219, support the AGN picture instead.

Because no core-jet emission was detected in the 5 GHz image, we infer the 1.5–5 GHz spectral indices to be steeper than −1.2, −1.1, and −1.0 for the core, J1 and J2, respectively, assuming three times the rms noise in the 5 GHz image as an upper limit to the flux density of these components. The errors in the spectral index values are ≈20%. It is known that VLBI "cores" in some Seyferts do not exhibit flat or inverted spectra and can have steep spectra instead (e.g., Roy et al. 2000; Kharb et al. 2010; Orienti & Prieto 2010; Bontempi et al. 2012; Panessa & Giroletti 2013). This is unlike the cores in powerful radio sources, which are suggested to be the unresolved bases of relativistic jets.

We remade the VLBA images of KISSR 1219 using different weighting schemes by varying the ROBUST parameter, in order to check for the compactness of the core. We found that the core, or rather a compact portion of it, remained intact and centrally concentrated in all the images, even as the jet components (J1 and J2) got largely resolved out in ROBUST −5 (pure uniform weighting) images. This result supports the suggestion that the brightest feature in the bottom panel of Figure 1, is indeed the radio core, albeit with partially resolved jet emission in close proximity, which makes the overall core spectrum steep.

The peak intensity of the 1.5 GHz VLA radio core is ∼3.35 mJy beam⁻¹; the total flux density of the entire core-jet structure is ∼2.33 mJy. The final rms noise in the 1.5 GHz VLA image is ∼100 μJy beam⁻¹. Assuming equipartition of energy between relativistic particles and the magnetic field (Burbidge 1959), we obtained the minimum pressure, and the particle energy (electrons and protons) at minimum pressure using the relations in O'Dea & Owen (1987). The total radio luminosity was estimated assuming that the radio spectrum extends from (ν_l=) 10 MHz to (ν_u=) 100 GHz with a spectral index of α = −1. It was assumed that the ratio of the relativistic proton to relativistic electron energy was unity. Table 2 lists equipartition estimates for plasma filling factors of unity and 0.5 (e.g., Blustin & Fabian 2009). The total energy in particles and fields, E_tot, is estimated as E_tot = 1.25 × E_min, while the total energy density, U_tot, in erg cm⁻³ is =E_tot(ϕV)⁻¹. The minimum energy magnetic fields are of the order of ∼1 mG on the VLBA scales and ∼30 μG on the VLA scales, while the total energy E_tot is ≈4 × 10⁵³ erg and ≈7 × 10⁵⁵ erg, respectively. The lifetime of electrons in the radio component undergoing both synchrotron radiative and inverse Compton losses on CMB photons was estimated using the relation in van der Laan & Perola (1969). These are of the order of 10,000 years on the VLBA scales and ∼4 Myr on the VLA scales.

If we assume that the kiloparsec-scale emission is emission from a broad AGN outflow in KISSR 1219, then we can estimate the time-averaged kinetic power (Q) of this outflow following the relations for radio-powerful sources listed by Punsly & Zhang (2011; see also Willott et al. 1999): we derive F₁₅₁ using the 1.4 GHz flux density from FIRST and a jet/lobe spectral index of −0.8, and obtain a kinetic power of Q ≈ 2.5 × 10⁴¹ erg s⁻¹. This Q value is typical of outflows in low-luminosity AGN (e.g., Mezcua & Prieto 2014).

We did not see any connection between the VLA jet and the diffuse emission ≈63 arcsec (46 kpc) to the east of KISSR 1219, and at a P.A. similar to the radio jet (Figure 4). Instead, we find a weak DSS optical source but a strong WISE infrared source (in bands 1 and 2 at 3.4 μm and 4.6 μm, respectively) in the center of the diffuse emission. We therefore conclude that this diffuse emission is unrelated to KISSR 1219. It comes from a possibly higher redshift radio galaxy.

We can rule out the radio structure on kiloparsec scales to be emission from the spiral arms of the galaxy, on the basis of its one-sidedness. Because the galaxy is nearly face-on, it is difficult to explain why the counter spiral arm is not emitting radio waves. The straight VLA radio feature, which is exactly coincident with the VLBA core-jet structure is consistent with it being a synchrotron jet from the AGN in KISSR 1219.

In order to test if the one-sided VLA jet in KISSR 1219 could be a result of relativistic beaming effects, we estimated the jet-to-counter-jet surface brightness ratio (R_J) at a distance of 6'' (=4.37 kpc) from either side of the core (∼0.46 mJy beam⁻¹ on the jet side versus ∼0.19 mJy beam⁻¹ on the counter-jet side). For the jet structural parameter of p = 3.0 (continuous jet with α = −1; Urry & Padovani 1995), R_J is 2.5. We can estimate a similar R_J on VLBI scales; R_J = 8.5 for the jet component (J1) at a distance of ∼40 mas (=30 pc) from the radio core. Assuming that the torus half-opening angle is ∼50° (e.g., Simpson et al. 1996), and therefore the orientation of KISSR 1219 should be greater than ∼50° for it to be classified as a Seyfert 2, we can derive lower limits on the jet speeds on both parsec- and kiloparsec scales, in order to produce the observed R_J values. The lower limits of jet speeds are ∼0.55c and ∼0.25c on the VLBA and VLA scales, respectively. Because these are viable values of jet speeds in Seyfert galaxies (Ulvestad et al. 1999), we cannot rule out Doppler boosting (and dimming) effects being responsible for the one-sided jet structure in KISSR 1219.

If the missing counter-jet emission on the other hand, is a result of free–free absorption, the required electron densities of the ionized gas on both parsec- and kiloparsec scales can be estimated, using the relations, $\mathrm{EM}=3.05\times {10}^{6}\,\tau \,{T}^{1.35}\,{\nu }^{2.1}$, and ${n}_{e}=\sqrt{\mathrm{EM}/l}$ (Mezger & Henderson 1967). Here EM is the emission measure in pc cm⁻⁶, τ is the optical depth at frequency ν in GHz, T the gas temperature in units of 10⁴ K, n_e the electron density in cm⁻³, and l is the path length in parsecs. In order to account for the observed jet-to-counter-jet surface brightness ratios of 8.5 on parsec-scales and 2.5 on kiloparsec scales, the optical depth at 1.5 GHz should be at least ∼2 and ∼1, respectively.¹⁴ For a gas temperature of 10⁴ K and a path length of 1 pc and 100 pc for the VLBA and VLA jets, respectively, EM of ≈1.4 × 10⁷ pc cm⁻⁶ and 7.1 × 10⁶ pc cm⁻⁶, and n_e of ≈4000 and 250 cm⁻³ are required for free–free absorption, on the VLBA and VLA scales, respectively. Such ionized densities can indeed come from gas clouds in the narrow-line region (see also Section 3.6 ahead). However, the volume filling factor of NLR clouds is small—of the order of 10⁻⁴ (e.g., Alexander et al. 1999)—making them unlikely candidates for absorbers, especially on kiloparsec scales where the jet width becomes larger (>1 kpc). Ionized gas in giant H ii regions with n_e ∼ 100–1000 cm⁻³ and lifetimes ∼10⁷ years could, in principle, also be the candidate media for free–free absorption (Clemens et al. 2010). The volume filling factor of this gas is also small, ≥0.2 (e.g., Walterbos & Braun 1994). We conclude that Doppler boosting/dimming effects are a more favorable explanation for the one-sided jet structure observed in KISSR 1219.

The bulge stellar velocity dispersion as derived by our line fitting is σ_⋆ = 142.0 ± 6.9 km s⁻¹. Following the ${M}_{\mathrm{BH}}\mbox{--}{\sigma }_{\star }$ relation for late-type galaxies by McConnell & Ma (2013), $\mathrm{log}\left(\tfrac{{M}_{\mathrm{BH}}}{{M}_{\odot }}\right)=8.07+5.06\mathrm{log}\left(\tfrac{{\sigma }_{\star }}{200\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)$, the mass of the BH in KISSR 1219 turns out to be (2.1 ± 1.5) × 10⁷ M_⊙. The Eddington luminosity (L_Edd ≡ 1.25 × 10³⁸ M_BH/M_⊙) is ≈2.6 × 10⁴⁵ erg s⁻¹. The bolometric luminosity estimated using the [O iii] λ5007 line flux is ∼6.19 × 10⁴³ erg s⁻¹. The Eddington accretion rate (≡L_bol/L_Edd) is ∼0.02, which is typical of low-luminosity Seyfert galaxies (Ho 2008). The mass accretion rate ${\dot{M}}_{\mathrm{acc}}={L}_{\mathrm{bol}}/\eta {c}^{2}$ for this galaxy is ∼0.01 M_⊙ yr⁻¹, assuming a mass-to-energy conversion efficiency factor η = 0.1.

If we assume that all of the radio luminosity is attributable to star formation, a star-formation rate (SFR) can be derived using the relation in Condon (1992). SFR (for stellar masses ≥5 M_⊙) was ∼4.8 M_⊙ yr⁻¹ for the 1.5 GHz VLA emission assuming a spectral index of −0.8. The SFR can also be estimated using the optical emission lines (Kennicutt 1998). The SFR calculated from all three components of the Hα line (=5.0 × 10⁴⁰ erg s⁻¹) turns out to be ∼0.40 M_⊙ yr⁻¹. This SFR is an order of magnitude lower than the SFR derived from the radio flux density on kiloparsec-scales. Therefore, the "radio excess" implied from these data supports the AGN origin of the radio emission in KISSR 1219.

Basic properties of the line emitting gas can be obtained from the measurements of the emission line fluxes. We have used the combined line flux of all the fitted Gaussian components in the following analysis. The [S ii] or [O ii] doublet line intensity ratios gives information about the average electron density n_e of the gas. Using the [S ii] λλ6716, 6731 line intensities (${R}_{[{\rm{S}}{\rm{II}}]}=1.162$) and following the standard calibration of Osterbrock (1989), the electron density of the gas is estimated to be ∼250 cm⁻³ at a temperature of 10⁴ K, consistent with the densities observed in Seyfert galaxies¹⁵ (Ho 1996). The FWHM([O i]) > FWHM([S ii]), is also consistent with the empirical trend seen in Seyferts (Ho 1996). This trend is an indication of the density stratification in the sense that dense material is closer to the center. Since [S ii] probes the low-density region, the density stratification suggests that this gas is located far from the center.

The ratio of Balmer lines Hα/Hβ or Balmer decrement, estimates the attenuation of the emission lines due to dust along the line of sight. In low-density regions of AGN, a value of 3.1 is often used (Osterbrock & Ferland 2006). Any deviation from this value is a signature of dust extinction. The ratio of Hα to Hβ for KISSR 1219 is 6.25 indicating the presence of dust distributed as a foreground screen. The optical extinction A_V toward this galaxy is obtained to be 2.2. We used the galactic extinction law based on Cardelli et al. (1989).

The gas mass of the NLR can be expressed as ${M}_{\mathrm{gas}}=\tfrac{{m}_{p}{L}_{{\rm{H}}\alpha }}{{n}_{e}{j}_{{{\rm{H}}}_{\alpha }}(T)},$ where m_p is the mass of the proton, n_e ∼ 250 cm⁻³ is the electron density, ${j}_{{{\rm{H}}}_{\alpha }}(T)=3.53\times {10}^{-25}$ cm³ erg s⁻¹ is the emission coefficient. The calculated ionized gas mass is ∼4.7 × 10⁵ M_⊙. If we assume that this gas is being pushed out, the dynamical time scale of the gas flowing out with an average outflow velocity of ∼300 km s⁻¹ to the edge of the NLR (in this case, >2.2 kpc, which is the size of the SDSS fiber) is t_dyn = r/v = 7 Myr. The corresponding mass outflow rate thus is ${\dot{M}}_{\mathrm{out}}={M}_{\mathrm{gas}}/{t}_{\mathrm{dyn}}=0.06$ M_⊙ yr⁻¹. The mass outflow rate is an order of magnitude lower than the current SFR in the galaxy. The gas outflow therefore does not seem to have a significant impact on the SFR of the galaxy.

Line ratio diagrams are powerful tools to investigate the physics of emission line gas. Using grids of theoretical models, one can estimate the parameters such as shock velocity, ionization parameter, or chemical abundance. The standard optical diagnostic diagrams include [O iii] λ5007/Hβ versus [N ii] λ6583/Hβ, λ5007/Hβ versus [S ii] λλ6716,6731/Hα, or [O i] λ6000/Hβ λ5007 versus [O iii] λ5007/Hβ (Veilleux & Osterbrock 1987). MAPPINGS III shock and photoionization modeling code has been used to predict the line ratios (Dopita & Sutherland 1996; Allen et al. 2008). We use the IDL Tool for Emission-line Ratio Analysis (ITERA; Groves & Allen 2010) for generating line ratio diagrams. These line ratios were studied for several models including dust-free AGN photoionization, dusty AGN photoionization, shock-only, precursor-only, and shock+precursor models (Groves et al. 2004; Allen et al. 2008).

Shocks can be a powerful source of generating ionizing photons. The kinetic energy of the shock compresses and heats the ambient gas, producing a significant fraction of the EUV and X-ray photon field via free–free emission. While slow moving gas just heats up the gas, fast radiative shocks can generate a stronger radiation field. The photoionization front of this strong radiation field increases rapidly and exceeds that of the shock front, thereby separating itself from the shock front. This photoionization front forms the "precursor" H ii region in the upstream gas. The emission in this precursor region dominates the emission in the shock-only region. A combination of the radiative shocks and the precursor region results in a stronger radiation field and a harder spectrum. Figures 5–7 show the standard optical line ratio diagrams. The line ratios of KISSR 1219 lie in the general area expected for Seyfert galaxies in the BPT diagram. The observed line ratios are reasonably well explained by both dusty AGN photoionization as well as fast radiative shocks (v_s ∼ 300 km s⁻¹, see Figure 5–7) with a precursor, likely driven by the AGN jet.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We note that in Figure 7, the observed [N ii] line ratio appears to not fit the dusty AGN photoionization model. This, however, can be attributed to the higher abundance of N (at least a factor of 2; Osterbrock 1989) and the observed [N ii] line ratio can be brought into agreement with the dusty AGN photoionization model by increasing the N abundance. In the shock+precursor model, however, the strength of the ionizing field does not depend significantly on the atomic abundance. Therefore, both the models are viable for explaining the line ratios in KISSR 1219.