Quasar feedback and the origin of radio emission in radio-quiet quasars

Abstract

1 INTRODUCTION

Black hole feedback – the strong interaction between the energy output of supermassive black holes and their surrounding environments – is routinely invoked to explain the absence of overly luminous galaxies, the black hole versus bulge correlations and the similarity of black hole accretion and star formation histories (Tabor & Binney 1993; Silk & Rees 1998; Springel, Di Matteo & Hernquist 2005; Hopkins et al. 2006). After years of intense observational effort, specific examples of black hole-driven winds have now been identified using a variety of observational techniques, both at low and at high redshifts (Nesvadba et al. 2006, 2008; Arav et al. 2008; Moe et al. 2009; Alexander et al. 2010; Dunn et al. 2010; Fabian 2012; Harrison et al. 2012).

How these outflows are launched near the black hole and established over the entire host galaxy remains a topic of active research. In particular, it is becoming clear that as the outflow impacts an inhomogeneous interstellar medium of the galaxy, the winds are expected to contain gas at a wide range of physical conditions, and these different ‘phases’ of the winds require different types of observations (Veilleux, Cecil & Bland-Hawthorn 2005). As a result, the determination of the physical parameters of these outflows – including such basic parameters as the mass, the momentum and the kinetic energies they carry – remains challenging. In the last few years, a lot of progress in this area has been made by observations of the coldest components of the outflows which are in the form of neutral or even molecular gas (Feruglio et al. 2010; Fischer et al. 2010; Rupke & Veilleux 2011; Sturm et al. 2011; Aalto et al. 2012; Rupke & Veilleux 2013a,b; Veilleux et al. 2013; Cicone et al. 2014; Sun et al. 2014). It remains unclear how common this cold component is in winds driven by a powerful active nucleus and how the mass, momentum and energy carried by the wind are distributed across the different phases.

Several years ago, we embarked on an observational programme to determine whether radio-quiet, luminous quasars have observable effects on their galaxy-wide environment. One of our lines of investigation is to determine the extent and kinematics of the warm (T ∼ 10⁴ K) ionized gas – the so-called narrow-line region of quasars. We simplify the observational task by looking at obscured quasars (Zakamska et al. 2003; Reyes et al. 2008) – those where the line of sight to the nucleus is blocked by intervening material, allowing us to study the distribution of matter in the galaxy unimpeded by the bright central source. We find strong evidence that ionized gas is extended over scales comparable to or exceeding that of the host galaxy; furthermore, it is kinematically disturbed and is not in equilibrium with the gravitational potential of the galaxy (Greene et al. 2009, 2011; Hainline et al. 2013).

More recently, we surveyed a sample of obscured radio-quiet quasars using a spectroscopic integral-field unit (Liu et al. 2013a,b). We found extended ionized gas encompassing the entire host galaxy (median diameter of nebulae of 28 kpc), suggestive of wide-angle outflows, and determined kinetic energies of these outflows to be well in excess of 10⁴⁴ erg s⁻¹, with a median 2 per cent conversion rate from the bolometric luminosity to the kinetic energy of warm ionized gas; more energy can be carried by other components. Furthermore, we identified several candidate objects where the wind has ‘broken out’ of the denser regions of the galaxy and is now expanding into the intergalactic medium, sometimes in bubble-like structures (Greene, Zakamska & Smith 2012; Liu et al. 2013b).

These observations demonstrate the presence of extended ionized gas in host galaxies of type 2 quasars, which is apparently out of dynamical equilibrium with the host galaxy and is likely in an outflow on the way out of the host. In this paper, we examine spectra of several hundred obscured quasars and we study the relationships between gas kinematics and other physical properties of these objects. In Section 2, we describe the sample selection, the data set and the measurements. In Section 3, we conduct kinematic analysis of the optical emission lines. In Section 4, we discuss the relationships between multiwavelength properties of quasars and kinematic measures of their ionized gas nebulae. In Section 5, we present composite spectra and discuss trends in weak emission lines. We present qualitative models for radio emission and emission lines in Section 6, and we summarize in Section 7.

We use an h = 0.7, Ω_m = 0.3, Ω_Λ = 0.7 cosmology throughout this paper. Sloan Digital Sky Survey (SDSS) uses vacuum wavelengths, but for consistency with previous literature we use air wavelengths in angstroms to designate emission lines. Wavelengths are obtained from NIST (Kramida et al. 2013) and Atomic Line List¹ and converted between air and vacuum as necessary using Morton (1991). Objects are identified as SDSS Jhhmm+ddmm, with full coordinates given in the catalogue by Reyes et al. (2008). We use ‘1D’, ‘2D’ and ‘3D’ abbreviations for one-, two-, and three-dimensional values.

2 DATA AND MEASUREMENTS

2.1 Sample selection and host galaxy subtraction

The obscured quasar candidates studied here were selected from the spectroscopic data of the SDSS (York et al. 2000) based on their emission-line ratios and widths to be the luminous analogues of Seyfert 2 galaxies (Zakamska et al. 2003); the most recent sample contains 887 objects at z < 0.8 (Reyes et al. 2008). Infrared observations demonstrate that these sources have high bolometric luminosities (up to 10⁴⁷ erg s⁻¹; Zakamska et al. 2004, 2008; Liu et al. 2009). Chandra and XMM–Newton observations show that they contain luminous X-ray sources with large amounts of obscuration along the line of sight (Zakamska et al. 2004; Ptak et al. 2006; Vignali et al. 2010; Jia et al. 2013). Hubble Space Telescope imaging and ground-based spectropolarimetry demonstrate the presence of scattered light – a classical signature of a buried broad-line active nucleus (Antonucci & Miller 1985; Zakamska et al. 2005, 2006). In other words, all follow-up observations thus far are consistent with these objects being luminous obscured quasars. Our estimate of the number density of these objects suggests that they are at least as common as unobscured quasars at the same redshifts and line luminosities (Reyes et al. 2008).

In this paper, we examine the kinematic structure of narrow emission-line gas in 568 objects (out of the entire sample of 887 by Reyes et al. 2008) selected to have [O iii] luminosities above 10^8.5 L_⊙. Their distribution in [O iii] luminosity/redshift space is shown in Fig. 1. At the median redshift of the sample presented here z = 0.397, the SDSS fibre (3 arcsec in diameter) covers the galaxies out to 8 kpc away from the centre.

Figure 1.

The distribution of the entire Reyes et al. (2008) sample of type 2 quasars in the redshift/[O iii] luminosity space (grey) and the 568 objects with kinematic analysis in this paper (black).

Open in new tabDownload slide

Ideally, we would like to measure the kinematics of the ionized gas relative to the host galaxy potential. The SDSS spectroscopic pipeline provides a high-quality redshift based on fits of observed spectra to a variety of library templates. In the cases of our objects, this pipeline latches on to the strong emission lines, so the redshifts are affected by narrow-line kinematics and may be offset from the redshifts of the host galaxies. Thus, our first step is to determine the host galaxy redshifts based on the absorption features produced in stellar photospheres.

Even though the quasars are obscured, the continuum from the host galaxy is very difficult to detect. One component of the continuum is due to the stars in the host galaxy. Furthermore, while the direct emission from the quasar is completely blocked in most cases, some quasar light reaches the observer after scattering off of the interstellar medium of the host galaxy and becomes an important continuum contribution at high quasar luminosities. Finally, in very luminous cases the Balmer continuum produced by free–bound transitions in the extended narrow emission-line gas is also seen (Zakamska et al. 2005), and the emission lines tend to confuse the search for stellar features in the continuum.

We use the stellar velocity dispersion code described in detail in Greene & Ho (2006) and Greene et al. (2009) to model the host galaxy continuum and to establish the systemic velocity. The continuum of each quasar is modelled as the linear combination of three stellar models plus a power-law component to mimic a possible scattered light contribution (Zakamska et al. 2006; Liu et al. 2009). For templates, we use Bruzual & Charlot (2003) single stellar population models rather than individual stellar spectra.

We first shift the spectra into the approximate rest frame by using the SDSS pipeline redshifts. Then we fit the continuum over the wavelength range 3680–5450 Å, allowing for the stellar models to have a velocity of up to 300 km s⁻¹ relative to the SDSS frame and to be broadened with a Gaussian function that represents the stellar velocity dispersion of the host. When the host galaxy is well detected, typical stellar absorption features visible in the spectra include Ca H+K, G band and Mg Ib lines. After this procedure, the entire host continuum is subtracted and the spectrum is shifted into the fine-tuned host galaxy frame.

We are able to identify some host galaxy features in 271 objects. In most of these cases, the absorption features are so weak that we do not consider the reported host velocity dispersions to be accurate. The greater benefit of the host subtraction procedure is that in these 271 objects we can analyse narrow-line kinematics relative to an accurately determined host frame. In the remaining objects we find no evidence for stellar features, so we subtract a featureless continuum. For the majority of objects, our workable wavelength range covers [O ii] λλ3726, 3729, [O iii] λ5007 and everything in between.

2.2 Fitting functions and non-parametric measurements

We aim to use non-parametric measures that do not depend strongly on the specific fitting procedure. We need robust measures or robust analogues of the first four moments of the line profile: typical average velocity, velocity dispersion, and the skewness and the kurtosis of the velocity distribution. We fit the profiles with one to three Gaussian components in velocity space, but in principle other fitting functions could be used. We assign no particular physical significance to any of the parameters of the individual components; rather, the goal is to obtain a noiseless approximation to the velocity profile.

We use relative change in reduced χ² values to evaluate which fit should be accepted; if adding an extra Gaussian component leads to a decrease in χ² of <10 per cent, we accept the fit with a smaller number of components. The single-Gaussian fit is accepted for 36 objects, a two-Gaussian fit is accepted for 132 and the remaining 400 objects are fitted with three Gaussians. Almost all objects that have high signal-to-noise ratio (S/N) observations have reduced χ² values that are too high to be statistically acceptable and thus would require either a larger number of components or different fitting functions to be fitted to statistical perfection. Fortunately, the non-parametric measures that we derive are rather robust: adding the third Gaussian component changes our second moment measure w₈₀ by less than 10 per cent in 83 per cent of objects. Examples of line fits are shown in Fig. 2. The objects selected for this figure have the top 10 highest values of w₈₀ (our analogue of the velocity dispersion defined below).

Figure 2.

Spectra of the [O iii] λλ4959, 5007 doublet in the 10 objects with the highest w₈₀ values (2314 ≤ w₈₀ ≤ 2918 km s⁻¹), with their multi-Gaussian fits. The two lines in the doublet are fit simultaneously, under the assumption that the kinematic structure of both lines is the same and that the ratio of amplitudes is 0.337. Dashed lines show the positions of v₁₀, v₅₀ and v₉₀.

Open in new tabDownload slide

Armed with fitting functions performed in velocity space f(v), we construct the normalized cumulative velocity distribution |$F(v) = \int _{-\infty }^vf(v^{\prime }){\rm d}v^{\prime }/\int _{-\infty }^{+\infty }f(v^{\prime }){\rm d}v^{\prime }$|⁠. Since the velocity profile is a noiseless non-negative function, F(v) is strictly monotonically increasing. We then determine the velocities at which 5, 25, 50, 75 and 95 per cent of the line flux accumulates. The median velocity v₅₀ is the solution of the equation F(v) = 0.5. The width comprising 90 per cent of the flux is w₉₀ = v₉₅ − v₀₅, the width at 80 per cent is w₈₀ = v₉₀ − v₁₀ and the width comprising 50 per cent of the flux is w₅₀ = v₇₅ − v₂₅. All these values have dimensions of velocity (km s⁻¹). For a Gaussian profile, the value w₈₀ is close to the conventionally used full width at half-maximum (FWHM; w₈₀ = 2.563σ = 1.088 FWHM; w₉₀ = 3.290σ). For a typical object in our sample (median w₈₀ = 752 km s⁻¹) the instrumental dispersion of the SDSS (σ_inst = 70 km s⁻¹) contributes only a few per cent to the line width (Liu et al. 2013b), and as the line profiles are typically non-Gaussian, we do not attempt to deconvolve the resolution except a couple of cases noted explicitly when we use Gaussian components individually.

