Discovery of a Spatially and Kinematically Resolved 55 kpc Scale Superbubble Inflated by an Intermediate-redshift Non-BAL Quasar

In order to understand the dynamics and main properties of the ionized outflows, we perform a kinematical analysis of the forbidden lines. We remove the BEL from the quasar nucleus, which is represented by double Gaussians in the spatial regions where the narrow emission lines are negligible. With the quasar contribution removed, we scrutinize the spatially resolved emission lines following the method described in Zhao et al. (2021).

We perform a two-step spectral fit to delineate the gas kinematics. First, we extract the spectrum in each spaxel and subtract the continuum using the interpolation method from wavelength on two sides of the [O iii] line, where the continuum is free of any line emission and artifacts. We use the Fe ii template from Tsuzuki et al. (2006) to subtract the Fe ii emission. Second, we assume that the [O iii] doublet is originating from the same upper level, and the intensity ratio is I(5007)/I(4959) ∼ 3. The lines are fitted with the same central velocity and velocity dispersion. The profile of the [O iii] λ5007 Å emission line is generally complex, so the [O iii] λ λ4959, 5007 doublet is fitted with a combination of multiple Gaussians by minimizing χ² using the Python package MPFIT. We fit the [O iii] profiles to no more than three Gaussians, following Liu et al. (2014).

The [O iii] nebula surrounding HE 0238−1904 is spatially resolved by our IFS observation. The [O iii] map of HE 0238−1904 is shown in Figure 1(a), where the false color is used to represent the intensity on a logarithmic scale. The surface brightness sensitivity (rms noise) of our [O iii] maps is approximately σ ∼ 3 × 10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻². We use a 1.5σ threshold to create this map.

Zoom In Zoom Out Reset image size

We obtain nonparametric measurements of the emission line profiles following the method described in Liu et al. (2013a). These include the following.

• 1.  
  Zeroth-moment map: intensity map of the [O iii] λ5007 line.
• 2.  
  First-moment map: v_med, median velocity map.
• 3.  
  Second-moment map: line width map, W₈₀, the velocity width of the line that encloses 80% of the total flux. This is defined as the difference between the velocities at 10% and 90% of cumulative flux: W₈₀ = v₉₀ − v₁₀.
• 4.  
  Asymmetry: $A=\tfrac{({v}_{90}-{v}_{\mathrm{med}})-({v}_{\mathrm{med}}-{v}_{10})}{{W}_{80}}$. With our definition, a profile with a significantly blueshifted wing has a negative A value, while a symmetric profile has A = 0.

We perform these nonparametric measurements on the best-fitting profiles. Figure 1 shows these parameters of the ionized gas derived from the fit of the [O iii] λ5007 line. The maps are obtained by selecting only pixels with a signal-to-noise ratio (S/N) ≥ 1.5.

The forbidden emission line [O iii] λ λ5007, 4969 doublet is adopted as a tracer of ionized outflows on large scales. Figure 1 shows the [O iii] intensity, velocity, velocity dispersion, and asymmetry maps. Our IFS data confirmed that the [O iii] λ5007 line is spatially and kinematically resolved. In the following, we discuss the velocity, velocity dispersion, and asymmetry of the ionized gas and compare them with previous works (e.g., Arav et al. 2013; Liu et al. 2013a; Carniani et al. 2015).

The velocity field of the [O iii] nebula is remarkably well organized (Figure 1(b)). Only the blue-side outflows are detected, whereas the red side is probably obscured by the host galaxy along the line of sight. For this reason, the [O iii] line profile is asymmetric, with a prominent blueshifted wing (Figures 1(d) and 2).

Zoom In Zoom Out Reset image size

We find the maximum W₈₀ value of ∼1600 km s⁻¹ in the center, which is comparable to that of known quasar outflows (e.g., Liu et al. 2013a; Carniani et al. 2015; Zakamska et al. 2016; Kubo et al. 2022) but considerably larger than the usual narrow lines in the quasar (Lonsdale et al. 1993). The W₈₀ value drops to ∼200 km s⁻¹ in the outer region. The asymmetry map shows regions with heavy blueshifted wings (<0) that are spatially associated with a high velocity dispersion (∼1600 km s⁻¹; see Figure 1(c)).

