MUSE Reveals Extended Circumnuclear Outflows in the Seyfert 1 NGC 7469

NGC 7469 is a nearby galaxy hosting a Seyfert type 1 AGN with a supermassive black hole with mass log₁₀M_BH = ${7.32}_{-0.10}^{+0.09}$ M_⊙ and bolometric luminosity L_bol = ∼10⁴⁵ erg s⁻¹ (Ponti et al. 2012). It is classified as a luminous infrared galaxy (LIRG) due to the starburst concentrated in its circumnuclear ring, triggered by interaction with the IC 5283 galaxy. The ring outer radius is 25 (∼830 pc) and the inner radius 07 (∼232 pc), with bright knots at ∼15 (∼500 pc; Genzel et al. 1995). This ring contains young (1–20 Myr) massive stars (Diaz-Santos et al. 2007).

Davies et al. (2004) and Izumi et al. (2015) found molecular gas structures in the central region using millimeter observations, including a circumnuclear disk (∼300 pc) that Izumi et al. (2020) revealed to be an X-ray dominated region produced by the AGN.

Blustin et al. (2007) reported X-ray spatially unresolved nuclear winds (warm absorbers) with LoSVs of −580 to −2300 km s⁻¹, later confirmed by Behar et al. (2017) and Mehdipour et al. (2018) at LoSVs of −400 to −1800 km s⁻¹ within a distance of 2–80 pc from the black hole. All these authors reported ultraviolet (UV) counterparts of the X-ray winds.

Müller-Sánchez et al. (2011) found a biconical outflow in the infrared coronal line [Si vi]λ1.96 μm using the Very Large Telescope (VLT)/Spectrograph for INtegral Field Observations in the Near Infrared (SINFONI) and Keck/Osiris IFS, extending up to 380 ± 25 pc from the AGN, with a maximum velocity of ∼130 km s⁻¹ at 220 ± 25 pc.

While this outflow was not reported in recent optical IFS observations with the Gran Telescopio CANARIAS (GTC)/Multi-Espectrógrafo en GTC de Alta Resolución para Astronomía (MEGARA) by Cazzoli et al. (2020), they described a non-rotational turbulent component that might be associated with it.

We report, for the first time, optical extended ionized outflows in the circumnuclear region of NGC 7469, based on the kinematics of the [O iii]λ5007 and Hβ emission lines. The comoving distance to NGC 7469 is 69.64 Mpc¹¹ (z = 0.01632, Keel 1996), adopting a standard ΛCDM cosmology (H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_Λ = 0.7 and Ω_M = 0.3).

NGC 7469 was observed in 2014 August 19, during the science verification run of the Multi-unit Spectroscopic Explorer (MUSE; Bacon et al. 2004) IFS instrument at the VLT of the European Southern Observatory (ESO; Chile). The pilot study of the All-weather MUse Supernova Integral-field of Nearby Galaxies survey (AMUSING; Galbany et al. 2016) and the AMUSING++ compilation¹² (López-Cobá et al. 2020) include these data. MUSE covers a field of view of $1\ {\mathrm{arcmin}}^{2}$ with a spatial sampling of 02 per spaxel (top-left panel, Figure 1). For NGC 7469 each spaxel presents a scale of 66.44 pc ($332.2\ \mathrm{pc}\,{\mathrm{arcsec}}^{-1}$). The seeing had a FWHM = 123 (∼409 pc). Data reduction followed the standard procedures, using the Reflex (Freudling et al. 2013) package and the MUSE pipeline (Weilbacher et al. 2014). We corrected for systemic velocity using the value by Keel (1996) (4898 ± 5 km s⁻¹), obtained by a combination of emission-line methods due to the absence of absorption lines in the central region.

Zoom In Zoom Out Reset image size

Beam smearing—the scattering of light from the spatially unresolved broad-line region (BLR) and inner narrow-line region (NLR) due to seeing—can lead to overestimation of the size of the extended narrow-line region (ENLR) and of the velocity of associated outflows in type 1 AGN (Husemann et al. 2016). This effect was corrected through the QDeblend3D software described in Husemann et al. (2012, 2014, 2016). We assume a surface brightness model for the host galaxy, fixing brightness, effective radius, Sérsic index, and axis ratio, with values from Bentz et al. (2009). QDeblend3D produces two data cubes (Figure 1): one contains the AGN continuum and emission from the BLR and the NLR (hereafter "AGN data cube"); the other contains the host galaxy stellar continuum, SF regions, and the ENLR (hereafter "Host data cube").