We can measure the asymmetry of the velocity profile relative to the median velocity by computing a dimensionless relative asymmetry R = ((v₉₅ − v₅₀) − (v₅₀ − v₀₅))/(v₉₅ − v₀₅). Negative values correspond to cases where the blueshifted wing of the line extends to higher velocities than the redshifted one, and positive values correspond to cases where the redshifted wing dominates. This measure is a non-parametric analogue of skewness and is equal to 0 for any symmetric profile (including a single Gaussian). Furthermore, we can measure the prominence of the wings of the profile, or the non-parametric analogue of the kurtosis, by computing r₉₀₅₀ ≡ w₉₀/w₅₀. For a Gaussian profile, this value is equal to 2.4389. Values higher than this indicate profiles with relatively more extended wings than a Gaussian function: for example, a Lorentzian profile f(v) = 1/(γ² + v²) (where γ is the measure of the profile width) has r₉₀₅₀ = 6.3138. Values lower than the Gaussian value indicate a profile with a stronger peak-to-wings ratio and are rarely encountered in our sample.

Finally, we compute the absolute asymmetry of the profile, which is A = (flux(v > 0)−flux(v < 0))/total flux. In terms of the normalized cumulative velocity distribution, A = 1 − 2F(0). This asymmetry is dimensionless and it is positive for profiles with more flux at redshifted wavelengths than at blueshifted wavelengths.

Values A and v₅₀ critically depend on an accurate determination of the host galaxy redshift, because this is what we use to fix the v = 0 point. If no absorption features in the composite stellar light of the host are detected, the redshift can only be determined from the emission lines themselves, which renders the absolute velocity and skewness meaningless. Values of R, w₉₀, w₈₀, w₅₀, and r₉₀₅₀ include only differences between velocities and do not hinge on the accurate determination of the host velocity.

2.3 Robustness of non-parametric measures

In this section, we evaluate the performance of the non-parametric measures. The theoretical advantage of the non-parametric measures is in their relative insensitivity to the fitting functions used. We test this assumption by repeating all the fits using a set of one, two or three Lorentzian (f(v) = 1/(γ² + v²)) profiles. The Lorentzian function has significantly more flux in the faint wings than does the Gaussian function, and this shape is not borne out in the observations of the line profiles in our sample. The quality of the fits with Lorentzian profiles is significantly poorer than of those with Gaussian profiles, and therefore all our final non-parametric measures are based on multi-Gaussian fitting as described in the previous subsection.

Nevertheless, we carry out the comparison between the non-parametric measures derived from the two methods. We find that both sets of fits yield nearly identical absolute asymmetries and median velocities v₅₀, which is not surprising because both these measures are most sensitive to the correct identification of the line centroid. All other measures (the widths, relative asymmetry and r₉₀₅₀) are strongly correlated between the two sets of fits, but the specific values are systematically different. The line width w₈₀ as measured from the Lorentzian fits is about 25 per cent higher than that from the Gaussian fits; the relative asymmetry is significantly weaker as measured by Lorentzian profiles than the one measured by the Gaussian ones; and r₉₀₅₀(Lorentzian) is approximately equal to r₉₀₅₀(Gaussian)+2. All these differences are as expected from the fitting functions with different amounts of power in the extended wings.

The conclusions we derive from this comparison are two fold. First, since the multi-Lorentzian fits are not only statistically but also visibly inferior to the multi-Gaussian ones, the real systematic uncertainty on the w_80,90 – the key measurements discussed in this paper – is significantly smaller than the 25 per cent difference between line widths calculated from these two methods. This is very encouraging. (For the majority of objects, w₈₀ is accurate to 10 per cent or better, as measured from the comparison of non-parametric measures derived from two-Gaussian and three-Gaussian fits.) Secondly, we confirm that the non-parametric measures are relatively robust: although the Lorentzian profiles do not yield statistically good fits, they nevertheless give reasonable estimates of the non-parametric measures.

We perform an additional test to determine the effect of the S/N of the spectra on our measurements. A narrow Gaussian emission line with a weak broad base observed with a high S/N is represented by two Gaussians in our multi-Gaussian fit, and its non-parametric measures include the power contributed by the weak broad base. On the contrary, if the same object is observed in a lower quality observation, the weak base is not necessarily recognized as such because the χ² of the two-Gaussian fit may be indistinguishable from the one-Gaussian one and thus the latter will be preferred.

We conduct the following Monte Carlo test to explore the effect of noise on our measurements. We take eight of the highest S/N objects in our catalogue, four with w₈₀ > 1000 km s⁻¹ and four with w₈₀ < 500 km s⁻¹. We then downgrade the quality of these spectra by adding progressively higher Gaussian random noise to the original (essentially noiseless) spectra and conduct all our multi-Gaussian and non-parametric measures in the manner identical to that used for real science observations. The results are shown in Fig. 3. Both absolute asymmetry and v₅₀ are relatively insensitive to the noise and do not show systematic trends. The line width w₈₀ does not depend on the noise for the four objects with relatively narrow (w₈₀ < 500 km s⁻¹) lines. But for the broad-line objects the measured w₈₀ noticeably declines as the quality of observations worsens, reflecting the ‘missing broad base’ phenomenon. Measurements with peak S/N > 10 are relatively safe from this phenomenon: only one of the eight objects shows a noticeable decline of w₈₀ at S/N ≲ 30. The relative asymmetry and the kurtosis-like r₉₀₅₀ quickly drop to single-Gaussian values as the S/N decreases below 20 or so.

Figure 3.

The dependence of non-parametric measures on the peak S/N of the [O iii] emission line for eight spectra whose quality is progressively artificially downgraded by adding Gaussian random noise.

Open in new tabDownload slide

3 KINEMATIC ANALYSIS OF INTEGRATED SPECTRA

3.1 Outflow signatures

Interpreting the line-of-sight gas kinematic measurements in terms of physical 3D motions of the gas is highly non-trivial. This is true even when spatial information is available, for example, via integral-field observations (Liu et al. 2013b), but is even more so when we have to rely only on one spatially integrated spectrum. The reason is that if an object exhibits a spherically symmetric optically thin outflow, then its emission-line profiles are symmetric and peaked at zero velocity, since there is a large amount of gas moving close to the plane of the sky. Therefore, there is simply no ‘smoking gun’ outflow signature in the emission-line profile of such a source. If an outflow has non-zero optical depth, then its presence can be inferred by its absorption of the background light at wavelengths blueshifted relative to the quasar rest frame, as happens in quasar absorption-line systems (Crenshaw, Kraemer & George 2003; Arav et al. 2008), but this is not our case. Therefore, proving that a given line profile is due to gaseous outflow is often difficult and relies on indirect arguments (Liu et al. 2013b).

To investigate the relationship between the observed line-of-sight velocity dispersion of the lines and the typical outflow velocity, we consider for the moment a spherically symmetric outflow with a constant radial velocity v₀. Because different streamlines have different inclinations to the line of sight, the observer sees a range of velocities – in this simplest case of the spatially integrated spectrum, the line profile is a top-hat between −v₀ and v₀, so that w₈₀ = 1.6v₀. In Liu et al. (2013b), we consider an outflow with a power-law luminosity density and calculate the velocity profiles in a spatially resolved observation, finding a typical w₈₀ ≃ 1.3v₀ in the outer parts of the narrow-line region. If a flow consists of clouds moving on average radially with velocity v₀, but also having an isotropic velocity dispersion σ, then the observed w₈₀ can be approximately calculated as |$w_{80}\simeq \sqrt{(1.6v_0)^2+(2.42\sigma )^2}$|⁠, as long as v₀ and σ are not too dissimilar.

Another simple case is a radial flow, in which at every point there are clouds with a range of radial velocities. As an example, we consider the results of the 2D simulations of quasar feedback by Novak, Ostriker & Ciotti (2011), in which at every distance from the quasar the higher density regions (‘clouds’) have a wide distribution of radial velocities (Novak, private communication) which ranges from 0 to ≳ 1000 km s⁻¹, with a median among all clouds of v₀ = 220 km s⁻¹. We use the velocities of clouds straight from these simulations and several different luminosity density profiles to produce mock emission-line profiles.

The resulting profiles are shown in Fig. 4; they are insensitive to the adopted luminosity density profile because the velocity distribution of clouds does not change appreciably as a function of distance in these simulations. The profiles are peaked at zero velocity because there is a large population of clouds with small radial velocities, while the velocity dispersion of clouds is neglected in our calculations. If the clouds have an isotropic velocity dispersion which does not vary with distance, a more faithful profile can be obtained by convolving the profiles in Fig. 4 with a Gaussian, which would make the profiles broader and less peaky. The measured velocity width of the profile is w₈₀ ≃ 1.4v₀, similar to the scaling obtained in other simple cases, despite the broadness of the velocity distribution in this example.

Figure 4.

Mock emission-line velocity profiles constructed from the cloud velocity distribution from Novak et al. (2011), using six different emissivity profiles (from linearly declining, to Gaussian, to flat, to centrally tapered power law). All six curves are essentially on top of one another because the velocity distribution of the clouds hardly varies with the distance from the quasar in these simulations. The actual median radial velocity of the clouds is 220 km s⁻¹, and we measure w₈₀ = 310 km s⁻¹ and w₉₀ = 470 km s⁻¹ for the simulated profile.

Open in new tabDownload slide

To sum up, (i) there is no tell-tale outflow signature in an optically thin, spherically symmetric outflow; (ii) deviations from spherical symmetry (and moreover from axial symmetry) are required to produce asymmetric line profiles; and (iii) the velocity width of the emission line can be used to estimate the outflow velocity, w₈₀ ≃ (1.4–1.6) × v₀. The most natural way in which the symmetry may be expected to be broken is due to dust obscuration, either by dust embedded in the outflow itself or by dust concentrated in the galactic disc. In either case, the redshifted part is more affected by extinction, and thus excess blueshifted emission is considered a sufficient indicator of an outflow (Heckman et al. 1981; De Robertis & Osterbrock 1984; Whittle 1985a; Wilson & Heckman 1985). In such case, the apparent w₈₀ decreases typically by ≲ 30 per cent if the extinction is ≲ 2.5 mag and concentrated in a disc (Liu et al. 2013b), reducing w₈₀/v₀ by the same amount. Thus, for a given w₈₀, an asymmetric profile indicates a somewhat higher v₀ than a symmetric one.

Double-peaked profiles are expected in some geometries for a bi-conical outflow or more complex outflow kinematics (Crenshaw & Kraemer 2000), but can also be due to the rotation of the galaxy disc or two (or more) active nuclei in a merging system of galaxies, each illuminating its own narrow-line region. Distinguishing these possibilities usually requires follow-up observations at high spatial resolution, and the relative frequency of these scenarios remains a matter of debate (Comerford et al. 2009; Liu et al. 2010; Shen et al. 2010; Fu et al. 2012; Barrows et al. 2013; Blecha, Loeb & Narayan 2013), but it appears that outflows dominate over dual active nuclei. It is likely that complex outflow kinematics is responsible for the majority of split-line profiles in our sample, and we conduct non-parametric kinematic measurements of such objects in the same way we do for the rest of the sample and include them in all our analyses.

3.2 Analysis of [O iii] kinematics