In Figure 1(d), we show the map of the asymmetry parameter A. The asymmetry parameter A is uniformly negative in the bright central part of HE 0238−1904, indicating heavily blueshifted wings in the line profiles. This is the telltale signature of an outflow that may be proceeding in a symmetric fashion but whose redshifted part is obscured by the material in the host galaxy or near the nucleus (Whittle 1985). In the fainter outer region where the peak S/N of even the brightest emission line [O iii] is just a few, typically one Gaussian component is sufficient to fit the line profile. Therefore, the asymmetry parameter tends to be at zero values (A = 0).

Motivated by the different kinematic components apparent in the velocity map, we extract and fit the [O iii] lines in 12 regions. Figures 1(b) and 2 show the position and spectra of these extraction regions, referred to as A–L. Regions A and B correspond to the position of the quasar itself and its immediate vicinity, while C–L correspond to the outer regions. The spectral fits reveal that the width of the line profile in the center is broader than that of the outer regions. The most significant asymmetry in the line profile is present in the center. Consequently, the smooth morphology of the [O iii] nebula, the velocity, the high velocity dispersions of the gas, and the blueshifted asymmetry map all suggest that we have detected ionized outflowing gas in HE 0238−1904.

In Figure 3, we show values of W₈₀ in all spaxels as a function of projected distance from the center. The radial profiles of W₈₀ are almost flat at projected distances R ≲ 8 kpc and appear to decrease at larger radii. This is different from previous results found in other quasars based on long-slit and IFS observations, which reported flat W₈₀ profiles (Greene et al. 2011; Liu et al. 2013a). One of the concerns with the measured decline in W₈₀ is that the broad component can no longer be identified in the outer parts, and the W₈₀ measurement is due to the narrow component, hence the declining W₈₀. In this case, the W₈₀ parameter would be almost constant across the nebulae in most cases, perhaps declining slightly toward the outer parts. This is not consistent with our W₈₀ profile. The possible origin of the rapid decline in W₈₀ is discussed in Section 5.1.

Zoom In Zoom Out Reset image size

The ionized gas extends to the southeast of the nucleus (see Figure 1), where three blobs (regions II–IV; Figure 1(a)) and a nuclear region (region I) are present. Both the velocity and the velocity dispersion of regions II–IV are similar (v ∼ 700 km s⁻¹, W₈₀ ∼ 400 km s⁻¹).

The electron density can be derived from [O ii] I(3729)/I(3726) ratios. For this purpose, we stack the spectra taken from regions I–IV and fit the [O ii] doublet emission lines by fixing the kinematics of the two lines (Figure 4). Assuming an electron temperature of 10,000 K, we estimate the electron density to be n_e = 422 ± 355, 208 ± 54, and 264 ± 99 cm⁻³ in regions I, III, and IV, respectively. The electron density in region II is n_e ≤ 155 cm⁻³, which is the upper limit at 3σ significance. This result is close to that of other quasar outflows (a few hundred per cubic centimeter; Nesvadba et al. 2006, 2008). Note that the [O ii] doublet is blended together in the center and difficult to fit, so we stack and fit the spectra from the relatively outer region around the nucleus (region I). Hence, the actual electron density can be higher in the nuclear region.

Zoom In Zoom Out Reset image size

There are two possibilities for the origin of the ionized gas. The gas can be either a nuclear outflow launched by AGN or a tidal tail due to galaxy interaction. First, the extended ionized gas is ubiquitous in quasars and generally driven by radiation pressure, so the ionized gas in HE 0238−1904 could well be the same. On the other hand, one could suspect that this spatially extended feature with a low velocity dispersion (regions II–IV) is in fact a tidal tail resulting from interaction with a lower-mass companion. We analyze all spectra of the companion galaxies around HE 0238−1904 in the field of view, and the difference between galaxy systemic velocity and this extended emission line gas velocity is found to be at least several hundred kilometers per second. This excludes the possibility of ionized gas being a tidal tail because the velocity difference between companion galaxies and tidal tails is generally ≲100 km s⁻¹ (e.g., Fu et al. 2021).