We subtract a synthetic stellar continuum template from each spectrum in the Host data cube, to obtain "pure emission" spectra. We create the templates (after correcting for Galactic extinction using dust maps produced by Schlegel et al. 1998) using the Starlight population synthesis code (Cid Fernandes et al. 2005). We use the base set of 150 simple stellar populations¹³ selected by Asari et al. (2007)—25 ages (1 × 10⁶–1.8 × 10¹⁰ yr) and six metallicities (0.005–2.5 Z_⊙)—with the Chabrier (2003) initial mass function. Emission lines are masked out. Intrinsic extinction is corrected in the process using the Cardelli et al. (1989) reddening law (R_V = 3.1).

We study the line profiles of [O iii]λλ4959,5007 and Hβ in a box 5'' (1.611 kpc) wide, centered on the AGN, that contains the SF ring and corresponds to the field studied by Diaz-Santos et al. (2007; Figure 1). The line profiles of many spaxels are asymmetric and broadened at their bases, suggesting the presence of winds. We characterize these features using two approaches (see below).

The results of [O iii]λ4959 and [O iii]λ5007 are consistent, including the expected one-third flux ratio. We only report here the results for [O iii]λ5007 (hereafter [O iii]), due to its higher S/N.

Following Harrison et al. (2014), each emission line is fitted with three Gaussian components (Figure 2(a)), whose sum creates a synthetic line profile. To fit the Gaussians we use the Levenberg–Marquardt algorithm (Marquardt 1963), implemented with the IDL MPFIT libraries (Markwardt 2009). Spectra were interpolated, resulting in ≈50 data points per emission line, versus ≈10 in the original data. This increases the chance of obtaining a good solution by decreasing the solution space, and reduces computing time. We reject all Gaussian components with S/N < 3.

Zoom In Zoom Out Reset image size

We avoid physically interpreting the Gaussians. Instead we characterize the kinematics by deriving quantities from the cumulative function of the whole synthetic line profile: the offset velocity Δv (the mean of the velocities at the 5th and 95th percentiles) the width at 80% of the flux W₈₀, and the FWHM. Hence, Δv measures the asymmetry of the line profile related with gas motion on the line of sight; the sign indicates the direction. The W₈₀/FWHM ratio is the relative broadening at the base of the line.

López-Cobá et al. (2020) used an alternative non-parametric approach; however, the current approach is accurate enough for our goal.

We fit two Gaussian components to the [O iii] line (Figure 2(b)) following the approach by Woo et al. (2016) and Karouzos et al. (2016). The Gaussian component closest to the rest-frame velocity (hereafter, central component), is related with gas dominated by the host galaxy gravitational potential, rotating in the galactic disk; the second Gaussian component accounts for the outflowing gas as a whole (hereafter, blueshifted component; see Figure 2(b)).

We apply the same fitting method described in Section 2.4. We measure the kinematic parameters of each Gaussian component: the LoSV from the Doppler shift of the Gaussian peak with respect to the rest frame, and the velocity dispersion (σ) from the FWHM—corrected for the instrumental width (FWHM_inst ∼ 158 km s⁻¹ for [O iii]).

We estimate uncertainties through Monte Carlo simulations, following Lenz & Ayres (1992), iterating 1000 times. For the non-parametric approach, the mean uncertainties across the studied field are ∼3.9, ∼6.7, and ∼6.6 km s⁻¹ for Δv, W₈₀ and the FWHM of [O iii], respectively. For Hβ the corresponding values are ∼2.9, ∼6.5, and ∼6.6 km s⁻¹.

For the two-Gaussian components approach, the mean uncertainties in LoSV and σ are, respectively, ∼60 and ∼32 km s⁻¹ for the [O iii] blueshifted component and ∼30 and ∼6.6 km s⁻¹ for the central component.

Figure 3 shows the maps of Δv, W₈₀ and W₈₀/FWHM for [O iii] and Hβ. As a reference for the position of the AGN and the massive SF regions of the ring, the maps are overlapped to the Hubble Space Telescope (HST) Advanced Camera for Surveys (ACS) F330W near-UV image.¹⁴ Projections of the edges of the [Si vi]λ1.96 μm outflow by Müller-Sánchez et al. (2011) are shown, with the blueshifted cone pointing west and the redshifted cone pointing east.

Zoom In Zoom Out Reset image size

The [O iii] results reveal the existence of two outflow kinematic regimes, located in different regions. In Figure 3(a), in a region labeled as "A," a strong [O iii] blueshifted asymmetry (Δv up to −310 km s⁻¹, reaching ∼531 pc (∼16) to the north) extends northwest of the center, in between the AGN and the massive SF regions of the ring. Approximately half of it extends within the limits of the projected west [Si vi]λ1.96 μm cone, while the rest lies outside to the north.