In Fig. 5 we present the results of the non-parametric measurements of the [O iii] line in our sample of type 2 quasars. Both the relative asymmetry R and the absolute asymmetry A demonstrate slight preference for negative values, i.e. for blue excess. The blue asymmetries indicate that there is at least some outflow component in the [O iii]-emitting gas in type 2 quasars. The sample means and standard deviations for these values are A = −0.03 ± 0.15 and R = −0.08 ± 0.15. Line widths, relative asymmetries, kurtosis parameters and (to a lesser extent) median velocities are correlated with one another, in the sense that objects with broader lines also have a more pronounced blue excess (negative R), higher r₉₀₅₀ and more negative v₅₀ (for the latter, the sample mean and standard deviation are −14 ± 75 km s⁻¹). Thus, as we discussed in the previous section, it is plausible that the [O iii] velocity width can serve as a proxy for the outflow velocity. Similar relationships have been reported by other authors in type 1 quasars, e.g. by Steinhardt & Silverman (2013), who find a stronger blueshift of [O iii] relative to broad Hβ and [O ii] λλ3726, 3729 as a function of [O iii] width.

Figure 5.

Distributions of non-parametric kinematic measures and their mutual relationships. Points are colour-coded from light grey to dark grey by the peak S/N (lightest have S/N ≲ 20, grey have S/N from 20 to 50 and black have peak S/N ≳ 50). In all the plots involving absolute asymmetry A and median velocity v₅₀ only points with reliable host redshifts are displayed since both these values critically depend on the host redshift determination.

Open in new tabDownload slide

From these correlations, a picture emerges in which the outflow component, or at least the component of the outflow which is more affected by the obscuration, tends to be broad. We further explore this notion in Fig. 6, where we split the line profiles into a ‘broad’ and a ‘narrow’ component. For objects with two-Gaussian fits, the designation is straightforward, but the majority of objects require three Gaussians, in which case we pick the two most luminous ones and designate them ‘broad’ and ‘narrow’ according to their velocity dispersions. In Fig. 6, we show that indeed the broader of the Gaussian components is the one that tends to be blueshifted relative to the narrow ones. The narrow cores tend to be well centred in the host galaxy frame; the mean and standard deviation of the centroids of the narrower Gaussian components is 3 ± 150 km s⁻¹ for the 271 objects with accurately determined host redshifts. On the contrary, the broader components tend to be slightly blueshifted, both relative to the host galaxy frame (velocity centroid of −60 ± 210 km s⁻¹ in the 271 objects with accurate host redshifts) and relative to the narrow components (v_c,broad − v_c,narrow = −90 ± 270 km s⁻¹ for the entire sample of 568 sources).

Figure 6.

Velocity offset between the centroids of the broad component from the narrow component (black for profiles decomposed into two Gaussians, grey for the two Gaussians that dominate the flux in three-Gaussian decompositions). v_c is for the velocity centroids. Broad components tend to be blueshifted relative to the narrow ones. The velocity dispersions of individual Gaussian components have been corrected for the instrumental resolution (σ_inst = 70 km s⁻¹ subtracted in quadrature).

Open in new tabDownload slide

Because the narrow component is well centred in the host frame, it is tempting to postulate that the narrow Gaussian component tends to be produced by gas in dynamical equilibrium with the host galaxy, e.g. in rotation in the galaxy disc, and is simply illuminated by the quasar, whereas the broad component is due to the outflow. However, we hesitate to make this inference, as there is no particular reason to assign any physical meaning to the individual parameters of the Gaussian components. We again draw a lesson here from Fig. 4, where the mock emission-line profile is due entirely to the outflowing clouds and can be decomposed into several Gaussian components, none of which correspond to the gas in rotation in the host galaxy.

We go further in Fig. 7, where we show that neither the overall line width nor the width of the narrower Gaussian component show any correlation with the stellar velocity dispersion. Since even among the 271 objects where the host galaxy was detected many of the stellar velocity dispersions are rather poorly determined for the reasons discussed in Section 2.1, for this figure we use only the better determined stellar velocity dispersions from Greene et al. (2009). That the overall line width (left-hand panel) shows no relationship with stellar velocity dispersion is not surprising if most of the line width arises due to the outflow. But the narrow cores do not appear to show any relationship with the stellar velocity dispersion either. As a result, it seems likely to us that virtually none of the [O iii]-emitting gas in type 2 quasars is in dynamical equilibrium with the host galaxy. This is in contrast to the situation in lower luminosity active galaxies in which [O iii] width strongly correlates with galaxy rotation and/or bulge velocity dispersion (Wilson & Heckman 1985; Whittle 1992; Nelson & Whittle 1996; Greene & Ho 2005). In such objects, it appears that the gas motions are in accord with the gravitational forces in the galaxy, and the gas is simply illuminated and photoionized by the active nucleus to produce the narrow-line region.

Figure 7.

Taking just the objects with well-determined stellar velocity dispersions from Greene et al. (2009), we plot the overall width of the [O iii] emission line in the left-hand panel and the dispersion of just the narrow component in the right-hand panel as a function of the stellar velocity dispersion (grey points for the narrower of the two dominant components in a three-Gaussian fit, black points for the narrower of the two components in a two-Gaussian fit). The left-hand panel is similar to fig. 7 of Greene et al. (2009) even though the exact non-parametric measures of dispersion used in that paper were defined differently. The error bars on w₈₀ reflect non-Gaussianity of lines and are calculated by converting w₈₀ from w₅₀ and w₉₀ assuming a Gaussian profile; thus, for a Gaussian profile they would be zero. Stellar velocity dispersions and σ_narrow have been corrected for instrumental dispersion; since lines are non-Gaussian, w₈₀ values are plotted as observed (typical correction is a few per cent).

Open in new tabDownload slide

The relative shifts between the narrower and the broader components could arise if the narrow component is produced on all scales in the host galaxy, where it is less likely to suffer from strong extinction, whereas the broader component is produced closer to the nucleus, where it is more likely to be affected by extinction. This is consistent with an outflow that is driven close to the nucleus and then gradually slowed down by the interactions with the interstellar medium (Wagner, Umemura & Bicknell 2013). Furthermore, this picture is consistent with the apparent decline of the line width in the outer parts of the outflow seen in the integral-field unit observations of type 2 quasars (Liu et al. 2013b), although the effect is small, with velocity width declining only by 3 per cent per projected kpc.

The sample mean and standard deviation of line width is w₉₀ = 1230 ± 590 km s⁻¹, with a median of 1060 km s⁻¹, minimum of 370 km s⁻¹ and maximum of 4780 km s⁻¹ (w₈₀ = 880 ± 430 km s⁻¹, median 752 km s⁻¹, min 280 km s⁻¹, max 2918 km s⁻¹), much higher than that of local ultraluminous infrared galaxies (ULIRGs; median w₉₀ ≃ 800 km s⁻¹) and especially of those without a powerful active nucleus in their centre (median w₉₀ ≃ 600 km s⁻¹ for pure starbursts; Hill & Zakamska 2014). For ULIRGs, the line widths and the outflow velocities strongly correlate with the power source (higher for active galaxies, lower for starbursts; Rupke & Veilleux 2013a; Hill & Zakamska 2014), and the majority of the objects in our sample show line widths consistent with quasar-driven outflows, as expected. In Liu et al. (2013b), we estimated that for gas discs rotating in the potential of the most massive galaxies line widths do not exceed w₈₀ ≃ 600 km s⁻¹. Thus, the line-of-sight gas velocities that we see in our sample are too high to be confined even by the most massive galaxy potential, and this gas cannot be in dynamical equilibrium with the host galaxy.