In addition, the intensity ratios of the emission lines facilitate our analysis of the physical conditions of the ionized gas. Specifically, we use the [O iii]/Hβ ratios to quantify the degree of ionization. To obtain higher S/Ns, we stack the spectra from regions II–IV shown in Figure 1(a). The Hβ and [O iii] spectra and the [O iii]/Hβ ratios in the extended regions (II–IV) are shown in Figure 5. In this object, Hβ almost follows the same spatial distribution as that of [O iii] because [O iii]/Hβ is almost a constant. The ratio reveals that [O iii]/Hβ is close to 10 in these three regions, implying a high-ionization state in general. Based on the BPT diagram (Baldwin et al. 1981), [O iii]/Hβ > 10 implies AGN dominance. The [O iii]/Hβ profile indicates that the AGN plays a dominant role in the ionization of the large-scale gas.

Zoom In Zoom Out Reset image size

In Figure 3, we find that the W₈₀ parameter is almost flat at projected distances R ≲ 8 kpc and appears to decrease rapidly at larger radii. Here we discuss the possible origins of the rapid decline of W₈₀ with increasing distance from the center. This decline of W₈₀ implying an apparent narrowing of the line profile is likely due to the fact that the outflow becomes more directional or collimated within these structures. We thus confirm the previous results based on UV analyses that indicated collimated outflow (Muzahid et al. 2012).

There are various mechanisms capable of establishing an apparently declining W₈₀ profile. First, in the central region, the wind expands in all directions from the quasar, whereas at larger distances, the opening angle of the outflow is decreased, perhaps because there are low-density regions along which the wind prefers to propagate. Also, the outflow in the central region may be experiencing large turbulent motions due to interaction with the ISM, but once they escape out of the inner galaxies, the flow becomes more organized and mostly radial. The episodic quasar outbursts may drive a shock wave through the ISM of the galaxy and clear out some of it (Novak et al. 2011). In any subsequent episodes, the wind suffers less resistance from the ISM and breaks out in these directions, producing the large-scale bubbles (Faucher-Giguère et al. 2012). Numerical simulations show that this phenomenon is expected for jet-driven winds (Sutherland & Bicknell 2007). We detected three spatially resolved blobs, which are likely part of the rim of the superbubble. In addition, the outflowing gas is launched and accelerated somewhere close to the quasar and proceeds ballistically. Thus, these outflowing gaseous blobs eventually slow down as they overcome the potential well of the host galaxy, producing a decline in W₈₀.

The ionized outflow traced by [O iii] emission extends to 50 kpc and beyond, and W₈₀ remains larger than 500 km s⁻¹ up to 20 kpc (Figure 3). Combined with the large-scale morphology and the W₈₀ profile, the extended gas we detected closely resembles a superbubble. Future high-quality soft X-ray observations are needed to fully investigate the origin of this ionized gas.

We find that the ionized gas southwest of the nucleus extends to a projected distance reaching 55 kpc (Figure 1). Some characteristic parameters of the ionized gas are further estimated. We adopt a blueshifted velocity of v₀ ∼ 690 km s⁻¹ and a distance to the galaxy nucleus of R₀ ∼ 55 kpc. We estimate a dynamical timescale to be t ∼ 7.8 × 10⁷ yr, i.e., the time required for the gas from the nuclear region to reach such a distance with an average velocity of v₀. The mass of the gas can typically be estimated using Hβ luminosity L_Hβ and electron density n_e (e.g Liu et al. 2013a; Harrison et al. 2014). The total mass of this ionized extended gas can be derived as

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{\mathrm{gas}}}{2.82\times {10}^{9}{M}_{\odot }}=\left(\displaystyle \frac{{L}_{{\rm{H}}\beta }}{{10}^{43}\,\mathrm{erg}\,{{\rm{s}}}^{-1}}\right){\left(\displaystyle \frac{{n}_{e}}{100\,{\mathrm{cm}}^{-3}}\right)}^{-1}.\end{eqnarray} \tag{ 1 }$$