This region also presents a broadening of W₈₀ > 600 km s⁻¹ and W₈₀/FWHM > 3 (Figures 3(b) and (c)). This high asymmetry and broadening are strong kinematic evidence of the presence of an outflow in the "A" region.

The existence of a second kinematic regime in the rest of the ring and the inner regions is disclosed by less prominent asymmetry (−200 ≤ Δv ≤ −100 km s⁻¹) and broadening (250 < W₈₀ < 500 km s⁻¹ and 1.5 < W₈₀/FWHM < 3). We ignore here the outer regions of the field due to their lower S/N.

The Hβ maps show only one outflow regime, similar to the slowest one observed in [O iii]. The Δv, W₈₀, and W₈₀/FWHM maps are shown in Figures 3(d), (e), and (f), respectively. Outside the "A" region, Hβ has similar but less pronounced Δv than [O iii], with most spaxels having Δv > −200 km s⁻¹.

Low W₈₀/FWHM < 2 ratios dominate the field except for the western SF ring. Despite reaching W₈₀ > 300 km s⁻¹ (with a maximum of ∼640 km s⁻¹), they keep W₈₀/FWHM < 3. High FWHM values keep low W₈₀/FWHM ratios, possibly due to the presence of turbulence or shocks, not only in the wind but also in the bulk of the gas, affecting the whole line profiles.

Figure 4(a) shows the LoSV map of the [O iii] blueshifted component. The two outflow regimes found in Section 3.1 are confirmed, with different LoSV ranges of the blueshifted component: high velocities (LoSV < −400 km s⁻¹) are found in the "A" region; lower velocities (LoSV > −400 km s⁻¹) are found around the rest of the field. The high LoSV spaxels in Figure 4(a), cover a smaller area than that of the high Δv spaxels in Figure 3(a). The fastest blueshifted component (shown at Figure 2(b)) reaches LoSV = −715 km s⁻¹, at ∼240 pc (∼07) northwest of the center. Velocity dispersion values of σ (>300 km s⁻¹) do appear in the the "A" region (Figure 4(b)), similar to W₈₀, but they extend further over the western massive SF regions.

Zoom In Zoom Out Reset image size

The velocity–velocity dispersion (VVD) diagram (Karouzos et al. 2016; Woo et al. 2016) provides a straightforward visualization of the different kinematic regimes (Figure 4(c)). Once the central and blueshifted components are plotted in the VVD, the former is clearly is separated from the latter in the parameter space. Following the interpretation of Woo et al. (2016) and Karouzos et al. (2016), the more extreme kinematics of the blueshifted component could be evidence of outflows (not dominated by the gravitational potential of the host galaxy), while the positive and negative LoSV values of the central component suggest motion in a galactic disk.

The high LoSV spaxels in the "A" region of Figure 4(a) (with a median LoSV of −581 km s⁻¹ and median σ of 300 km s⁻¹) are clearly separated from the bulk of the blueshifted components (with median LoSV and σ of −284 and 188 km s⁻¹, respectively), which is consistent with the existence of the two kinematic regimes in the blueshifted component.

Note that if the mean LoSV of the central component was assumed as systemic reference frame, the outflow velocity would be 134 km s⁻¹ slower. However, the existence of both outflow regimes would still hold solidly, with the slower gas mean LoSV ∼ −150 km s⁻¹. In fact, more than one Gaussian is needed to fit the line profiles, as shown by the non-parametric approach.

The BPT diagnostic diagram (Baldwin et al. 1981) informs on the excitation mechanisms, provided that only Gaussian components with similar kinematics, i.e., tracing the same bulk of gas, are considered. The blueshifted component kinematics of [O iii] and Hβ are inconsistent for many spaxels, and shall not be combined in the same diagram. Due to blending of the shifted components, we could unambiguously identify only the peaks of the central components in the Hα–[N ii] complex and use them for BPT diagnostics.

Figure 5 shows the BPT-NII diagnostics—based on the [O iii]λ5007/Hβ and [N ii]λ6584/Hα line ratios—of the central component. The mean flux uncertainties (in ×10⁻¹⁷ erg cm⁻² s⁻¹ units) are ∼0.58 for [O iii], ∼0.52 for Hβ, ∼0.76 for [N ii]λ6584, and ∼3.4 for Hα.

Zoom In Zoom Out Reset image size

AGN excitation dominates the northeast quadrant up to ∼500 pc from the center, while a combination of SF and transition-object-like (TO) excitation extends across the rest of the field (low ionization nuclear emission line region (LINER)-like excitation is scarce). The "A" region with the fastest [O iii] outflow presents a combination of AGN, TO and SF excitation, the latter over the massive SF regions of the ring. We use this result as exploratory; direct extension to the blueshifted component should be avoided. More detailed analysis will be presented in a future paper.