The range of [O iii] luminosities in our sample is not all that large, with 90 per cent of sources between |$\log (L{\rm [O\,\small {III}]}/\rm{L}_{{\odot }}) = 8.5$| and 9.5, and the remaining 10 per cent sources spanning the higher decade in luminosity. There is some tendency of objects with more pronounced outflow signatures (higher width, higher r₉₀₅₀) to have higher [O iii] luminosity (Spearman rank correlation coefficient r_S ≃ 0.19 for both relationships, probability of the null hypothesis of uncorrelated data sets P_NH = 10⁻⁵), however, no correlations are seen between L[O iii] and absolute asymmetry, relative asymmetry or median velocity (Table 1). Furthermore, any correlations between L[O iii] and kinematic measures can be strongly affected by the cutoff in the [O iii] luminosity distribution which is due to our sample selection (⁠|$\log (L{\rm [O\,\small {III}]}/\rm{L}_{{\odot }})\ge 8.5$|⁠. We report the significance of correlations at face value, without trying to account for the effects of the luminosity cutoff.

3.3 Fainter lines

We compare the kinematic structure of the brightest lines, [O iii] λ5007, [O ii] λλ3726, 3729 and Hβ. It has long been known that line kinematics often vary as a function of the ion ionization potential or the line critical density (Whittle 1985b). On the basis of largely anecdotal evidence, we previously established that in type 2 quasars with highly asymmetric or split-line [O iii] λ5007 profiles the [O ii] λλ3726, 3729 profiles seemed less complex (Zakamska et al. 2003). The difficulty of this analysis is illustrated in Fig. 3: because [O iii] λ5007 is by far the brightest line, it is much easier to miss a weak broad component in [O ii] than in [O iii], giving the impression that the [O ii] is narrower or lacks kinematic structures present in [O iii]. To remedy this problem, in what follows we use only the emission lines detected with peak S/N > 10 (about 230 objects), and we prefer S/N > 20 (about 140 objects).

In Fig. 8, we show the comparison between [O iii], [O ii] and Hβ line widths. Hβ is slightly (8 per cent) systematically narrower than [O iii] on average. Much of the difference is likely attributable to the S/N effect described in Section 2.3, since it is easier to miss a weak broad component in a noisy profile of Hβ than in the much higher S/N profile of [O iii]. There are only a few cases where the ‘by eye’ examination of the [O iii] and Hβ profiles superposed on one another reveals that the [O iii] is genuinely significantly broader than the noisier Hβ (Fig. 9).

Figure 8.

Comparison between non-parametric width measurements for [O iii] λ5007, Hβ and [O ii] λλ3726, 3729. In light grey are sources with estimated peak S/N of Hβ (left) and [O ii] (right) between 10 and 20; in black are sources with S/N > 20. Because of the minimal S/N requirement, 223 sources appear in the left-hand panel and 237 in the right. For [O ii], the top of each bar corresponds to the non-parametric measure of the width of the entire non-deblended doublet, whereas the bottom of the bar includes the correction for the doublet splitting according to equation (1).

Open in new tabDownload slide

Figure 9.

[O iii] (dotted grey), Hβ (solid grey) and [O ii] (black) line profiles for the 10 objects with the broadest [O iii] which also have S/N([O ii]) > 10. The [O ii] profiles are as observed, without any accounting for the doublet nature of the profile; at these extreme line width the doublet splitting ∼220 km s⁻¹ is not noticeable. Both Hβ and [O ii] are capable of showing asymmetries and complicated profile features. The kinematics of Hβ tends to be less extreme than that of [O iii], and of [O ii] even less so, with Hβ typically seen in between [O iii] and [O ii] profiles.

Open in new tabDownload slide

In contrast, [O ii] is noticeably narrower than [O iii] for high-width objects (Fig. 8). The average doublet-corrected width of [O ii] over all 237 objects with peak S/N[O ii]>10 is only 5 per cent smaller than that of [O iii]; however, when we consider only objects with w₉₀[O iii]>1500 km s⁻¹, the difference in width increases to 27 per cent. Overall while [O ii] sometimes displays asymmetries and complicated profiles, it does not show the extremely broad features seen in [O iii], with Hβ demonstrating kinematics that are intermediate between [O ii] and [O iii].

We conduct a similar examination of He ii λ4686 and [O iii] λ4363 profiles, but only a handful of objects with S/N > 10 in these lines have relatively broad [O iii]. In these objects, [O iii] λ5007, He ii λ4686 and [O iii] λ4363 kinematic structures look, within the uncertainties, consistent with one another.

Line flux ratios involving [O iii] λ5007, Hβ, [O ii] λλ3726, 3729, He ii λ4686 and [O iii] λ4363 (Osterbrock & Ferland 2006; Liu et al. 2013a) are diagnostic of the ionization state and temperature of the gas. We look for trends in these flux ratios as a function of all kinematic parameters and luminosity of [O iii]. Unfortunately, the variations in these line fluxes are subtle enough that high S/N values in the fainter lines are required to measure their fluxes to the required accuracy. In particular, if it is the weak broad components that vary as a function of kinematics and/or luminosity, an S/N of several tens in these lines would be required to detect these trends, as shown in Fig. 3. As a result, we do not find any definitive trends in any of these ratios even when restricting the analysis to the small number of the highest S/N objects. Instead, we perform such measurements in Section 5 using composite spectra.

4 KINEMATICS AND MULTIWAVELENGTH PROPERTIES

4.1 Kinematic indicators and radio emission

The exact shape of the radio luminosity function of active nuclei remains a topic of much debate, including whether there is a very broad distribution of intrinsic luminosities or whether there is a true dichotomy in this property (Kellermann et al. 1989; Xu, Livio & Baum 1999; Ivezić et al. 2002; Jiang et al. 2007; Condon et al. 2013). In any case, the objects on the high end are traditionally called ‘radio loud’ and are sometimes very extended in the radio band. In these cases, collimated relativistic jets propagate out to several hundred kpc from the host galaxy, and the radio emission of these sources is dominated by the lobes where the jet energy finally dissipates. The majority of our matches are weak (a few mJy) point sources at the 5 arcsec resolution of the FIRST survey. Since the optical broad-band emission of type 2 quasars is a poor proxy for their luminosity, we use the distribution of our sources in the [O iii] luminosity/radio luminosity plane to define the radio-loud/radio-quiet boundary (Xu et al. 1999) and find that about 10 per cent qualify as classical radio-loud sources, and a similar fraction are significantly ( ≳ 10 kpc) extended (Zakamska et al. 2004). For Seyfert galaxies, Ho & Peng (2001) find that some radio-loud candidates can be missed when using the spatially integrated optical or infrared luminosity in defining the radio-loud/radio-quiet boundary and suggest using the nuclear properties in such definitions. Luckily in our case the quasars are extremely luminous and dominate the bolometric output, so our definition is unlikely to be affected by this bias.

We find a strong correlation between the [O iii] line width and the radio luminosity (Fig. 10; Table 1). For nearby lower luminosity active galaxies, similar relationships were previously reported by many authors (Heckman et al. 1981; Wilson & Heckman 1985; Whittle 1992; Nelson & Whittle 1996) and by Veilleux (1991b) whose sample is shown in Fig. 10 for comparison. While there is a small tail of objects with very high radio luminosities, most of the correlation is due to the cloud of points at νL_ν[1.4 GHz] = 10³⁹-10⁴¹ erg s⁻¹. Although these radio luminosities seem high by comparison to those of the local active galaxies (e.g. the red points from the Veilleux 1991b sample), we need to keep in mind that the [O iii] luminosities of our type 2 quasars are at the extreme end of the luminosity distribution (right-hand panel of Fig. 10). Thus, in the space of [O iii] versus radio luminosities (Xu et al. 1999) most of the type 2 quasars in our sample follow the radio-quiet, rather than the radio-loud, locus. Recently, Mullaney et al. (2013) also found a trend of increasing line width with radio luminosity using composite spectra of type 1 quasars with a wide range of [O iii] luminosities which overlaps with ours. Similarly, Spoon & Holt (2009) reported a correlation between outflow signatures in mid-infrared emission lines and radio luminosities among ULIRGs. Because these objects tend to be more dust obscured on galaxy-wide scales and have high rates of star formation which contributes to their radio emission, a direct comparison between our sample and theirs is complicated, but at face value the radio luminosities of the two samples do overlap.

Figure 10.

Left: radio luminosities (k-corrected to the rest-frame 1.4 GHz) versus the velocity width w₉₀ of the [O iii] emission. Crosses and upper limits show 568 type 2 quasars in our kinematic analysis. The 568 type 2 quasars in our kinematic analysis are colour-coded by robustness of line fits (light grey, light blue for peak S/N([O iii]) < 20, black, dark blue for higher S/N), with crosses for radio detections and arrows for radio upper limits. Orange circles show Seyfert galaxies from Veilleux (1991b); NVSS is used to determine the radio fluxes of these sources. The thick solid red lines show median radio luminosities in bins of line width. The top line assumes that radio non-detections are close to the limit of the survey, the bottom line assumes that non-detections are 10 times fainter than the limit of the survey, and the middle line assumes that the non-detections are two times fainter than the limit of the survey. The dotted line shows νL_ν[1.4 GHz] |${\propto}\;w_{90}^2$|⁠; it is not a fit to the data. Right: there is a slight tendency for objects with broader [O iii] to also have higher [O iii] luminosity, but because of the cutoff of the [O iii] luminosity the significance of this relationship is difficult to establish.

Open in new tabDownload slide

Taking all sources and radio upper limits at face value (i.e. assuming that the non-detections are close to the survey limit) yields a Spearman rank correlation coefficient of r_S = 0.29. Excluding radio-loud sources using their distribution in the [O iii]–radio plane (Zakamska et al. 2004) yields r_S = 0.33. If the non-detections are significantly below the survey limit, then the correlation is even stronger because there are hardly any non-detections on the high-w₉₀ end of the diagram. In Section 6.1, we demonstrate that the radio fluxes of objects not detected by FIRST are likely within a factor of 2 of the FIRST survey limit. If we suppress all upper limits by a factor of 2 (but keep radio-loud sources in), the correlation has r_S = 0.34.

We also find an anticorrelation between radio luminosity and absolute asymmetry, r_S = −0.37, in the sense that objects with stronger blue asymmetry tend to have stronger radio emission. A similar relationship exists between radio luminosity and the median velocity v₅₀. In all these cases the null hypothesis (that the two data sets are uncorrelated) is rejected with P_NH < 10⁻⁵. Radio luminosity is not correlated with the other kinematic measures of outflow activity (relative asymmetry and kurtosis r₉₀₅₀). For comparison, we also show the weaker trend between w₉₀ and [O iii] luminosity reported in Section 3.2 in the right-hand panel of Fig. 10.

4.2 Kinematic indicators and infrared luminosity

We cross-correlate the entire sample of SDSS type 2 quasars (Reyes et al. 2008) against the Wide-field Infrared Survey Explorer (WISE) catalogue within 6 arcsec. Out of 887 objects in the catalogue, 876 objects have matches in W1 (3.6 μm) and W2 (4.5 μm); 829 objects in W3 (12 μm); and 773 objects in W4 (22 μm) with S/N above 2.5. The 11 objects without W1 matches are visually examined; in almost all cases there is an actual detection at the position of the quasar, but it is blended with a brighter nearby object and is thus not reported in the catalogue. We interpolate between the WISE fluxes using piece-wise power laws to calculate νL_ν at rest-frame 5 and 12 μm and the index between these two, νL_ν ∝ λ^β (higher index means redder spectral energy distribution). As the majority of the sources are well above the detection limit for the survey, the analysis of the WISE matches is not affected by non-detections to the same extent as the analysis of radio emission in the previous section.

The mid-infrared luminosities of type 2 quasars in our sample strongly correlate with their radio luminosities, [O iii] velocity widths and [O iii] luminosities (Fig. 11). The relationships between infrared luminosities, radio luminosities and [O iii] luminosities in low-luminosity active galaxies have been pointed out by many authors before, most recently by Rosario et al. (2013). For direct comparison with their work, we show the locus of normal star-forming galaxies and the Rosario et al. (2013) line separating two branches of active nuclei in the left-hand panel of Fig. 11. Using νL_ν at 5 μm instead results in similar relationships albeit with somewhat larger scatter.

Figure 11.

WISE monochromatic k-corrected luminosities as a function of radio luminosities (left), [O iii] velocity widths (middle) and [O iii] luminosities (right). In the left-hand panel, the blue points are those with upper limits on radio fluxes, whereas the black points are radio detections. The red ellipse is the locus of star-forming galaxies from Rosario et al. (2013), and the dashed line is their separation line between infrared bright (above the line) and infrared faint (below the line) branches of active nuclei. In the middle panel, grey points are for objects with peak S/N[O iii] < 20.

Open in new tabDownload slide

Mid-infrared luminosities νL_ν[12 μm] are correlated with all kinematic measures, in the sense of higher mid-infrared luminosity in objects with stronger outflow signatures, with significance ranging between P_NH = 0.02 and <10⁻⁵ (Table 1). The strongest correlation is with w₉₀ (r_s = 0.44, P_NH < 10⁻⁵). There is a hint (P_NH = 0.04) that the mid-infrared luminosity correlates positively with positive values of v₅₀, while correlating negatively with the negative values, suggesting that either a strong blue asymmetry or a strong red asymmetry may be a sign of an outflow. Objects with higher w₉₀ tend to have redder mid-infrared spectral energy distributions (higher β), with r_S = 0.23 and P_NH < 10⁻⁵.

Rosario et al. (2013) point out that the Seyfert galaxies in their sample lie almost exactly on top of the locus of the star-forming galaxies in the radio/infrared diagrams. These authors conclude that only 15 per cent of the infrared flux of these objects is due to the active nucleus and that the correlation between infrared and radio fluxes seen among Seyfert galaxies is simply a reflection of the standard radio/infrared correlation due to star formation. Type 2 quasars appear to lie on the luminous extension of the locus of the star-forming galaxies, and thus it is tempting to postulate that the same arguments apply in our sample, except the star formation rates of the host galaxies must be much higher than those seen in Seyferts by Rosario et al. (2013).

However, this explanation is unlikely to extend to the objects in our sample. The mid-infrared colours and fluxes of type 2 quasars at these luminosities are dominated by the quasar, not by the host galaxy (Lacy et al. 2004; Stern et al. 2005; Zakamska et al. 2008). Thus, the strong correlation between radio and mid-infrared in this regime (and the excess of the radio emission over the amount seen in nearby star-forming galaxies) suggests that the radio emission in radio-quiet quasars is related to the quasar activity, not to the star formation in its host.

Our sample has 54 objects in common with Jia et al. (2013) who analysed XMM–Newton and Chandra snapshots of a large sample of obscured quasars deriving their X-ray luminosities, spectral slopes and amount of intervening neutral gas absorption. These objects were either targeted by X-ray observatories or serendipitously lie in the fields of view of other targets. Similarly to Veilleux (1991b), we do not find any correlations between any of the kinematic indicators and any of the X-ray spectral fitting parameters. In particular, there is no correlation between the optical line width and the absorption-corrected (intrinsic) X-ray luminosity. Removing the 24 Compton-thick candidates (in which obtaining intrinsic X-ray luminosities is particularly difficult) still reveals no relationship between X-ray parameters and gas kinematics.

One possibility is that the lack of such correlation implies the lack of strong influence of X-ray emission on the launching of the winds. Another possibility (which we find more likely) is that the existing X-ray observations of obscured quasars are not yet of sufficient quality to probe this relationship. The uncertainties in the intrinsic X-ray luminosities are rather high, because the observed fluxes need to be corrected for intervening absorption. As a result, even the correlation between X-ray and mid-infrared luminosities – which are both supposed to be tracers of the bolometric luminosity – is rather weak (P_NH ≃ 0.02). Among the 30 Compton-thin sources, the X-ray to mid-infrared luminosity ratio has a dispersion of 0.7 dex, similar to the value reported for local Seyfert 2 nuclei (LaMassa et al. 2011) who argue that correcting for intrinsic absorption is difficult even when high-quality X-ray observations are available.

We perform a simulation in which we randomly draw 30 points from the infrared versus kinematics correlation in Fig. 11, middle. We find that the correlation is still detected with P_NH ≲ 0.01 significance. But if we add a Gaussian random variable with a 0.7 dex dispersion to the log of the infrared luminosity, the correlation is no longer detected. Thus either an intrinsic dispersion or observational uncertainties (related to difficulties of correcting for intervening absorption) of this magnitude are sufficient to destroy a correlation in a sample of 30 objects (the number of Compton-thin obscured quasars with available X-ray luminosities). It will be interesting to probe the connections between the ionized gas kinematics and ultraviolet, optical and X-ray luminosities in type 1 quasars, where correcting for intervening absorption is not a significant problem and where the relative strengths of these correlations could elucidate the primary driving mechanism of the ionized gas outflows.

5 COMPOSITE SPECTRA

5.1 Constructing composites

To further test the trends we find in Sections 3 and 4 and to study weak emission features, we produce several sets of composite spectra. We choose a quantity that is easily measurable in every object (e.g. [O iii] line width in our first example) and bin the sample into five equal-size bins (114 objects) in this quantity. We then arithmetically average all host-subtracted spectra within each bin. This allows us to obtain high S/N composites while being able to tease out the dependences on the chosen parameter. The composites produced in five bins of the [O iii] velocity width are presented in Fig. 12 and the composites in five bins of [O iii] luminosity in Fig. 13. Below we also discuss composites made in bins of infrared and radio luminosity, although these are not shown.

Figure 12.

Composite in five bins of [O iii] w₈₀ width (from light to dark w₈₀ is increasing) normalized to the peak flux density of the [O iii] flux. Composite is straight-up error weighted mean in the host frame (whether or not the host subtraction is reported to be accurate). The bins are w₈₀ = 284-546 km s⁻¹, 546-673 km s⁻¹, 673-856 km s⁻¹, 856-1153 km s⁻¹ and 1153-2918 km s⁻¹. Most of the line peaks line up from one composite to the next, but the peaks of [S ii] λλ4069, 4076 and [N i] λλ5198, 5200 show an increase. [FeV ii] λ5158 appears to increase as well.

Open in new tabDownload slide

Figure 13.

Composite in five bins of [O iii] luminosity (from light to dark L[O iii] is increasing) normalized to the peak flux density of the [O iii] flux. The bins are log L[O iii]/L_⊙ = 8.50-8.61, 8.61-8.78, 8.78-8.97, 8.97-9.23, 9.23-10.41. Many of the line peaks decrease in the red→ blue sequence, but the line widths do not vary drastically, reflecting the weakness of the relationship between [O iii] luminosities and widths.

Open in new tabDownload slide

Even though the spectra are already host subtracted, for accurate measurements of very weak lines we need a more accurate continuum subtraction. We select a dozen continuum-dominated wavelength intervals and spline-interpolate between them to produce the model of the continuum which is then subtracted from the composite spectrum. Continuum subtraction is the dominant source of systematic uncertainty in measuring the weakest lines.

We then construct the velocity profile of the narrow-line region from the [O iii] λλ4959, 5007 lines and we use this profile to fit about 15 emission lines. Several doublets which have wavelength separations that are too small to be resolved in our spectra are fitted with fixed ratios between the components, as follows: [O ii] λλ3726, 3729 with a 1:1 ratio, [S ii] λλ4069, 4076 with a 3:1 ratio (the doublet structure of this line is clearly visible in the lower w₈₀ composites) and [N i] λλ5198, 5200 with a 1:1 ratio. For each emission feature, given the velocity profile from the [O iii] line, there is only one adjustable parameter – its amplitude. The fit is linear in all 15 amplitudes.

A few features are close blends: He i λ3889 is blended with Hζ (8→2 transition) and [Ne iii] λ3968 is blended with Hε (7→2 transition). In both cases, the non-hydrogen emission makes a larger contribution to the blend, but not an overwhelmingly dominant one. To measure He i λ3889 and [Ne iii] λ3968, we estimate the Balmer decrement from Hβ, Hγ and Hδ and use the derived extinction values to estimate Hε and Hζ, assuming Case B recombination (Osterbrock & Ferland 2006). We then subtract the extrapolated Hε and Hζ fluxes from the corresponding blends to obtain He i λ3889 and [Ne iii] λ3968 fluxes separately.

The values of extinction we find using Balmer decrement and the Small Magellanic Cloud extinction curve from Weingartner & Draine (2001) are in the range A_V = 1.0–1.5 mag, in agreement with our previous estimates for type 2 quasars (Reyes et al. 2008) and with typical values in the literature for narrow-line regions of Seyfert galaxies (Bennert et al. 2006). Extinction values are higher (A_V ≃ 1.5 mag) for the two highest width composites than for the other ones (A_V ≃ 1.0–1.1 mag), somewhat reminiscent of the results by Veilleux (1991a) who finds higher extinction values for objects with stronger line asymmetries. As a function of luminosity, extinction decreases steeply and monotonically, from A_V = 1.7 to 0.9 mag. Because it is not clear that the extinction values derived from Balmer decrements apply to all other emission features (which may originate in a different spatial region), we do not apply extinction correction to any measurements, unless explicitly stated otherwise.

The optimal composite-making practices depend on the goal of the composite (Vanden Berk et al. 2001; Lacy et al. 2013). The typical per-spectroscopic-pixel errors of the SDSS spectra are the result of the overall plate reductions and thus are not too dissimilar from one object to the next, so in an error-weighted average all weights would be roughly the same. Therefore, we chose to use simple arithmetic averages to produce our composites (Lacy et al. 2013). It is encouraging that the [O iii]/[O ii] line ratios measured from the composites (log [O iii]/[O ii] = 0.67 − 0.74, depending on [O iii] width) are consistent with the one measured from individual spectra (mean = median and sample standard deviation 0.70 ± 0.23) and [O iii] λ4959/[O iii] λ5007 is within 1 per cent of its theoretical value in all composites. This means that our composite-making procedure does a reasonable job of preserving line ratios, both for cases when the intrinsic distribution of the line ratio is very narrow (e.g. [O iii] λ4959/[O iii] λ5007) and when it is fairly broad (e.g. [O iii]/[O ii]).

5.2 Wolf–Rayet features

One striking result is the appearance of a broad complex around He ii λ4686 in the composite with the highest w₈₀. This feature is somewhat reminiscent of the broad emission signatures of Wolf–Rayet stars (Brinchmann, Kunth & Durret 2008; Liu et al. 2009) – massive young stars which can radiatively drive outflows from their photospheres with velocities reaching 2000 km s⁻¹. A detection (or lack thereof) of young stellar populations is of particular importance in our work, since supernova explosions are capable of driving powerful outflows. Thus, a significant star-forming population could be responsible for at least some of the outflow activity, and therefore this feature deserves particular scrutiny.

The emission-line region photoionized by the quasar produces a multitude of forbidden and recombination (‘nebular’) features in the wavelength range between 4600 and 4750 Å (Liu et al. 2009). The broad lines produced by the Wolf–Rayet stars at these wavelengths are N v λ4613, N iii λ4640, C iii/iv λ4650 and He ii λ4686 (Brinchmann et al. 2008). Thus the challenge is to disentangle the two sets, made even more difficult by the fact that at the high velocities of the ‘narrow’ lines in our highest w₈₀ objects the velocities in the quasar-ionized set are not dissimilar from the velocities in the Wolf–Rayet set.

In Fig. 14 we zoom in on the wavelength range of interest. We find that only four nebular emission features ([Fe  iii] λ4658, He ii λ4686, [Ar iv] λ4711 and [Ar iv] λ4740) are sufficient to explain most of the observed flux. For each composite, we assume that these four features are associated with the narrow-line region of the quasar and have the same kinematic structure as Hβ. We take the Hβ velocity profile from the same composite and we adjust its amplitude to fit each of these four features. We see no evidence that the relative fluxes of the four lines vary from one composite to the next; thus the models shown in Fig. 14 have [Fe  iii]/He ii fixed at 0.15 and each [Ar iv]/He ii fixed at 0.20, making it essentially a one-parameter fit (the overall ratio of He ii to Hβ). No broad component in He ii is required to produce an adequate fit. There is some evidence in the last composite of some emission filling in between the two argon features, possibly due to [Ne iv] λ4725.

Figure 14.

Zoom in on the wavelength region that may contain Wolf–Rayet signatures for the five composite spectra made in bins of [O iii] line width, which is increasing from left to right. In each panel, the grey histogram shows the composite spectrum and the red curves (light grey) show Hβ velocity profiles from the same composite scaled to match the amplitudes of [Fe  iii] λ4658, He ii λ4686, [Ar iv] λ4711 and [Ar iv] λ4740. In all five panels, [Fe  iii]/He ii = 0.15 and each [Ar iv]/He ii = 0.2. The blue curve (solid black) shows the sum of the red profiles. The light blue curve (dotted) in the last panel shows the putative N iii λ4640 component.

Open in new tabDownload slide

Although most of the ‘broad’ feature turns out to be due to a superposition of relatively broad nebular lines, in the last panel we see excess emission centred at around 4640 Å. This is close to the wavelength of one of the more prominent Wolf–Rayet features, N iii λ4640. To estimate the flux of this feature, we model it using the same Hβ velocity profile (although in this case there is no particular physical reason to do so, since the N iii line is expected to reflect the kinematics of the stellar winds, rather than be produced in the extended low-density medium of the host galaxy and reflect the kinematics of quasar-ionized gas). We find a reasonable fit with N iii/He ii = 0.15, but we rather doubt the identification of this feature as N iii associated with the Wolf–Rayet stars in the host galaxy. The feature is noisy and is somewhat offset to the blue from the nominal expected centroid. When all spectra are co-added together into a single composite, the evidence for N iii disappears.

We assume for the moment that the N iii feature is in fact detected and estimate its median luminosity. To this end, we use the observed N iii/[O iii] ratios in the composite and the median [O iii] luminosity of the objects in the composite log L[O iii]/L_⊙ = 8.91 to find L(N iii) ≃ 3 × 10⁶ L_⊙. We then use N iii fluxes and starburst models from Schaerer, Contini & Kunth (1999) and Leitherer et al. (1999) to estimate star formation rates of ∼6 M_⊙ yr⁻¹. This is significantly lower than the estimates of star formation rates in type 2 quasar hosts by Zakamska et al. (2008) (a few tens  M_⊙ yr⁻¹). The difference between the two methods is not alarming because there are still significant discrepancies between modelled and observed fluxes of Wolf–Rayet features, especially as a function of metallicity (Brinchmann et al. 2008). But the low fluxes of Wolf–Rayet features in our composites – if they are detected at all – further reinforce our understanding that the star formation rates in type 2 quasar hosts are not adequate for producing the observed radio emission. At 10  M_⊙ yr⁻¹, using calibrations by Bell (2003) and Rosario et al. (2013) we find that typical star-forming galaxies would produce νL_ν[1.4 GHz] = 2.5 × 10³⁸ erg s⁻¹ and νL_ν[12 μm] = 2.9 × 10⁴³ erg s⁻¹, more than an order of magnitude below the values seen in our sample.

5.3 Line ratios as a function of line width and line luminosity

Fig. 15 summarizes our analysis of the line ratios as a function of [O iii] width and luminosity. We report the trends we see among the composite spectra at face value, but we caution that many of the lines we discuss are too weak to be detected in individual objects. As a result, in most cases it is not possible at the moment to evaluate whether the trends we see in the composite spectra accurately represent the trends in the population; a couple of exceptions are noted below.

Figure 15.

Line ratios as a function of line width (left) and [O iii] line luminosity (right) calculated from the composite spectra. On the horizontal axis, we use the median w₈₀ (in km s⁻¹) or the [O iii] line luminosity of the spectra contributing to the particular composite. Some curves are colour-coded (or shown with different line types) in no particular order to help distinguish them on the plot, and the bottom six curves are offset from each other by a small amount in the horizontal direction for the same reason; no actual change of line width or line luminosity is present.

Open in new tabDownload slide

For every composite, we compute errors in the line ratios using bootstrapping: we randomly resample (100 times) with replacement the set of spectra that contribute to each composite, recompute the composite, and recompute the line ratios. The error bars shown in Fig. 15 encompass 68 per cent of the resampled line ratios and represent the statistical error; the systematic uncertainties (those due to continuum subtraction, due to fitting procedure and due to assuming the same kinematic shape for all emission lines) are not included. Typical error bars for the line ratios of the brighter lines such as [O ii] and [Ne iii] are about ∼0.03 dex, reaching 0.2 dex for fainter lines such as [FeV ii]. For brighter lines the dispersion of line ratios within the population can be measured from individual objects and is typically ∼ 0.3-0.4 dex (Stern & Laor 2012); thus our bootstrapping procedure confirms that the error in the mean is consistent with |${\sim } 1/\sqrt{N}$| of the dispersion within the population (N = 114 is the number of the spectra in each composite). For fainter lines measurement uncertainties dominate over the statistical ones.

We find an increase of [S ii] λλ4069, 4076 (and to a lesser extent [N i] λλ5198, 5200) as a function of [O iii] width, which is our proxy for outflow velocity. Over the range of [O iii] widths that we probe, these lines increase, respectively, by factors of 2.6 and 1.6 (or by 0.4 and 0.2 dex) relative to [O iii]. Line ratios that vary with kinematics are often a tell-tale sign of shock excitation in starburst galaxies and low-ionization nuclear emission-line regions (Veilleux et al. 1994, 1995). [S ii] and [N i] trace warm weakly ionized gas phase, which is common in shocks that penetrate deep into clouds and that can produce extended low-ionization regions (Tielens 2005). The increase in the relative prominence of these lines suggests an increase in shock-ionization contribution, perhaps in direct response to an increase of the quasar wind velocity.

The [O iii] λ4363/[O iii] λ5007 ratio is ≃ 1.7 × 10⁻² and does not vary in a regular fashion with either [O iii] width or luminosity. This value is significantly higher than that predicted by photoionization models (Villar-Martín et al. 2008), and is higher still if corrected for extinction, suggesting some contribution from shocks. Combined shock- and photoionization models (Moy & Rocca-Volmerange 2002) can help explain this ratio, but because photoionization makes the dominant contribution to the [O iii] emission the lack of dependence on line kinematics is not too surprising.

Both [S ii] λλ4069, 4076 and especially [N i] λλ5198, 5200 appear to decline as a function of [O iii] luminosity, by factors of 1.3 and 2.4 (0.1 and 0.4 dex), respectively. A similar decline of stronger low-ionization lines such as [N ii] λ6583, [S ii] λλ6716, 6731 and [O i] λ6300 relative to [O iii] is seen in type 1 quasars by Stern & Laor (2012). Unfortunately, these particular lines are not accessible in our higher redshift objects, but applying their scalings to our lines and the range of [O iii] luminosities probed by our data, one would predict a decline by factors of 1.3–1.8, consistent with our observations. This decline suggests that the low-ionization regions are being destroyed as L[O iii] increases, likely as a result of the increase in the bolometric luminosity and thus the availability of ionizing photons. The same effect could also be produced by the increase in the opening angle of quasar obscuration, which would increase the size of the photoionized regions and lead to obliteration of the low-ionization regions by direct quasar radiation, but the strong positive [O iii]–infrared correlation suggests that the bolometric luminosity increase is the more likely driver for the observed correlations.

The [O ii]/[O iii] ratio decreases by a factor of 1.4 (0.14 dex) as a function of [O iii] luminosity, which was previously reported by Kim, Ho & Im (2006) and Silverman et al. (2009), with our values being in agreement with theirs. Quasar photoionization models produce poor fits to this ratio (Villar-Martín et al. 2008). Explanations usually involve a combination of high-ionization regions photoionized by the quasar and low-ionization regions photoionized by star formation (Kim et al. 2006), so the decrease in [O ii]/[O iii] with [O iii] luminosity then indicates the decrease in the relative contribution of the lower ionization regions, which is similar to the behaviour of [S ii] and [N i] lines. However, unlike [S ii] and [N i], [O ii]/[O iii] decreases with [O iii] line width (an effect tentatively seen by Steinhardt & Silverman 2013). We speculate that this trend may be an indirect hint that the quasar-driven wind of increasing velocity (as measured by [O iii] width) suppresses star formation in the host galaxy, a phenomenon which on occasion can be observed in action in individual sources (Cano-Díaz et al. 2012).

We do not find any noticeable trends in line ratios with radio luminosity, except that the lines are getting broader as the radio luminosity increases. This confirms our previous finding of the correlation between [O iii] width and radio luminosity (Fig. 10), but suggests that radio emission is not one of the principal drivers of the state of the ionized gas. Furthermore, there is no evidence that Wolf–Rayet features appear in higher radio luminosity objects, which is in qualitative agreement with our previous conclusion that radio emission is not dominated by the star formation in the host. We also examine composite spectra made in bins of infrared luminosity and find their behaviour similar to those made in bins of [O iii] luminosity, except that the width of the lines is a significantly stronger function of infrared than [O iii] luminosity, as we already know from Figs 10 and 11.

In the highest line width composite, the profile of [O ii] is visibly narrower than that of [O iii], and indeed we find w₉₀[O iii]=2132 km s⁻¹ (w₈₀ = 1498 km s⁻¹), whereas that of [O ii] (corrected for doublet structure using equation 1) is 1410 km s⁻¹ (1040 km s⁻¹). This confirms the trend we see in Fig. 8 and may mean that the [O ii]/[O iii] ratios calculated from the composite spectra are overestimated (as they are computed with a fixed velocity profile of [O iii]). For any lines other than [O ii], even with the high S/N of the composite spectrum it is not clear whether the kinematic structures are consistent from one line to the next. There is a hint that [N i] λλ5198, 5200 is narrower than [O iii] as well. We do not detect any obvious overall blueshifts of any features relative to [O iii].

Another curious tidbit is the increase in the He i λ3889/He i λ4471 ratio with [O iii] width. This ratio ranges from 1.2 to 2.1, which is significantly lower than the standard Case B values (Osterbrock & Ferland 2006). If we correct them for the extinction, we find values in the range from 1.8 to 2.6 (still somewhat lower than Case B) and increasing monotonically with w₈₀ (in the last w₈₀ bin, He i λ3889/He i λ4471 is likely underestimated because He i λ3889 starts to blend with [Ne iii] λ3869), without any trend with L[O iii]. It is possible that He i is affected by self-absorption (Monreal-Ibero et al. 2013) which declines as a function of line width.

5.4 Coronal lines

We detect two very high ionization features, [Fe V ii] λ3759 and λ5158, which show a tendency to increase with [O iii] width. Such lines, termed ‘coronal’, arise from ions with high-ionization potential, 99 eV in the case of Fe V ii. The most frequently discussed mechanism for their production is excitation by quasar photoionization fairly close to the nucleus, akin to an ‘intermediate’ region between the very compact broad-line region and a more extended narrow-line one (Grandi 1978; Gelbord, Mullaney & Ward 2009; Mullaney et al. 2009; Müller-Sánchez et al. 2011). In photoionization models, high electron and photon densities are typically required to produce coronal lines, which explains why they arise close to the nucleus, e.g. at sub-pc scale in the case of a Seyfert galaxy modelled by Mullaney et al. (2009). These authors propose that the inner edge of the obscuring material is the plausible location of the coronal emission regions. Similarly, the rapid variability of coronal line emission in another Seyfert galaxy suggests that the size of the emission region is just a few light years (Komossa et al. 2008).

Such small distances present a challenge to the study by Rodríguez-Ardila et al. (2011) who find that the coronal lines are detected at similar rates in type 1 and type 2 objects, indicating that they originate outside the obscuring material, although Rodríguez-Ardila et al. (2011) still require very high densities, n_e = 10⁸-10⁹ cm⁻³. If coronal lines originate somewhere near the inner edge of the obscuring material, depending on the orientation we may see them in an occasional type 2 object, but overall one would expect to see lower detection rates of coronal lines in type 2 objects than in type 1s in this scenario.

It would be surprising to see a region close to or inside of the obscuring material in type 2 quasars, and even more surprising to see its flux correlate with the kinematics of much lower ionization, much more extended gas. Some recent models demonstrate the possibility of producing coronal lines on the surface of narrow-line region clouds via photoionization from the quasar (Stern, Laor & Baskin 2014), and while this approach solves the problem of visibility of coronal lines in type 2 quasars, it may not explain the dependence on line width. Therefore, we suggest that instead of circumnuclear photoionization in type 2 quasars some or most of the coronal lines are due to shocks, possibly those produced when the radiatively accelerated wind from the quasar runs into dense clouds.

We draw inspiration from the observations by Mazzalay et al. (2013) who find that the extended ( ≲ 170 pc) coronal line emission in NGC 1068 is likely dominated by shock excitation. These authors find that coronal line ratios are more easily explained by shock models than by photoionization models. Furthermore, they find that the kinematic structure of coronal lines is similar to that of several lower ionization features. Most intriguingly, the spatial and kinematic distribution of coronal lines is similar to that of very low ionization lines, such as [Fe ii], characteristic of shocks propagating into partly ionized medium. The latter point is particularly interesting as we observe the simultaneous rise of coronal lines and [S ii] and [N i] in our sources.

In NGC 1068, coronal lines are co-spatial with the radio jet, and thus Mazzalay et al. (2013) suggest that jet-driven shocks may play a role in line excitation. This is also a possibility for our sources and for type 1 quasars (Husemann et al. 2013; Mullaney et al. 2013), especially in light of the fact that radio emission is strongly correlated with line kinematics (Section 4.1). However, we do not favour this scenario. While radio luminosity increases with [O iii] widths, and so do the fluxes of [FeV ii], [S ii] and [N i], there is no evidence that the latter lines increase in the composites made in bins of radio luminosity. What this suggests to us is that the primary driver for all these correlations is not the power of the radio emission, but the kinematics of the gas. Higher outflow velocity implies higher velocities of shocks driven into dense clouds, which enhances coronal emission, while the radio emission follows as discussed in Section 6.1.

What sort of conditions would be required to produce coronal lines by driving low-density wind into high-density clouds? The typical shock velocities inside the shocked clouds implied by the observed coronal line ratios in NGC 1068 are v_shock = 300–1000 km s⁻¹ for cloud densities n_cloud = 100–300 cm⁻³ (Mazzalay et al. 2013). To produce the same values in our picture, the incident wind velocity must be even higher, |$v_{\rm wind}\sim v_{\rm shock}\sqrt{n_{\rm cloud}/n_{\rm wind}}$|⁠, where n_wind is the much lower density of the volume-filling component of the wind, so the required v_wind can be several thousand km s⁻¹. Such high velocities are not achievable at large distances in low-luminosity active galaxies, where the radiatively driven wind would get quenched by the interaction with the interstellar medium, but could be typical of quasar-driven winds (Moe et al. 2009; Faucher-Giguère & Quataert 2012; Zubovas & King 2012) even on galaxy-wide scales (Arav et al. 2013; Liu et al. 2013b).

6 DISCUSSION

6.1 The origin of the line kinematics/radio/infrared relationships

The tantalizing relationship between radio and infrared luminosities and the [O iii] line kinematics may contain interesting clues about the structure and the driving mechanism of ionized gas outflows in quasars, but its interpretation is not straightforward. For quite some time, it was thought that the extended ionized gas around active galaxies was much more likely to be found around radio-loud objects than around radio-quiet ones (Stockton & MacKenty 1987). This was in line with theoretical ideas: powerful relativistic jets can inflate overpressured cocoons which in turn sweep up galactic and intergalactic medium plausibly resulting in narrow-line emission (Begelman & Cioffi 1989).

With more sensitive observations, extended ionized gas emission has now been detected around radio-loud and radio-quiet quasars alike, indicating that gas nebulae around both types have roughly similar sizes (Liu et al. 2013a,b; Liu, Zakamska & Greene 2014), which reach 10 kpc or greater at the highest luminosities (Greene et al. 2011; Hainline et al. 2013, 2014; Liu et al. 2013a; Harrison et al. 2014). The major difference appears to be that the nebulae around radio-loud objects tend to be more disturbed (lumpy) or more elongated, whereas the nebulae around radio-quiet objects are smooth and featureless (Liu et al. 2013a). The gas kinematics are disturbed and show outflow signatures not only in the central brightest parts of the nebulae, but all the way out to several kpc from the nucleus (Greene et al. 2011; Liu et al. 2013b).

The recent discovery of the relationship between gas kinematics and radio emission within the radio-quiet population goes to the heart of the origin of the outflows in the majority of quasars, as well as of the origin of radio-quiet emission itself. Within the last year, there have been several new results on these issues (Condon et al. 2013; Husemann et al. 2013; Mullaney et al. 2013). In this section, we discuss pros and cons of various mechanisms proposed to explain the radio-quiet emission of quasars and the correlation between radio emission and gas kinematics.

The observed gas kinematics/radio correlation is reminiscent of the classical work on narrow-line kinematics by Veilleux (1991a,b,c) who found a correlation between [O iii] λ5007 line width and radio luminosity within a sample of 16 nearby Seyfert galaxies. In these and other nearby objects, high-resolution radio observations sometimes reveal a strong relationship between the position angle and morphology of the radio emission and the ionized gas emission, especially for the kinematically disturbed gas component (e.g. Veilleux 1991c; Bower et al. 1995; Capetti, Axon & Macchetto 1997; Falcke, Wilson & Simpson 1998; Capetti et al. 1999). In some cases the jet appears to be at least in part responsible for the physical conditions in the narrow-line region (e.g. Kraemer, Ruiz & Crenshaw 1998). Interestingly, Leipski et al. (2006) find a correlation between the size of the narrow-line region and the size of the radio source in nearby radio-quiet type 1 quasars. Thus, we first consider the possibility that the gas kinematics/radio relationship is established via direct interactions between the jet and the gas.

In sources where the jet and the ionized gas emission are aligned, it is tempting to postulate that radio jets are either driving the outflow or provide shock ionization, making them responsible for the morphology and/or excitation of the ionized gas, which may provide some basis for the radio/kinematics correlation. However, there may be a certain bias in reporting only those objects where high-resolution radio observations reveal a jet, especially if its position angle is not far off the orientation of the line emission. Far from every source shows evidence of jets in high-resolution radio observations (Veilleux 1991c; Ho & Ulvestad 2001). In large unbiased samples the alignment between jets and ionized gas is statistically weak at best (Privon et al. 2008), although some authors find a stronger relationship (Schmitt et al. 2003), especially in compact radio sources (Shih, Stockton & Kewley 2013). This remains a difficult issue because the morphology and the position angle of the radio emission and the ionized gas emission are a function of scale, and it is unclear which scales are relevant. Additionally, in principle the direction of photoionization by the active nucleus (which often determines the orientation of the nebula) and the jet (which may be responsible for the disturbed kinematics) can be misaligned, further complicating the picture.

We examined 20 objects from our sample that are located in Stripe 82 and were mapped by the Very Large Array by Hodge et al. (2011). These observations are deeper than FIRST (rms = 52 versus 130 μJy beam⁻¹) and have higher spatial resolution (FWHM = 1.8 versus5 arcsec), better matched to the sizes of the quasar host galaxies at the typical redshifts of our sample (1 arcsec at z = 0.5 is 6.1 kpc). Of the 20 sources, one object is a giant Fanaroff–Riley type II radio galaxy in both data sets; 13 are point sources both in FIRST and in Stripe 82; and the remaining 7 are not detected in FIRST, but are all detected in Stripe 82 as point sources with fluxes between 0.4 and 0.7 mJy, i.e. just a factor of a couple below the catalogue limit of the FIRST survey. None of the 19 objects showed resolved morphology in the higher resolution observations. Among type 2 quasars, Lal & Ho (2010) find a mix of unresolved and slightly resolved sources at 0.8 arcsec resolution, but only a small fraction (∼1/7) with linear jet-like signatures. The radio emission of the ∼10 kpc ‘superbubble’ quasar outflow (Greene et al. 2012) is compact (<1 kpc; Wrobel, private communication). Overall we find it unlikely that large-scale ( ≳ 5 kpc) jets are responsible for the correlation between gas kinematics and radio emission we see in our sample.

A more promising hypothesis is that in the radio-quiet objects the radio emission is due to a compact jet ( ≲ 1 kpc) not resolved in the current radio observations which is driving the outflow of ionized gas (Spoon & Holt 2009; Husemann et al. 2013; Mullaney et al. 2013). Because the kinematics of the gas is disturbed on galaxy-wide scales (out to ∼10 kpc; Liu et al. 2013b), the effect of the jet on the gas would need to be indirect in this scenario: the jet would need to inject the energy into the interstellar medium of the galaxy on small unresolved scales, with the energy then converted into a wide-angle outflow that engulfs the entire galaxy as observed. Although the observational test of this hypothesis seems straightforward (every radio-quiet quasar should have a jet), in practice testing this paradigm is difficult, in part because it is not clear that every extended radio structure represents a jet. Furthermore, radio emission of radio-quiet active galaxies tends to be compact (Nagar et al. 2003), so it is often necessary to zoom in all the way into the central several pc in order to resolve the radio emission (Nagar, Falcke & Wilson 2005). On those scales, even when a linear structure seems unambiguously jet like, it is not clear that it is sufficient to power an energetic galaxy-wide outflow.

A completely different approach is followed by Condon et al. (2013) who argue, on the basis of the shape of the radio luminosity function, that the radio emission in radio-quiet quasars is mostly or entirely due to star formation in their host galaxies. In our objects, in favour of this hypothesis is the strong correlation between the infrared and the radio fluxes, which furthermore appear to lie on the extension of the locus of normal star-forming galaxies (Fig. 11). However, at the luminosity of the radio-quiet subsample of log (νL_ν[1.4 GHz], erg s⁻¹) = 40.0 ± 0.7 (mean = median and standard deviation), the median star formation rate suggested by the classical radio/star formation correlation (Helou, Soifer & Rowan-Robinson 1985; Bell 2003) would be an astonishing 400  M_⊙ yr⁻¹. Although the star formation rates of type 2 quasar hosts are among the highest in the population of active galaxies, in the ballpark of a few tens  M_⊙ yr⁻¹ (Zakamska et al. 2008), they fall short of the values required to explain the observed radio luminosity, as was also noticed by Lal & Ho (2010). Furthermore, as was already noted above, the mid-infrared fluxes of our objects shown in Fig. 11 are dominated by the emission of the quasar-heated dust. Thus, both radio luminosity and mid-infrared luminosity in Fig. 11 are unlikely to be due to star formation.

Although compact jets cannot be ruled out by the existing data, we put forward an alternative suggestion for the origin of the radio-quiet emission. We propose that luminous quasars radiatively accelerate winds (Murray et al. 1995; Proga, Stone & Kallman 2000) which then slam into the surrounding medium and drive shocks into the host galaxy. This shock blast initiated by the quasar-driven wind accelerates particles, just like what happens in a supernova remnant, which then produce synchrotron emission (Stocke et al. 1992; Jiang et al. 2010; Faucher-Giguère & Quataert 2012; Zubovas & King 2012). To estimate the energetics of the wind required in this scenario, we assume for the moment that the efficiency of converting the kinetic energy of the outflow L_wind into radio synchrotron emission is similar in starburst-driven and quasar-driven winds. In a galaxy forming ψ  M_⊙ yr⁻¹ worth of stars, the kinetic energy of the starburst-driven wind is 7 × 10⁴¹ψ erg s⁻¹ (Leitherer et al. 1999). The same galaxy produces νL_ν[1.4 GHz] = 2.5 × 10³⁷ψ erg s⁻¹ worth of radio emission (Bell 2003), with efficiency of 3.6 × 10⁻⁵ of converting kinetic energy into radio luminosity.

Applying the same efficiency to our objects, we find that to reproduce our median radio luminosity of 10⁴⁰ erg s⁻¹ we need L_wind = 3 × 10⁴⁴ erg s⁻¹. The bolometric luminosities of our objects are poorly known, but using scalings presented by Liu et al. (2009) we can estimate that for the median [O iii] luminosity of our sample (log L[O iii]/ L_⊙ = 8.9) the bolometric luminosity is 8 × 10⁴⁵ erg s⁻¹. Thus, the median fraction of the bolometric luminosity converted to the kinetic luminosity of the wind is 4 per cent, in rough agreement with our previous estimates based on observations of spatially resolved winds (Liu et al. 2013b). These ideas are further supported by our observation that νL_ν[1.4 GHz] has a close-to-quadratic dependence on [O iii] width (Fig. 10), suggesting that the relationship between the kinetic energy of the ionized gas and radio luminosity is linear.

In such heavily obscured quasars we cannot directly investigate the effect of the optical luminosity or even X-ray luminosity on gas kinematics or radio emission, which is an interesting direction to pursue using samples of type 1 quasars. In these objects, there is a known correlation between radio luminosity and their accretion power derived from ultraviolet luminosity (Falcke, Malkan & Biermann 1995). While this correlation is usually interpreted as a physical connection between accretion and resulting jets, an alternative explanation is that the accretion power leads to winds which in turn generate radio emission as described above without involving jets.

6.2 Radiatively driven quasar winds

In Fig. 16 we present a schematic of the components of the narrow-line regions of quasars discussed in this paper. The quasar is in the centre and is characterized by its bolometric luminosity L_bol. Obscuring material is distributed anisotropically around it (dark grey ‘torus’), with Ω being the opening angle of the obscuration which determines the directions that can be photoionized by the direct quasar radiation. The relationship between the two values L_bol and Ω is not well known. Demographic studies of quasars continue to disagree on whether the obscured fraction stays constant or decreases as a function of luminosity (Treister, Krolik & Dullemond 2008; Lawrence & Elvis 2010), but it is likely that at a given luminosity there is a wide distribution of Ω (Zakamska et al. 2006).

Figure 16.

A cartoon of our model for the narrow-line region of a quasar (central five-point star). The quasar is surrounded by circumnuclear obscuring material (dark grey ‘torus’) whose opening angle is Ω. The wind propagates more isotropically than the ionizing radiation from the quasar which is concentrated in bi-cones marked ‘photoionized’. Optical emission lines are produced by dense clouds (small circles).

Open in new tabDownload slide

Another component of the model is the quasar-driven wind. The wind may be initially driven anisotropically, e.g. close to the equatorial plane of the accretion disc (Murray et al. 1995; Proga et al. 2000), and when it is first launched near the quasar, it is very fast, with velocities up to ∼0.1c. The wind runs into the interstellar medium close to the nucleus and interacts with it, producing shock waves that propagate through the interstellar medium of the galaxy. If this medium is clumpy, the propagation of the wind proceeds along paths of least resistance, and the shape of the large-scale wind is determined largely by the distribution of the interstellar medium (Wagner et al. 2013), rather than by the initial anisotropies. Thus, the wind is shown propagating in all directions, even those that are affected by the circumnuclear obscuration.

The wind velocity v_wind likely varies as a function of the location in the galaxy and is not well-known observationally. In our data, the line-of-sight velocity dispersion of the optical emission lines is characteristic of the velocities of the clouds v₀ (light grey), not the velocity of the wind. The acceleration of clouds and their survival as they are impacted by the wind is a topic of active research (Mellema, Kurk & Röttgering 2002; Cooper et al. 2009; Alūzas et al. 2012; Barai, Proga & Nagamine 2012; Faucher-Giguère, Quataert & Murray 2012), but it seems natural that the velocity of the wind should be at least as high as the velocity of the clouds that it accelerates (v_wind ≳ v₀) and that they should be strongly related. Theoretical arguments suggest wind velocities of v_wind ≳ 1000 km s⁻¹ (King 2005; Faucher-Giguère & Quataert 2012; Zubovas & King 2012), in rough agreement with our inferred values of v₀.

We do not observe v₀ directly; rather, we infer these values from the line-of-sight velocity width of the emission lines. The increase in shock-diagnostic lines [S ii], [N i] and possibly [FeV ii] as a function of [O iii] width provides evidence that the observed velocity width is a good proxy for the outflow velocity. Furthermore, the [O iii] width positively correlates with mid-infrared luminosity (Fig. 11), which we take to suggest that the outflow is ultimately driven by the radiative pressure near the quasar (Murray et al. 1995; Proga et al. 2000). After the winds are launched, they interact with the interstellar medium in complex ways, which perhaps explains why the correlation between infrared luminosity and [O iii] velocity width (which our observations probe on galaxy-wide scales; Liu et al. 2013b) has a large scatter.

7 CONCLUSIONS

7.1 Gas kinematics and outflows

In this paper we study the kinematics of ionized gas emission in 568 luminous obscured quasars from the sample of Reyes et al. (2008), primarily using their [O iii] λλ4959, 5007 emission lines. For every object, we determine a set of non-parametric measures of line asymmetry, velocity width and shape. We find that objects with blueshifted emission and line asymmetries are more prevalent in our sample, which we take as a signature that outflows in which the redshifted part is affected by dust extinction are common in our sample. The velocity widths we see in our sample (median w₈₀ = 752 km s⁻¹ and w₉₀ = 1060 km s⁻¹, max w₈₀ = 2918 km s⁻¹ and w₉₀ = 4780 km s⁻¹) are much higher than those in starburst galaxies (median w₉₀ ≃ 600 km s⁻¹; Rupke & Veilleux 2013a; Hill & Zakamska 2014).

Using the conversion between the line-of-sight velocity width of the line and the outflow velocity from Section 3.1, we can estimate outflow velocities as v₀ ≃ w₈₀/1.5. We find that for about half of the objects in our sample, v₀ ranges from ∼500 to ∼2000 km s⁻¹. It is likely that in every object a range of cloud velocities is present; these estimates should be thought of as the median velocities in each source. From the observations presented here, it is not known which spatial scales dominate these estimates, but from our spatially resolved observations (Liu et al. 2013b) it is clear that at least some of the clouds maintain similar velocities all the way out to ∼10 kpc from the quasar.

Neither the overall width nor the width of the narrower component in the multi-Gaussian line decomposition correlate with the velocity dispersion of the host galaxy. We therefore find no evidence for significant amounts of gas in dynamical equilibrium with the host (e.g. a component associated with the rotating galaxy disc). The higher velocity gas tends to be on average blueshifted by ∼100 km s⁻¹ relative to the lower velocity gas. This may mean that the higher velocity gas is found on somewhat smaller spatial scales which are more prone to dust extinction. This conclusion is in line with the slight observed decrease of line-of-sight velocity dispersions of the gas as a function of the distance from the nucleus seen in spatially resolved observations of quasar winds (Liu et al. 2013b). The effect is very subtle (3 per cent decrease per each projected kpc), and overall high gas velocities are maintained over the entire host galaxy (Liu et al. 2013b). Other strong lines, such as [O ii] and Hβ, also show outflow signatures, but significantly more mild than those seen in the [O iii] line. In the objects with the most extreme [O iii] kinematics, [O ii] is much narrower than [O iii], with the kinematics of Hβ being in between the two.

The increase in shock-diagnostic lines [S ii] λλ4069, 4076 and [N i] λλ5198, 5200 as a function of [O iii] width strongly support the use of [O iii] width as a proxy for the outflow velocity. [O ii] declines both with [O iii] width and [O iii] luminosity, suggesting that the [O ii]-emitting regions are overionized by photoionization and quenched by increasing shock velocity. We find evidence for shock-ionization contribution to coronal lines [FeV ii].

7.2 Outflows and radio emission

[O iii] velocity width correlates strongly with radio luminosity. We suggest that the radio emission is a by-product of the outflow activity, with particles accelerated on the shock fronts as the radiatively driven quasar wind propagates into the interstellar medium of the host galaxy [this mechanism was first hypothesized by Stocke et al. (1992) to explain the low levels of radio emission in broad absorption-line quasars]. The median radio luminosity in our sample, νL_ν[1.4 GHz] = 10⁴⁰ erg s⁻¹, requires kinetic energy of the outflow of L_wind = 3 × 10⁴⁴ erg s⁻¹ if the efficiency of conversion in quasar-driven winds is similar to that in supernova-driven winds. The bolometric luminosities of our objects are poorly known, but we estimate that the median value for our sample is 8 × 10⁴⁵ erg s⁻¹ and thus the kinetic luminosity of the wind is 4 per cent of the bolometric power. This idea is further reinforced by the approximately quadratic dependence νL_ν[1.4 GHz] ∝ (width[O iii])², which implies a linear relationship between radio luminosity and kinetic energy of the outflow. [O iii] velocity is positively correlated with the infrared luminosity, which suggests that the ionized gas outflow is ultimately driven by the radiative pressure close to the accreting black hole.

Another possible explanation for the [O iii] width/radio correlation is that the mechanical energy of the relativistic jet (which has not yet broken out of its host galaxy) is used to heat an overpressured cocoon (Begelman & Cioffi 1989) which then launches the wind of ionized gas (Mullaney et al. 2013). As long as the jet is still compact, it appears as an unresolved radio core in FIRST observations. We use the scalings between the core radio luminosity and the jet kinetic energy by Merloni & Heinz (2007) to estimate how much kinetic power would be required to produce the observed median radio luminosity. We find L_jet = 3 × 10⁴⁴ erg s⁻¹, a value which is almost identical to the wind kinetic power obtained in the previous paragraph, even though they were estimated using completely different methods. This is likely not a coincidence: the efficiency of conversion of mechanical luminosity into synchrotron radiation is determined by the fraction of energy that can be converted into relativistic particles on shock fronts, which is likely independent of the origin of these shocks.

The similarity of energy requirements underscores the difficulty of distinguishing between these two mechanisms. In the jet scenario one expects to see a jet in the high-resolution radio observations, whereas in the wind scenario radio emission is more diffuse and is present everywhere where shocks are propagating; however, in practice the collimated part of the jet does not dissipate energy very efficiently and so can be hard to detect. In both models, the morphology of the ionized gas depends much more strongly on the distribution of the interstellar medium than on the exact driving mechanism (Gaibler, Khochfar & Krause 2011; Wagner, Bicknell & Umemura 2012; Wagner et al. 2013). Thus, the morphology of ionized gas is not necessarily a useful clue.

Radio spectral index could be a useful measurement for identifying recent or on-going particle acceleration (flatter synchrotron spectra for freshly accelerated, more energetic particles), but again particle acceleration may be happening in both scenarios. Interestingly, Lal & Ho (2010) find flatter spectral indices in type 2 quasars than expected for jets in the standard geometry-based unification models, so a combination of spectral and morphological investigations in the radio may be worth pursuing further.

One uncomfortable consequence of the jet scenario is that it requires for every radio-quiet quasar to have a powerful jet, with only a minority of them being active long enough to to break out of the galaxy (otherwise there would be too many extended radio sources). Another problem is that the jet scenario provides no ready explanation for the correlation between [O iii] width and infrared luminosity, unless jets are also directly connected to the accretion power (Falcke et al. 1995).

In the wind scenario, if the initial distribution of matter is roughly spherically symmetric, so will be the ionized gas emission and the radio emission. In a disc galaxy, once the size of the wind reaches the vertical scaleheight of the disc it propagates largely along the path of least resistance perpendicular to the disc (Wagner et al. 2013), forming two symmetric bubbles on either side. Ionized gas emission is concentrated in shells and filaments, whereas the radio emission is filling the bubbles. In low-resolution data, one would see both the line emission and the radio emission oriented roughly in the same direction, so distinguishing this mechanism from jet-induced outflow requires high-quality observations.

Such bubbles are directly seen in the radio in external galaxies (Cecil et al. 2001; Hota & Saikia 2006) and in our own Milky Way, where the bubbles are known as ‘the microwave haze’ (Finkbeiner 2004; Su, Slatyer & Finkbeiner 2010). These structures are also seen in X-rays which often closely trace the morphology of the radio emission and/or the ionized gas filaments (Cecil et al. 2001; Croston et al. 2008; Wang et al. 2009). In our Galaxy not only X-rays are observed (Sofue 2000), but the structures as seen in gamma-rays and are known as ‘the Fermi bubbles’ (Su et al. 2010). For comparison, the mechanical luminosity necessary to inflate the bubbles in NGC 3079 (Cecil et al. 2001) is ∼30 times lower than our median L_wind, whereas that inferred for ‘the Fermi bubbles’ in the Milky Way is 10⁴-10⁵ times lower than our L_wind. It is interesting that all of the examples above host jets (Cecil et al. 2001; Croston et al. 2008; Wang et al. 2009; Su & Finkbeiner 2012), but it remains unclear whether this is universally true and whether they are contributing most of the required power.

7.3 Feedback in low- and high-luminosity active galaxies

Our investigation focuses on the most luminous type 2 quasars at z ≲ 1. In these sources, we find strong evidence for quasar-driven winds on galaxy-wide scales, for radio emission associated with these winds, and for ionized gas velocities in excess of escape velocity from the galaxy. In lower luminosity active galaxies, many previous studies demonstrated that ionized gas is in dynamical equilibrium with the host galaxy and that radio emission is likely a bi-product of star formation processes. Thus, quasar-driven feedback may be present above some threshold in luminosity and absent below this threshold. We estimate this threshold by observing that in Fig. 11, our sources separate well from star-forming galaxies at about νL_ν[12 μm] = 2 × 10⁴⁴ erg s⁻¹. Using bolometric corrections from Richards et al. (2006), Liu et al. (2009) and Liu et al. (2013b), we estimate the corresponding bolometric luminosity to be L_bol = 3 × 10⁴⁵ erg s⁻¹ (uncertainty ±0.4 dex).

Given the crudeness of this estimate, we consider it to be consistent with the value of L_bol = 2.4 × 10⁴⁵ erg s⁻¹ (uncertainty ±0.3 dex) suggested by Veilleux et al. (2013) based on incidence of molecular outflows in ULIRGs (Veilleux et al. 2013). In principle the threshold value should be dependent on the depth of the potential (and therefore the stellar velocity dispersion), as well as the amount of gas that needs to be accelerated (Zubovas & King 2012), and it will be interesting to see whether these ideas are borne out in future analyses.

NLZ would like to acknowledge useful conversations with M. J. Collinge (who suggested the diffuse nature of radio emission) and S. Tremaine (who suggested the analogy with Fermi bubbles), as well as with B. Groves, J. Krolik, D. Kushnir and J. Ostriker. The authors are grateful to the anonymous referees, as well as to H. Falcke, C. Norman, G. Richards, H. Spoon, J. Stern, J. Stocke, S. van Velzen and S. Veilleux for constructive comments during the review process. NLZ is thankful for the continued hospitality of the Institute for Advanced Study (Princeton) where part of this work was performed.

REFERENCES