We take the electron density of n_e ∼ 400 cm⁻³ in region I because the bulk of the ionized gas is concentrated in the central region. We find L_Hβ ∼ 3.02 × 10⁴² erg s⁻¹ and M_gas ∼ 2.1 × 10⁸ M_⊙. Combining with the average velocity of v ∼ 690 km s⁻¹, the total kinetic energy of the gas can be estimated as

$$\begin{eqnarray}&&{E}_{\mathrm{kin}}=\displaystyle \frac{1}{2}{M}_{\mathrm{gas}}{v}_{\mathrm{gas}}^{2}=1.0\times {10}^{57}\,\mathrm{erg}.\end{eqnarray} \tag{ 2 }$$

We can also estimate that the mass outflow rate $\dot{M}$ is 2.7 M_⊙ yr⁻¹, and the kinetic energy rate is ${\dot{E}}_{\mathrm{kin}}\sim$ 4.1 × 10⁴¹ erg s⁻¹. The momentum flux of the outflow (${\dot{P}}_{\mathrm{outflow}}=\dot{M}\times v$) is 1.2 × 10³⁴ dynes, or $\mathrm{log}(c{\dot{P}}_{\mathrm{outflow}}/{L}_{\odot })=10.9$. The ratio of the momentum flux of the outflow to the AGN radiation is only 0.002.

For effective feedback, theoretical modeling predicts that ∼0.5%–5% of the AGN's bolometric luminosity is converted to the kinetic energy of the outflow (e.g., Hopkins & Elvis 2010). The kinetic luminosity of the outflow detected at 55 kpc is insignificant, only 0.0002% of its bolometric luminosity (L_bol = 10^47.2 erg s⁻¹; Arav et al. 2013). This is much lower than the previously estimated ${\dot{E}}_{\mathrm{kin}}/{L}_{\mathrm{bol}}\sim 1 \%$ from absorption on a 3 kpc scale (Arav et al. 2013) or the ∼0.05%–0.1% estimated for a sample of z ∼ 2.4 quasars (Carniani et al. 2015). Carniani et al. (2015) fit a log-linear relation between $\dot{M}$ and L_bol for the ionized outflows detected on kiloparsec scales in their quasar sample. Adopting their Equation (10), the expected $\dot{M}\sim 30{M}_{\odot }$ yr⁻¹ for HE 0238−1904 agrees fairly well with the derived $\dot{M}\sim 40{M}_{\odot }$ yr⁻¹ in Arav et al. (2013). In contrast, our estimated mass outflow rate at large scales (55 kpc) in HE 0238−1904 is an order of magnitude lower. The kinetic power of the outflow (≪0.1% L_bol) implies that it is no longer an important contributor to feedback at this scale.

Using the HST/COS and FUSE UV spectra of HE 0238–1904, Muzahid et al. (2012) reported a detection of outflowing gas in multiple absorption lines: O vi, Ne viii, and Mg x. They identified evidence for a similar covering factor in several absorption components that kinematically spread over ∼1800 km s⁻¹, suggesting a collimated outflow. This is consistent with our expectation from the radial variation of W₈₀.

Determining the location of the absorption outflow is usually challenging. Nevertheless, assuming a spherical geometry, Muzahid et al. (2012) constrained the radial distance of the absorbing gas to be R ∼ 90 pc based on the photoionization modeling result. The detailed UV absorption line analysis by Arav et al. (2013) was able to robustly derive a distance of R ∼ 3 kpc from the outflow to the nucleus using absorption troughs of O iv and O iv*. As revealed by the MUSE data, the projected distance of the outflow from the AGN is R ∼ 55 kpc for the blueshifted side. This is significantly further from the previous locations of the absorbing gas. In addition, both the electron density and the kinematics of these outflows determined from IFS of emission lines are different from those derived from previous UV absorption line analyses (e.g., Arav et al. 2013). The velocity of the absorption line (∼5000 km s⁻¹) is much higher than the velocity of the emission line (∼700 km s⁻¹). This strongly indicates that the outflows detected using the absorption and emission lines are clearly not the same component but likely stratified components of different spatial scale and velocity in the ionized phase outflow.