© ESO 2020

1. Introduction

Quasars represent some of the most energetic sources in the Universe that may regulate the gas flows in and out their host galaxies. A manifestation of the impact that super-massive black holes (SMBHs) may have on the galaxy are the well-established scaling relations between the black hole and host galaxy properties such as the black hole mass, M_BH, versus galaxy mass, M_*, (e.g. Magorrian et al. 1998; Läsker et al. 2016; Schutte et al. 2019) and the M_BH versus stellar velocity dispersion, σ_*, relations (e.g. Gebhardt et al. 2000; Batiste et al. 2017; Caglar et al. 2020).

One promising physical phenomenon to link the growth of the SMBH and the evolution of its host is that of fast winds (> 1000 km s⁻¹) launched from the accretion disc of active galactic nuclei (AGN, e.g. King 2003; Begelman 2003; Menci et al. 2008; Zubovas & King 2012; Faucher-Giguère & Quataert 2012; Choi et al. 2014; Nims et al. 2015; Hopkins et al. 2016). These winds are hypothesised to shock against the surrounding gas and drive outflows that propagate out to large distances from the AGN, heat the interstellar medium (ISM), and potentially eject large amounts of gas from the system (e.g. Ishibashi & Fabian 2016; Zubovas 2018). Such fast winds are now observed in a vast number of AGN host galaxies at both low and high redshift, in different phases of gas- neutral phase using sodium absorption lines (e.g. Krug et al. 2010; Rupke & Veilleux 2011; Cazzoli et al. 2016; Concas et al. 2019; Roberts-Borsani 2020), cold molecular gas phase, using different transitions of CO, HCN and [C II], for instance (e.g. García-Burillo et al. 2014; Tadhunter et al. 2014; Feruglio et al. 2017; Aladro et al. 2018; Michiyama et al. 2018; Zschaechner et al. 2018; Aalto et al. 2019; Husemann et al. 2019; Cicone et al. 2020; Veilleux et al. 2020), warm and hot molecular gas phase, using transitions in the mid- to near-infrared (Veilleux et al. 2009; Davies et al. 2014; Hill & Zakamska 2014; Riffel et al. 2015; Emonts et al. 2017; May et al. 2018; Petric et al. 2018; Riffel et al. 2020) and ionised gas phase observed using the rest-frame optical emission lines such as the forbidden transitions of [O III] λ5007 (e.g. Harrison et al. 2014; Kakkad et al. 2016; Zakamska et al. 2016; Fiore et al. 2017; Venturi et al. 2018; Baron & Netzer 2019; Coatman et al. 2019; Förster Schreiber et al. 2019). Due to the higher surface brightness of the ionised gas traced by the forbidden transition [O III] λ5007 relative to other optical transitions (e.g. [O II] λ3727, [S II] λ6716), outflows in this phase can be studied in detail for a large number of galaxies. In AGN, these forbidden ionised transitions such as [O III], [N II], and [S II] are emitted from the extended narrow line region (ENLR), making these transitions ideal to trace kiloparsec-scale ionised gas outflows from the AGN (e.g. Bennert et al. 2002; Hainline et al. 2014; Dempsey & Zakamska 2018).

It is particularly important to study the impact that such galactic-scale AGN-driven outflows may have on their host galaxies at z ∼ 2 where both the volume-averaged star formation rate and the black hole growth rate peak (e.g. Shankar et al. 2009; Madau & Dickinson 2014; Curran 2019; Wilkins et al. 2019; Tacconi et al. 2020). Tremendous progress has been made through integral field unit (IFU) spectroscopy, which provides spatially resolved information on the structure and extension of the outflows (e.g. Riffel et al. 2013; McElroy et al. 2015; Thomas et al. 2017; Revalski et al. 2018; Davies et al. 2019; Radovich et al. 2019). Compared to the classical narrow band imaging and/or slit spectroscopy, IFU spectroscopy allows us to identify the emission from the host galaxy by subtracting the contribution from the AGN. A few IFU studies have also claimed the presence of “negative” as well as “positive” feedback in the presence of outflows, that is, ionised outflows suppressing as well as enhancing star formation within the AGN host galaxies (e.g. Cano-Díaz et al. 2012; Cresci et al. 2015; Carniani et al. 2016; Maiolino et al. 2017; Gallagher et al. 2019); however, the full intepretation of these results can be complicated by effects such as dust obscuration (e.g. Whitaker et al. 2014; Brusa et al. 2018; Scholtz et al. 2020). Most of the literature is however focused on a limited number of targets, which have been pre-selected to have a higher probability of showing the presence of outflows such as targets selected based on their colours, high Eddington ratio, or high luminosity (e.g. Brusa et al. 2015; Perna et al. 2015; Kakkad et al. 2016; Bischetti et al. 2017; Perrotta et al. 2019). The lack of studies in a wider parameter space of the properties of the AGN and their host galaxies has so far prevented the assessment in a systematic way of a possible trend between the activity of the AGN itself (e.g. quantified by its luminosity) and the presence of outflows, and, consequently, whether all outflows have an impact on their host galaxies. We are therefore in need of an unbiased sample where a wider range in Eddington ratio and/or bolometric luminosity is used for follow-up outflow studies.

We are currently in an era of large IFU surveys of galaxies (Sánchez et al. 2012; Bundy et al. 2015; Bryant et al. 2015; Stott et al. 2016; Förster Schreiber et al. 2018; Wisnioski et al. 2019; den Brok et al. 2020). Such statistical samples are now able to eliminate the selection biases from previous studies and give an overall picture of the ISM dynamics in the presence of both star formation and AGN processes. Recent work targeting star forming galaxies at high redshift, such as the SINS/zC-SINF survey (e.g. Förster Schreiber et al. 2018; Davies et al. 2019) with SINFONI indicate that the majority of galaxies (∼70%) show ordered disc rotation, while the rest of the targets either show a turbulent disc structure or the presence of outflows in addition to the disc rotation. The mass loading factor, which is the ratio between the outflow mass and the star formation rate, is correlated with the level of star formation within the host galaxies. The outflow fraction in the SINS/zC-SINF survey is similar to that of the KMOS-3D survey (Förster Schreiber et al. 2019), where 30% of the 599 observed targets show the presence of outflowing gas, inferred as being driven by both star formation and AGN processes. Star formation-driven winds, on average, show low mass loading factors when compared to the AGN-driven outflows. Coupled with the reported results in the KROSS survey (e.g. Swinbank et al. 2019), a fraction of the total gas mass is believed to escape in the low mass galaxies, while all of the outflowing gas is retained in galaxies at high masses. IFU surveys of star forming galaxies come to a common conclusion that the outflow velocities are typically enhanced in the presence of an AGN.

Among AGN surveys, Leung et al. (2019) target optically-selected AGN using single-slit spectroscopy from the MOSDEF survey with bolometric luminosities in the range 10⁴⁴ − 10⁴⁷ erg s⁻¹ at z ∼ 1.4–3.8 and find outflows in 17% of their sample of ∼160 AGN. Moreover, it is claimed that the ionised gas mass outflow rates correlate positively with the luminosity of the AGN, but do not depend on the galaxy stellar mass, which is in contrast with the findings of Förster Schreiber et al. (2018) concerning star forming galaxies. Active galactic nuclei from the MOSDEF survey have higher incidence of outflows compared to redshift and stellar mass matched star forming galaxies. Although MOSDEF uses single-slit spectroscopy, which, as mentioned earlier, has its own limitations, a similar result was found by Harrison et al. (2016) with the KASHz survey targeting X-ray selected AGN (2–10 keV) at z ∼ 1.1–2.5 with KMOS. The more luminous KASHz AGN are found to be more likely to host outflows with velocities > 600 km s⁻¹, and these ionised gas velocities are ten times more prevalent in the AGN host galaxies than star forming galaxies at similar redshifts and similar Hα luminosity distributions.

Most of the current AGN surveys at z >  1 described above are based on seeing-limited observations. Consequently, both host galaxies and the possible outflows are usually spatially unresolved, which prevents one from obtaining a complete understanding of how extended these outflows are. In addition, the presence of unresolved outflows implies that further assumptions have to be made (e.g. on the scale and morphology), which results in larger uncertainties in derived quantities such as mass loading factor and outflow kinetic power. Therefore, the way forward is to use adaptive optics (AO) assisted observations to resolve smaller physical scales also at z >  1 (e.g. Davies et al. 2020a). A higher spatial resolution enables one to study the impact of AGN from kiloparsec scales to Mpc scales. In this paper, we present the first results from the Type 1 AGN sample from the SINFONI survey for Unveiling the Physics and Effect of Radiative feedback¹ (SUPER² Circosta et al. 2018). The SUPER survey is designed to overcome the limitation of low spatial resolution in order to understand the true spatial extent of outflows in AGN host galaxies at high redshift. Apart from the extension of the outflows and the ionised gas, SUPER aims at answering some of the fundamental questions such as: how prevalent are ionised outflows in X-ray selected AGN host galaxies? Do the outflow properties such as velocity and mass outflow rates show a correlation with those of the AGN and of the host galaxy (e.g. AGN bolometric luminosity, host star formation rate)? And if they do, what is the scaling relation between the corresponding parameters? Lastly, do these outflows have any effect on the host (e.g. its star formation rate or gas mass)?

This paper is arranged as follows: in Sect. 2, we describe the sample used and its characteristics, and in Sect. 3 we discuss the observing strategy, followed by the data reduction procedure and stacking of data cubes in Sect. 4. We describe in detail the analysis procedure and the corresponding results obtained in Sect. 5, discuss the implications of these results and derive scaling relations in Sect. 6. Concluding remarks are reported in Sect. 7.

Throughout this paper, we adopt the following ΛCDM cosmology parameters: H_o = 70 km s⁻¹, Ω_M = 0.3, Ω_Λ = 0.7.

2. Sample description

The sample presented in this paper is derived from the SUPER survey, a large programme with the Spectrograph for INtegral Field Observations in the Near Infrared (SINFONI, Eisenhauer et al. 2003) mounted on the Cassegrain focus of Unit Telescope 4 (UT4) at the Very Large Telescope (VLT). We briefly discuss the overall properties of the SUPER sample in this section, in particular the Type 1 sub-sample used in this paper. For more details on sample selection and the survey characteristics, we refer the reader to Circosta et al. (2018).

The SUPER survey consists of a sample of 39 AGN selected in the X-rays (L_{2 − 10 kev} >  10⁴² erg s⁻¹) in the Chandra Deep Field South (e.g. Luo et al. 2017), COSMOS-Legacy (e.g. Civano et al. 2016), the wide area XMM-XXL (e.g. Georgakakis & Nandra 2011; Liu et al. 2016; Menzel et al. 2016), and Stripe 82 X-ray (e.g. LaMassa et al. 2016) surveys and from the WISE/SDSS selected Hyper-luminous quasar sample (e.g. Bischetti et al. 2017). The X-ray selection with the luminosity cut ensures a pure AGN selection as there is no contamination from the host galaxy and/or X-ray binaries at these energies (see e.g. Brandt & Alexander 2015; Padovani et al. 2017). The sample covers a redshift range of 2.1–2.5, which is the epoch of maximal activity of the volume averaged star formation in galaxies and the growth of black holes in the universe, making it ideal to study effects of radiative feedback from the black hole on the host galaxy (e.g. Madau & Dickinson 2014). Owing to the presence of ancillary multi-wavelength data sets, we are able to derive accurate measurements of the black hole and the host galaxy properties via spectral energy distribution fitting of UV-to-FIR photometry and X-ray spectral fitting (details in Circosta et al. 2018 for the analysis of the multi-wavelength data sets and (Vietri et al., in prep.) for the black hole mass estimations). Of the 39 targets, 22 were classified as Type 1 (56%), and the remaining 17 as Type 2 (44%), based on the presence or absence of broad emission lines such as MgII or CIV in the optical spectra. One target, cid_1205, which was previously reported as a Type 2 AGN in Circosta et al. (2018) is now classified as a Type 1 AGN based on the presence of broad line region (BLR) emission in Hα line in the K-band SINFONI spectrum. So, 23 galaxies (58%) from the SUPER sample are now classified as Type 1 AGN. The overall SUPER sample spans a wide range in AGN and host galaxy properties, which allows us to identify any existing correlation between the outflow properties derived from the SINFONI data and those derived from the multi-wavelength ancillary data set.

In this paper, we focus on the analysis of the H-band SINFONI data of the Type 1 targets from the SUPER survey. The selected sample in the context of other SUPER targets and on-going IFU surveys is shown as cyan star symbols in Fig. 1. These targets populate the moderate-to-high luminosity range and span ∼2.5 orders of magnitude in bolometric luminosity among the SUPER sample. We explore the following range of black hole and host galaxy properties for the spatially resolved data set presented in this paper: log M_*/M_⊙ ∼ 10.38–11.20, SFR < 94–686 M_⊙ yr⁻¹, log L_bol/[erg s⁻¹] ∼ 45.4–47.9, log M_BH/M_⊙ ∼ 8.3–10.7, and log N_H cm⁻² < 21.25–> 24.1. The properties of the individual targets obtained from the spectral energy distribution (SED) fits and the available optical data are reported in Table 1.

Fig. 1. Past and ongoing IFU surveys targeting the [O III] λ5007 emission line to study the properties of outflows in AGN host galaxies in the Lbol versus redshift plane. The grey data points represent surveys at z <  1 (Husemann et al. 2013, 2014, 2017; Liu et al. 2013; Harrison et al. 2014; Karouzos et al. 2016; Bae et al. 2017; Rupke et al. 2017), while the solid black symbols are collected from z >  1 studies (Nesvadba et al. 2006, 2017; Alexander et al. 2010; Harrison et al. 2012; Cresci et al. 2015; Carniani et al. 2015; Brusa et al. 2016; Kakkad et al. 2016; Vayner et al. 2017; Vietri et al. 2018). The blue shaded area covers the parameter space of the KASHz survey (Harrison et al. 2016) and the red shaded area that of the overall SUPER survey (Circosta et al. 2018). This paper presents an analysis on the Type 1 sub-sample of the SUPER survey, the data points of which are shown as cyan star symbols. Adapted from Circosta et al. (2018), to which we refer the reader for further details.

3. Observations

The SUPER SINFONI observations were carried out between November 2015 and December 2018 as part of the ESO large programme 196.A-0377 in both service mode and visitor mode. We used SINFONI for AO-assisted observations in the H-band (1.45–1.85 μm) to trace the rest-frame optical lines Hβ and [O III] λ5007, and the K-band (1.95–2.45 μm) to trace [N II], Hα and [S II] λ6716, 6731, all of which can be ionised by both AGN and star formation processes. Due to the lack of an appropriate natural guide star, most of the observations were carried out with the laser guide star in clear sky conditions. Hence, all but three targets were observed using the seeing enhancer mode, which offers improved image quality compared to the natural seeing, since no suitable bright tip-tilt star was available close to the chosen targets, given their location in deep fields. We used a plate scale of 3″ × 3″ with a spatial sampling of 0.05″ × 0.1″, which gets re-sampled to 0.05″ × 0.05″ in the final data cube. Three targets were observed in seeing-limited mode during the visitor-mode runs when the conditions were not ideal to close the AO loop. The plate scale for observations without AO is 8″ × 8″ with a spatial sampling of 0.25″ × 0.25″ in the final reduced data cube. The average spectral resolutions in the H-band and K-band are ∼3000 and ∼4000, respectively, which translates to a channel width of ∼2 Å and ∼2.5 Å, respectively.

Most of the targets in the SUPER survey are too faint for direct acquisition during the observation (23.49 <  H_mag <  15.72, 22.27 <  K_mag <  15.34), and therefore a blind offset from a nearby bright star was used. Before each science observation, a dedicated point spread function (PSF) star was observed to get an estimate of the image quality after the AO correction and compare it to the natural seeing at the observatory. The PSF observation lasted 30–60 s on the source along with a sky exposure of a similar duration. This PSF observation was, however, not strictly used to select or discard the observation blocks to be used for the final cubes, because occasionally the conditions varied during the one hour observation. In case the conditions degraded or improved significantly, the exposures were discarded or included accordingly. The AO-assisted observations gave an average image quality of up to ∼0.22″ (median = 0.3″) in H-band and ∼0.2″ (median = 0.3″) in K-band, inferred from the FWHM of the dedicated PSF star observation. The PSF of the observations for each band is reported in Table 2.

During science exposures, we used a dithering pattern where the target was moved within the SINFONI field of view (FoV) so the sky for a particular frame was obtained from the subsequent frame and vice-versa. In the case of extended targets, the sky subtraction following such a dithering pattern would lead to the subtraction of the object signal itself and therefore a dedicated sky exposure was taken in a pattern “O–S–O–O–S–O” (S = sky; O = Object). Each object and sky exposure was limited to ten minutes to minimise the variation of the infrared sky. The total on-source exposure time for the targets in either band (H or K) ranges between 1 hour and ∼6 hours. Observing patterns and exposure times in either band for each target are summarised in Table 2.

To correct for atmospheric absorption and to flux calibrate the final co-added science cube, telluric stars were observed with the same setup as the science observations within 0.2 air mass and 2 hours of the science observations. Each telluric star exposure lasted 2–3 s with number of integration of 5, along with a sky observation of similar exposure time as that of the star. The stars were selected to have a K-band magnitude between 7 and 8.5 and a stellar type of B2V, B3V, B4V, or B5V.

Out of the 23 Type 1 AGN in the SUPER sample, 21 were observed in both the H-band and K-band (Table 2), 18 were detected in the H-band, and all 21 were detected in the K-band. The AGN X_N_53_3 was detected neither in continuum nor in emission lines in the H-band, but detected in K-band. The AGN X_N_44_64, although detected in continuum in both bands, lacks emission lines in the H-band. Two targets (lid_206 and S82X2106) from the list of Circosta et al. (2018) could not be observed since part of the survey was executed in visitor mode rather than service mode as originally planned, and therefore we had a higher fraction of time lost for bad weather conditions.

4. Data reduction

We used the latest version of the ESO pipeline (3.1.1) to reduce the SINFONI data. The pipeline corrects for the presence of non-linear and hot pixels, flags the pixels that have flat lamp intensities higher than a given threshold, and performs a flat field correction. It also computes optical distortions and slitlet distances and performs wavelength calibration using exposures from the xenon+argon arc lamp in the H-band and the neon+argon arc lamp in the K-band. Science exposures, PSF, and telluric star observations were reduced using the recipe sinfo_rec_jitter, which outputs re-sampled data cubes of the individual exposures of the science frames as well as the PSF and telluric cubes corrected for the distortions, bad pixels and calibrated for wavelength using the above mentioned steps.

The sky subtraction was performed externally using the improved sky subtraction procedure described in Davies (2007). During this procedure, 10–30% of the pixels from an object-free region were used to create a model sky spectrum, which was shifted in wavelength space to match the wavelength axis of the object frames. The processed sky spectrum was then subtracted from the object exposure across the FoV.

To remove the telluric absorption features and to flux calibrate the data cubes, first the hydrogen features were removed from the observed telluric star spectrum, divided by a black body spectrum and normalised it to obtain the response function of the instrument. Both the science and the telluric cubes were divided by this response curve to correct for the telluric features. The spectrum extracted from the corrected telluric cube was then convoluted with the appropriate filter (H-band or K-band) from the 2MASS catalogue (Two Micron All-Sky Survey: Skrutskie et al. 2006) to get the required flux per unit count to be applied to the entire data cube.

Lastly, the flux-calibrated individual frames from multiple exposures were combined using the pipeline recipe sinfo_utl_cube_combine with a sigma clipping parameter (–ks_clip=TRUE), and by scaling the sky (using –scale_sky=TRUE) within individual exposures. It has been verified that the sigma clipping does not remove signal from the original target itself. By setting the sky scaling parameter, the spatial median of each exposure is subtracted from the contributing exposure to remove sky background that might not have been removed in the previous steps of the reduction. For observations within the same night, the offsets given by the header keywords CUMOFFSETX/Y are verified as reliable by comparing the stacked cube with manually calculated offsets. Observations taken on a different night might have shifts in the centroid of the image. We performed a 2D Gaussian fit to the combined cubes obtained from individual observing blocks during each night, and the difference in the centroid of the Gaussian fits gave the relative offsets between the observations from different nights. After determining these offsets, each contributing cube was aligned and co-added such that the intensity of a pixel is given by the weighted mean of the intensity of the corresponding overlapping pixels from the individual cubes, where the weight depends on the exposure time of the individual frames. Any residual cosmic ray signal within the final co-added cube is then removed using a sigma clipping procedure. The total on-source exposure time of each target in the H- and K-bands is reported in Table 2.

The WCS (world coordinate system) coordinates of the co-added cube resulting from the SINFONI pipeline are inaccurate. We therefore applied an astrometric correction registering the peak of the continuum emission from the AGN with the optical/near-infrared coordinates reported in Circosta et al. (2018). Furthermore, we noted that the SINFONI observations of target J1549+1245 show a systematic spatial shift in the centroid of the continuum location as a function of wavelength. This could be due to an inaccurate correction of the atmospheric dispersion and/or to a rotation in the grating between the time of the science observation and when the wavelength calibration was obtained. To correct for such spatial shifts, we measured the offset along the X and Y directions as a function of wavelength and the derived offset functions were used to re-align the cubes using the drizzle algorithm (Fruchter & Hook 2002).

5. Analysis and results

The overall analysis of the reduced SINFONI cubes consists of three steps: (1) derive the properties of the ionised gas by modelling the integrated spectrum; (2) determine if the [O III] emission is resolved after subtracting the AGN PSF and/or comparing the spatial profiles with a suitable PSF model; and (3) measure the extension of the ionised gas at various velocity slices to determine the extent of the outflows. For the purpose of this paper, we focused the analysis only on the H-band SINFONI data to trace the kinematics of the ionised gas using the forbidden [O III] transition. The results from the analysis of the K-band spectrum will be presented in a later publication.

5.1. Modelling the integrated spectrum

In this section, we describe the line fitting procedure on the integrated spectrum. While this paper focuses on the narrow line region (NLR) properties, a more in-depth discussion on the BLR properties will be presented in a forthcoming publication (Vietri et al., in prep.).

The integrated spectrum for each object was extracted from a circular aperture centred on the target, which includes at least ∼95% of the total emission. The aperture used to extract the spectrum in each object is reported in Table 3. The target centre was calculated using a 2D Gaussian fit on the H-band continuum image obtained by collapsing the cube over all the spectral channels. The error on the spectrum was estimated creating an rms spectrum obtained from an object-free region. The analysis of the H-band spectrum was restricted to the region spanned by Hβ, [O III] λ4959, and [O III] λ5007. The residual sky lines were masked from the spectrum during the fitting procedure.

We modelled the extracted spectra using the scipy.curve-fit package in python, which uses the principle of least squares to find the optimal set of parameters for a given fitting model. We used a simple linear model to fit the AGN continuum. The iron emission was modelled using observed FeII templates from the literature (Boroson & Green 1992; Véron-Cetty et al. 2004; Tsuzuki et al. 2006). The [O III] and Hβ emission lines were modelled with Gaussian functions, and the kinematic components of the two lines were coupled with each other. The number of Gaussian components for the [O III] emission line was restricted to two, and the addition of the second Gaussian depended on whether it minimised the reduced chi-square value of the overall model. We refer to these individual Gaussian components as “narrow” or “broad” according to the values of the line width (FWHM), which are left as free parameters in the fitting procedure. We do not associate a physical meaning to each single Gaussian component, and consequently to the terms “narrow” and “broad”, but we use a non-parametric approach in the paper to define velocity as described further down in this section. For Hβ, a third Gaussian component or a broken power law was required to reproduce the BLR emission³, which has been used to infer the black hole masses of the SUPER targets (more details in Vietri et al., in prep.). The line centroid and width of the narrow and broad components of [O III] λλ4959,5007 and Hβ are tied to each other, based on the assumption of a common origin for these emission lines. Furthermore, the emission line ratio [O III] λ5007:[O III] λ4959 is set equal to ∼3:1 based on theoretical values (e.g. Storey & Zeippen 2000; Dimitrijević et al. 2007). In order to estimate the uncertainty on the derived parameters, we created 100 mock spectra by adding rms noise to the modelled spectrum and repeated the line-fitting procedure on these mock spectra. The errors reported in Table 3 are the standard deviation for each parameter obtained with this procedure. Since we lack sufficient constraints to correct for dust reddening, we do not attempt to correct the emission line luminosities for extinction effects.

As an example of the fit performed on the integrated H-band SINFONI spectra, we show the object X_N_115_23 and the spectral model in Fig. 2. The integrated spectra for the rest of the Type 1 SUPER targets can be found in Appendix A. The line-fitting parameters of the H-band spectrum are reported in Table 3. Out of the 21 Type 1 AGN presented in Table 2, we were able to derive line properties for 19 objects. The AGN X_N_53_3 is not detected in continuum and it does not show any emission lines in the H-band either. Furthermore, X_N_44_64 is detected in the continuum but does not have any emission lines. Finally, J1441+0454 is well detected in Hβ, but the spectrum does not show the presence of [O III] emission at the expected observed wavelength. This object is one of the brightest AGN in our sample, and the possible lack of [O III] emission at high bolometric luminosities has been previously reported in the literature (e.g. WISSH quasars: Bischetti et al. 2017). On the other hand, while the line fitting suggests that the spectrum around the expected location of [O III] is dominated by strong iron emission, it also found a significant highly blueshifted (∼−3000 km s⁻¹) [O III] line with FWHM∼1956 km s⁻¹. Such extreme blueshift has also been observed before in the literature in high-redshift extremely red quasars (ERQs, see Perrotta et al. 2019). Due to the degeneracy between the iron component and the [O III] emission and the highly blended nature of the observed Hβ-[O III] emission (see Fig. A.3), we limited the analysis of the [O III] properties to the integrated spectum for this object and do not attempt to characterise the extended nature of the [O III] emission in the next sections.

Fig. 2. Integrated H-band spectrum of X_N_115_23 shown here as an example. The grey curve shows the observed spectrum, the red curve shows the reproduced overall emission line model, the blue dashed curve shows the iron emission, and the dashed green curves show the continuum emission and the individual Gaussian components (narrow, broad, and BLR) used to reproduce the profiles of various emission lines. The blue vertical lines indicate the location of Hβ, [O III] λ4959 and [O III] λ5007. The vertical grey regions mark the channels with strong sky lines that were masked during the fitting procedure. The X-axis shows the rest frame wavelength after correcting for the redshift of the target, and the Y-axis shows the observed flux density. The integrated spectra of the rest of the targets are shown in Appendix A.

We adopted non-parametric measures for the line properties, which has the advantage that the parameter values do not depend on the fitting function adopted (e.g. the number of Gaussian components), which may strongly depend on the signal to noise of the spectrum under investigation (see e.g. Zakamska & Greene 2014 and Harrison et al. 2014 for more details). In particular, we measure the velocity at the tenth percentile of the overall [O III] line profile (v₁₀) and the velocity width of the line that contains 80% of the line flux (w₈₀ = v₉₀ − v₁₀). For a Gaussian profile, the value of w₈₀ approximately corresponds to the FWHM of the emission line. We also computed the maximum velocity, v_max, following the definition by Rupke & Veilleux (2013), as the shift between the narrow and the broad Gaussian components of [O III] plus twice the sigma of the broad Gaussian. For fits with single Gaussian, we estimate v_max as twice the sigma of the broad Gaussian. The non-parametric measures of velocities and the line widths are reported in Table 3.

In the following, we use w₈₀ to identify AGN with clear signatures of outflows. Figure 3 shows the distribution of w₈₀ for the AGN sample from SUPER and KASHz surveys in red and black histograms, respectively. The KASHz sample is matched in redshift (z ∼ 2.0–2.5) to the SUPER sample. For comparison, we also plotted the w₈₀ distribution of low-redshift mass-matched sample of MANGA galaxies from Wylezalek et al. (2020), which is shown in blue in Fig. 3. The light blue vertical area denotes the w₈₀ value from a stacked spectrum of mass-matched and redshift-matched galaxies from the KLEVER survey (Curti et al. 2020), a high-redshift survey targeting star forming galaxies. The stellar mass of the overall SUPER sample has been used to mass-match the comparison samples, as most of the Type 1 AGN do not have reliable stellar mass measurements (Circosta et al. 2018).

Fig. 3. Top panel: distribution of the non-parametric velocity, w80 measured for [O III] λ5007 line, for AGN and star forming samples from different surveys. The AGN sample consists of targets from the SUPER survey in red (this paper) and targets from KASHz survey matched in redshift in hatched black (Harrison et al. 2016). The Type 1 sample used in this paper has a median hard X-ray luminosity of 1044.8 erg s−1 and the redshift-matched KASHz targets have a median X-ray luminosity of 1043.9 erg s−1. The blue histogram shows the w80 distribution of mass-matched low redshift star forming sample from Wylezalek et al. (2020), while the ocean-blue vertical bar shows the w80 value from the redshift and mass-matched star forming sample from the KLEVER survey (334 ± 35 km s−1, Curti et al. 2020). The dashed black-line at 600 km s−1 corresponds to the w80 value used in this work to define if a target hosts an AGN-driven outflow. Lower panel: inverse cumulative w80 distribution with the same colour coding as in the top panel. It is clear from both the panels that targets from the SUPER survey occupy the higher end of the w80 distribution. Based on the w80 criteria, all the targets in the SUPER survey show the presence of outflows, and ∼52% of the redshift matched targets in the KASHz survey show outflows.

From the histogram and the cumulative distributions in Fig. 3, it is clear that the AGN sample from the KASHz and SUPER surveys occupy the higher end of w₈₀ values compared to mass-matched samples at low as well as high redshift star forming galaxies. A w₈₀ value of 600 km s⁻¹ lies at the high end of the tail of the w₈₀ distribution for low-redshift star forming galaxies and it is well above the average values obtained from KLEVER. Consequentely, in this paper we consider a w₈₀ value greater than 600 km s⁻¹ as a signature of an AGN-driven outflow. We note that a similar cut in w₈₀ was also used in previous works (e.g. Harrison et al. 2016) and is a conservative estimate when compared to the cuts used in other works (e.g. 500 km s⁻¹ in Wylezalek et al. 2020).

Based on the above definition of outflows, all the observed Type 1 AGN from SUPER detected in [O III] show the presence of AGN-driven outflows (90% outflow detection rate if we include targets observed and not detected in [O III] ). The KASHz survey, on the other hand, reported ∼50% of their sample having w₈₀ larger than 600 km s⁻¹ (Harrison et al. 2016), therefore a lower fraction of AGN with clear outflow signatures according to the adopted definition. If we restrict the KASHz sample to the same redshift range of SUPER, 2 <  z <  2.5, the fraction of KASHz targets hosting outflows is at ∼52%, as shown in Fig. 3. The difference between the w₈₀ distributions for SUPER and KASHz surveys is probably due to the different luminosity range of the AGN sampled by these surveys. Figure 4 shows the relation between total [O III] luminosity and hard X-ray luminosity for the targets presented in this paper as well as for the redshift-matched AGN sample from the KASHz survey. We have eight objects in common between the two surveys (for both Type 1 and Type 2 targets) and the w₈₀ values derived from SINFONI and KMOS for these objects are perfectly consistent with each other. As can be clearly seen from this figure, there is little overlap between the luminosity ranges covered by the two surveys: ≈90% of the Type 1 sample presented in this paper has L_{[2 − 10 kev]} >  3 × 10⁴⁴ erg s⁻¹, while ≈90% of the KASHz AGN are at L_{[2 − 10 kev]} <  3 × 10⁴⁴ erg s⁻¹ (Fig. 3, blue data points). As already reported in previous works (e.g. Fiore et al. 2017) and presented for the SUPER Type 1 sample in Sect. 6, there is a positive correlation between the velocity associated with the outflow and the bolometric luminosity of the AGN. Therefore, the detection of a higher fraction of outflows in the SUPER sample compared to the KASHz survey could naturally be explained with the prevalence of the faster outflows at higher bolometric luminosity (see Fig. 3).

Fig. 4. Total [O III] luminosity versus hard X-ray luminosity (LX(2–10 keV) for the SUPER (red circles) and KASHz targets (blue data points), the latter from a similar survey targeting X-ray selected AGN at high redshift. The combined data set of KASHz and SUPER show a positive correlation between the two quantities, shown by the black line. The shaded region represents the 1σ scatter on the relation. KASHz targets occupy the lower X-ray luminosity range compared to SUPER targets, which may explain the difference in the distribution of w80 in the two surveys in Fig. 3.

Although we use non-parametric measures to characterise the line properties, we also report here the incidence of outflows in our sample based on the presence of a broad Gaussian component in the line fit. For the 21 Type 1 AGN presented in this paper, 14 targets require a broad component with widths ∼1000–3000 km s⁻¹ in their [O III] profile. We can therefore claim based on this measure that 67% (∼75% if we limit to the objects with detected [O III] ) of the SUPER Type 1 sample has fast ionised outflows. This outflow incidence fraction is much higher than some previous works for star forming as well as AGN host galaxies (e.g. Harrison et al. 2016; Förster Schreiber et al. 2018; Swinbank et al. 2019) and comparable to other lower redshift studies targeting AGN (e.g. Rakshit & Woo 2018; Davies et al. 2020b). As we mentioned earlier, this method of classification for outflows is highly dependent on the signal to noise of the spectra and models used for line fitting, and therefore these comparisons are subject to biases.

Apart from the detection of the bright emission lines already discussed in this section, we also detect Hβ in all Type 1 sources, except cid_1205, which shows a strong skyline at the location of Hβ emission. All targets but X_N_66_23 show the detection of a BLR Hβ component. The non-detection of the BLR component of Hβ in X_N_66_23 is possibly a consequence of the short exposure time leading to a low signal to noise (S/N) in the spectrum. Compared to the [O III] profile where the broad component was detected in ∼75% of the sample, 62% of the sample required an additional broad Gaussian component for Hβ line. The Hβ component flux is usually ten times fainter than that of [O III] line in Type 1 AGN (e.g. Leighly 1999; Rodríguez-Ardila et al. 2000), hence for the current sample the non-detection of broad Hβ could be simply due to the lower S/N. Although the kinematic components of the [O III] and the Hβ line are coupled to each other, the relatively low S/N of the Hβ line means that the addition of second Gaussian in Hβ profile does not change the χ² of the fit in some targets. Based on this consideration, we will use the results on the [O III] line profile to assess the incidence of AGN-driven outflows in our sample. We refer the reader to Vietri et al. (in prep.) for further discussion on Hβ line profile.

5.2. Extension of the [O III] emission line region

Before discussing the kinematics of the ionised gas, we now focus on the actual extension and morphology of the emission line region. First, we need to verify whether the [O III] emission is actually extended or unresolved, and where what may erroneously be interpreted as extended emission is nothing more than beam smearing effects (e.g. Husemann et al. 2016; Villar-Martín et al. 2016). We used two techniques to asses if the emission is truly extended. The first consists of comparing the curve of growth (COG) of the total [O III] emission of the AGN with that of an observed and modelled PSF. The second technique, which we call the “PSF-subtraction” method, consists of producing a residual [O III] map after subtracting the nuclear [O III] spectrum spaxel by spaxel following a 2D PSF profile.

To construct the curves of growth, we started by extracting the spectrum from an aperture of radius 0.1″ centred at the peak location of the AGN continuum emission and perform the line fitting as described in Sect. 5.1. From the best fit, we derived the flux value for the total [O III] emission, and repeat this procedure increasing the extraction radius in steps of 0.05″. Finally, we reconstructed the curve of growth by plotting the line fluxes derived as a function of the aperture radius. We only plotted the total [O III] flux and not the individual components in these plots as we do not give any physical significance to individual Gaussian components. The COG are derived for 11 targets (out of the 21 in the Type 1 sample observed) for which the S/N in the integrated spectrum extracted from their respective apertures (Table 3) is greater than 5. The H-band data cubes of cid_1205 and S82X2058 are contaminated by a bright stripe, due to which spatially resolved analysis of these targets is not performed.

We then need to have an accurate description of the PSF to compare with the [O III] emission. We had designed our SINFONI observations to have dedicated PSF star observations close in time and space for each single science observations, and we can therefore repeat the COG procedure (described above for the AGN) for the dedicated PSF star, in order to derive the PSF profile as a function of the distance from the centre. We were not able to construct the growth curve of PSF observation for J1333+1649 and J1549+1245 as the PSF star was at the edge of the data cube. Nevertheless, since the AGN in this paper are Type 1, we have access to the spatially unresolved BLR emission as traced by Hβ, which we can use as an alternative method to trace the PSF profile (the only exception being X_N_66_23, for which we do not have a clear detection of the Hβ line). Out of the eight objects for which we can use both methods to trace the PSF, both the dedicated PSF star and the Hβ BLR component give consistent PSF profiles for five objects, so we were able to use both to trace the spatial resolution of the data cubes. For the remaining three objects (X_N_115_23, cid_346, S82X1905), the PSF profile as traced by the dedicated star observation is narrower than those derived from the Hβ profile. These differences between the PSF traced by the BLR emission and the PSF star could be due to the difference in the AO correction, which might change in long exposure observations. Also, as the PSF star observations were performed at the beginning of the ∼1 hour-long science OBs, for these three objects the conditions got relatively worse during the science observations. Therefore, we considered the PSF profile derived from the BLR to be the correct one to compare since they trace the conditions simultaneously with the science observations. Finally, we estimated errors on the COGs by repeating the fitting procedure 100 times on mock spectra obtained by adding rms errors derived from an object-free location of the spectrum on the original models. The COGs for the 11 targets with an S/N >  5 in each aperture are presented in Fig. 5.

Fig. 5. Curves of growth for the Type 1 sample as described in Sect. 5. The solid black line shows the curve from the image obtained by collapsing the PSF-observation cube at the location of [O III] λ5007 channels. The Hβ BLR component is shown as a dashed black curve with the grey area showing the 1σ uncertainly. The BLR component is unresolved and serves as a proxy to determine the PSF during the observations which is also compared with the dedicated PSF star observation. The red curve shows the curve of growth for the total [O III] emission, meaning narrow and broad components. All the uncertainties are at 1σ levels. The BLR component is not plotted for X_N_66_23 as it remains undetected. For the rest of the Type 1 sample in the SUPER survey, the S/N is not enough to perform such a curve-of-growth analysis.

We can now compare the COG of each single target with the PSF. As discussed above, whenever possible we considered the PSF curve obtained from the unresolved BLR the best representations of the observing conditions during the science observations. From Fig. 5, we find that for seven targets (≈63% of the targets for which this analysis was performed) there is a 2σ excess in the total [O III] emission compared to the PSF at the maximum radius sampled by the COG, and therefore we claim that in these objects the [O III] is spatially resolved. As was done also in previous works, we were able to compare the half-light radii (the radius containing 50% of the flux) of the [O III] and the PSF emission reported in Table 4 to asses the spatially resolved nature of the [O III] emission. The half-light radii were estimated using a linear spline to interpolate between the data points of the COG. For five objects (X_N_66_23, X_N_115_23, cid_346, J1549+1245, and S82X1905), the [O III] half-light radius is larger than the PSF half-light radius within ≈2σ. For the remaining targets, the COG of the [O III] emission and the half light radius are consistent with the PSF within 2σ. In particular, for the two targets X_N_81_44 and J1333+1649, the half-light radius method may not be sensitive to the fainter very extended emission beyond the PSF, contrary to the overall COG and PSF method described below.

As mentioned at the beginning of this section, we used a second method to asses the spatial extension of the ionised gas emission, namely the PSF-subtraction method described in Carniani et al. (2015). The basic principle behind the PSF-subtraction method is that if the observed emission is unresolved and/or is a result of beam smearing from the PSF of the AGN, the spectrum at any given distance from the location of the AGN is the same as the spectrum at the location of the AGN, except for a scaling factor across the spectrum (e.g. Jahnke et al. 2004). In this case, the BLR component of Hβ smears the emission across the FoV and is therefore AGN dominated. If the emission does not mimic the spectrum at the location of AGN, then the observed spectral line is resolved. Following this principle for our sample, we first modelled the spectrum extracted at the AGN location (within a radius of 0.1″), which was determined using the centre of the continuum in the H-band cube. We refer to this spectrum as the “nuclear model”, which we then subtracted from every pixel across the SINFONI FoV, after allowing a variation in the overall normalisation factor of the spectrum while keeping the rest of the kinematic components fixed. This is followed by creation of channel maps at the location of [O III] λ5007 emission in the nuclear-spectrum-subtracted cube. The channels used to create the residual maps were 6 Å wide, except cid_346, for which a channel width of 10 Å was used. Following the principle described above, a noisy residual map would indicate an unresolved emission, while a map with non-zero residual structures would indicate a resolved emission. Similarly to the case of the COG, we restricted this analysis to targets with an S/N larger than 5 in the integrated spectrum (extracted from the respective apertures reported in Table 3). We also restricted the analysis to targets that show BLR Hβ emission, as it is used to scale the overall normalisation of the nuclear spectrum. Since we did not detect a BLR component in X_N_66_23, we performed the PSF-subtraction analysis for ten SUPER targets. The results of this method are shown in Fig. 6. The left panels in these figures show the residual maps obtained after collapsing the nuclear-spectrum-subtracted cubes along the [O III] channels, while the right panels show the residual spectrum extracted at the region of excess emission. From Fig. 6, we see excess residuals at a level of > 2σ for six SUPER targets: X_N_81_44, X_N_115_23, cid_346, J1333+1649, J1549+1245, and S82X1905.

Fig. 6. Residuals obtained using “PSF-subtraction” method across the SINFONI FoV as described in Sect. 5 for the Type 1 sample. Left panel: residual image obtained by collapsing the channels at the location of [O III] λ5007 after subtracting the emission from the AGN. The colour map is on a linear scale with the yellow regions showing excess emission (i.e. > 3σ). North is up and east is to left. Right panel: spectrum extracted at the location of excess emission, with an aperture of 0.1″ × 0.1″ shown by the red square in the left panel. The vertical grey lines in the residual spectrum shows the location of the [O III] λ4959 and [O III] λ5007, respectively, while the channels used to create the residual map on the left are indicated with a magenta box. Non-zero residuals in the images show that [O III] emission is resolved, while a noisy pattern is indicative of unresolved emission, which are, respectively, marked in the left panels. The AGN X_N_66_23 is not included in this analysis due to non-detection of the BLR component.

In summary, from the joint analysis based on the COG and PSF-subtraction method, we identify seven sources from our sample that show extended [O III] emission. Further analysis on the spectra across the SINFONI FoV is therefore restricted to the seven targets, which show extended gas emission with at least one of the methods described above.

5.3. Modelling ionised gas kinematics across the FoV

In the following section, we describe the pixel-by-pixel analysis for the H-band raw data cubes of the seven targets resolved from the methods described in Sect. 5.2. For the modelling of the spectrum across every pixel, we used the parameters obtained in the integrated spectrum as a prior. Apart from the constraints set while modelling the integrated spectrum, we also fixed the centroid and the width of the unresolved BLR component of Hβ and only allowed a variation in its peak. In the H-band, the Gaussian parameters of other emission lines are allowed to vary. The emission line widths (FWHM) were also kept greater than ∼100 km/s to avoid any spurious fit to a sky line. The validity of the line fitting across each spaxel was checked by subtracting the emission line model from the raw data and making a collapsed map along the channels spanned by the [O III] λ5007 emission. Due to the low S/N in each spaxel in the cases of X_N_66_23 and cid_346, nearby spaxels were averaged within a radius of 0.1″ to produce binned spectra. This results in a higher S/N of the spaxels without compromising the respective resolution of the observations.

As a result of the pixel-by-pixel analysis, various line and continuum maps are created, which we describe below. These maps are shown in Figs. 7 and 8. All the maps are smoothed to the spatial resolution of the observations derived from the PSF and an S/N cut of 2 is employed in the [O III] profile. First, we created broad-band continuum maps, masking the emission lines, and used a 2D Gaussian fit to determine the photo-centroid of the resulting image, which we used to mark the AGN position.

Fig. 7. Flux and velocity maps of targets showing resolved/extended emission inferred from the COG and PSF-subtraction method (from top to bottom: X_N_81_44, X_N_66_23, X_N_115_23, and cid_346). Column 1: background image shows the flux map of total [O III] emission (narrow+broad). Column 2: [O III] line centroid map. Column 3: v10 map. Column 4: w80 map. In all the maps, the black star denotes the centre position of the H-band continuum, as a proxy for AGN location. North is up and east is left.

Fig. 8. Same as Fig. 7 for J1333+1649, J1549+1245, and S82X1905.

The total [O III] line flux maps are shown in the first column of Figs. 7 and 8, while the line centroid maps obtained from the peak of the total [O III] profile are presented in the second column. A smooth gradient is observed in some of the line centroid maps, for example in X_N_81_44, X_N_66_23, cid_346, and S82X1905, which could be indicative of a rotating disc in the host galaxy or a bipolar ionised outflow. The observed extreme high velocities in the w₈₀ maps (∼1000 km s⁻¹) for most of these galaxies suggest that the kinematics of the ionised gas are also likely affected by the central AGN. For several of our targets, we do detect extended features in emission, which are not produced by wings of the PSF as demonstrated by the analysis performed in the previous sections (e.g. X_N_115_23, CID_346, J1333+1649, J1549+1245, S82X1905). A possible explanation for these irregular features in the gas kinematics could be the presence of a merger event between two galaxies. As far as a major merger is concerned, we do not see clear evidence of multiple galaxies in the available ancillary data. We would, therefore, favour alternative scenarios in which the extended features could be interpreted as an inflow or outflow. Inflows are usually detected in absorption and are predicted to show lower velocities than those we measured for the SUPER sample (e.g. Bouché et al. 2013). Also, since models predict a low covering fraction for the inflows (e.g. Steidel et al. 2010), we are inclined to interpret these features as outflows. In order to further investigate the kinematic features of the ionised gas in these AGN, we also generated the maps with the spatial variation of v₁₀ and w_80, which are shown in the third and the fourth columns in Figs. 7 and 8. The inspection of both maps reinforces our conclusion that most of our objects do not present the kinematic signatures of an undisturbed rotating disc: the v₁₀ map shows regions with extremely blueshifted velocities of more than −1000 km s⁻¹, with a w₈₀ that peaks at the location of the maximum velocity shift. Although a proper modelling of the velocity field of these galaxies is beyond the scope of this paper, we discuss here the morphology and size of the outflows’ signatures.

Following the discussion of the integrated spectra in Sect. 5.1, we use the cut w₈₀> 600 km s⁻¹ to identify regions dominated by outflowing ionised gas. Projected velocities (e.g. v₁₀ maps in the third column in Figs. 7 and 8) are more sensitive to the galaxy inclination, while the w₈₀ maps are more likely to be less sensitive to inclination (e.g. see discussion in Harrison et al. 2012). For each source, we measured the maximum projected spatial extent of the w₈₀> 600 km s⁻¹ component of the [O III] line profile (we name this quantity D₆₀₀, following Harrison et al. 2014). From the D₆₀₀ values in the w₈₀ maps, we find that the outflows are extended from ∼2 kpc up to 6 kpc (see Table 4).

From the w₈₀ maps in Figs. 7 and 8, it is clear that the ionised outflows detected in these AGN present a diversity of projected structures. In five sources (X_N_81_44, X_N_115_23, cid_346, S82X1905) there is a velocity gradient along a particular axis, while for two objects (J1333+1649 and J1549+1245) the maps are more spherically symmetric. Some of the kinematic structures present in the w₈₀ maps could be consistent with a biconical outflow, where the inclination of the outflow with respect to the disc of the galaxy and the line of sight of the observer could explain the observed velocity maps. This means a high inclination of the outflow with respect to the line of sight results in a clear velocity gradient along the outflow axis, with possibly both the blueshifted and redshifted part of the outflowing gas detected and a low inclination of the outflow with respect to the line of sight would produce a more uniform and spherically symmetric velocity field. A proper kinematic modelling would be required to constrain the ionised gas kinematics, and we will present this analysis in an upcoming paper on the full survey sample.

5.4. Radial distance of high-velocity ionised gas using spectroastrometry

The analysis presented in Sect. 5.2 is useful in determining the presence or absence of extended ionised gas. However, the methods described there do not determine the bulk velocity of the observed extended gas, that is, if the extended component of the gas is in outflow. Conventionally, the velocity maps derived in Sect. 5.3 have been used to quantify the extension of the outflowing gas. These maps can only be used for targets that are well extended beyond the width (FWHM) of the observed PSF. For marginally resolved or unresolved targets, the velocity maps can be affected by PSF smearing. For the latter case, the spectroastrometry technique (Carniani et al. 2015) is useful to determine the radial distance of the bulk⁴ of the gas moving at certain velocity from the AGN location.

The principle of the spectroastrometry technique is that if the gas is moving at a distance R from the AGN position at a velocity different to the systemic velocity of the host galaxy (systemic here is defined by the position of the peak of the [O III] line in the spectrum), then the photo-centroid of the line at that velocity is also shifted by a distance R from the AGN position. In this method, we first collapse the original data cube along channels containing AGN continuum and the photo-centroid of the resulting image is determined with a 2D Gaussian fit, similar to the method described in Sect. 5.3. The photo-centroid of the continuum is used as an indicator of the AGN location. We then used the spectrum models derived in Sect. 5.3 to subtract all the components from the raw spectrum across every pixel, except those of the [O III] λ5007. The resulting cube containing only the [O III] λ5007 emission is then collapsed along wavelength channels with [O III] emission in bins of ∼200 km s⁻¹ such that we achieve an S/N of ∼2 in the [O III] flux for each spectral channel. Similarly to the continuum image, we then calculated the photo-centroid of each of the images obtained from the collapsed channels using a 2D Gaussian fit. The velocity bins of ∼200 km s⁻¹ ensures there are minimal errors in determining the centroid due to the spectral resolution of the instrument. The distance between the photo-centroid in each velocity channel and the AGN location was then plotted as a function of the centre velocity of the channels. We would like to point out that spectroastrometry does not calculate the size of the outflowing gas but the spatial offsets between ionized gas clouds at different velocities with respect to the systemic velocity of the host galaxy. The offsets can be measured at scales smaller than the spatial resolution of the observations, and therefore the method can provide the distance of the bulk of the gas moving at a given velocity in unresolved or marginally resolved targets.

We apply the spectroastrometry method to all 11 of the SUPER Type 1 AGN that have S/N >  5 in the [O III] line of the integrated spectrum. The results from the spectroastrometric analysis are shown in Fig. 9. The blue data points show the measured offsets from the AGN location, marked by the vertical dashed line at R = 0. The errors shown in the plot are the 1σ uncertainty obtained from the output of the two-dimensional Gaussian fitting procedure. For each target, we also produced the corresponding maps for each velocity channel included in the spectroastrometry analysis. Figure 10 shows the channel maps for X_N_115_23 as an example, while the channel maps for the rest of the targets have been moved to Appendix A.

Fig. 9. Results of spectroastrometry applied on the Type 1 SUPER targets. The X-axis represents the distance between the continuum centre (vertical dashed line) and the position of the photo-centre of the [O III] line at the corresponding velocity channel given in the Y-axis. The error bars represent 1σ uncertainty on the values. The red arrow marks the cloud with velocity greater than 600 km s−1 (300 km s−1 in the [O III] line profile) whose centroid is at the maximum distance from the AGN.

Fig. 10. 1.5″  × 1.5″ [O III] λ5007 channel maps of X_N_115_23 at different velocity slices, after subtracting the Hβ, [O III] λ4959 and iron models from the raw cube. Each velocity slice is 200 km s−1 wide and the displayed value is the centre velocity of the respective channel. Black cross marks the location of the H-band continuum peak, used as a proxy for the AGN position, and the cyan star marks the location of the [O III] centroid determined from a 2D Gaussian fit. North is up and east is to the left.

As presented in Sect. 5.1, we used a cut of w₈₀ >  600 km s⁻¹ to identify AGN with ionised outflows, and in the previous section we used the w₈₀ maps (for the object with spatially resolved emission) to characterise the maximum extension of the outflowing gas (D₆₀₀), similarly to previous studies (e.g. Harrison et al. 2014; Cresci et al. 2015). We can now use the spectroastrometry analysis to characterise (also for the objects that are not spatially resolved) at what distance the bulk of the outflowing gas is located. We marked the velocity channel with |v|> 300 km s⁻¹ with the maximum radial distance from the AGN location using a red arrow in the plots in Fig. 9. The corresponding velocity (v_o) and distance (R_o) are reported in Table 4. We inferred a maximum spatial offset of 0.30″ and a median offset of ∼0.1″ between the AGN location and bulk of the outflowing gas. This translates into a maximum physical physical distance of ∼2.2 kpc (0.8 kpc median), which is smaller than the spatial resolution of the observations in the H-band. These numbers are consistent with Carniani et al. (2015), where the maximum offset of outflowing gas is ≲2 kpc. We note that there could be multiple gas clouds at velocities above 600 km s⁻¹. We only highlight the ones at the maximum distance from the AGN. Therefore, the bulk of the high-velocity outflowing gas (|v|> 300 km s⁻¹) is generally concentrated within ∼2 kpc of the AGN location in the majority of the high-redshift AGN host galaxies. We would like to point out that the actual extent of the outflowing gas might be larger than the values obtained from the spectroastrometry, which aims to identify the distance from the bulk of the gas moving at a given velocity.

There are a couple of targets with notably interesting features from the spectroastrometry analysis. In the case of cid_346, the spectroastrometry method reveals the presence of gas moving at ∼600 km s⁻¹ at a distance of ∼0.3″ equivalent to a physical distance of ∼4 kpc from the AGN location. Similarly, we also find that the bulk of the gas moving at ∼300 km s⁻¹ is located at ∼4 kpc from the centre in the case of S82X1905. We argue that a possible explanation for such extended redshifted emission is that the receding part of the outflow is not obscured by the dust of the host galaxy in these two cases, possibly due to the larger extent of the outflow itself.

We also note from the spectroastrometry maps that the bulk of the low-velocity [O III] component, that is, the non-outflowing component of the [O III] emission is not necessarily emitted at the same location as that of the AGN in a few targets. This is expected, as even in low-redshift galaxies the [O III] emission is dominated in the extended NLR in the form of ionisation cones, which are not co-spatial with the AGN location.

We now summarise the different methods used to determine the extension of the ionised gas and the outflows associated with the gas. Using COG and PSF-subtraction methods, we compared the spatial distribution of the ionised gas with the observed PSF and inferred whether the ionised gas is extended or not. The COG and PSF-subtraction methods rely on robust PSF measurements either from the BLR emission, or, if the BLR emission is not available, a dedicated PSF observation. Using either of these methods, we find that the ionised gas is extended in X_N_81_44, X_N_66_23, X_N_115_23, cid_346, J1333+1649, J1549+1245, and S82X1905. In these targets with extended ionised gas, we constructed the flux and velocity maps and calculated the D₆₀₀ value, which is the maximum projected spatial extent of w₈₀ > 600 km s⁻¹. The D₆₀₀ parameter quantifies the spatial extent of the outflow associated with the ionised gas, which is consistent with the outflow definition we adopt throughout this paper, that is, w₈₀ > 600 km s⁻¹. The D₆₀₀ parameter is in the range of ≈1.5–6.5 kpc from the AGN location. The D₆₀₀ parameter should only be determined for targets that show the presence of extended ionised gas emission beyond the observed PSF. Lastly, we determined the radial distance at which the bulk of the ionised gas (weighted by luminosity) is moving at a velocity greater than 600 km s⁻¹ using the spectroastrometry technique. We find that for most galaxies, the bulk of the high-velocity gas (|v| > 600 km s⁻¹) is contained within ∼1 kpc of the AGN location. The spectroastrometry technique has the advantage that it can be used for targets that are unresolved or marginally resolved. However, the technique does not provide the true spatial extent of the ionised gas, but the centroid of the bulk gas moving at a specific velocity.

5.5. Mass outflow rates

To understand the impact that AGN feedback may have on their host galaxies, it is important to estimate the mass of the ouflowing gas and the kinetic energy associated with it. Although the ionised gas phase is believed to trace only a fraction of the total outflow mass (e.g. Cicone et al. 2018), the energy associated with these outflows can provide valuable information on its driving mechanism (e.g. Fiore et al. 2017; Brusa et al. 2018; Husemann et al. 2019; Shimizu et al. 2019).

An accurate estimation of the mass outflow rate for high-redshift galaxies is a major challenge due to the limited spatial resolution achievable with current instrumentation and the uncertainties on the physical properties of the outflowing gas. As a consequence, one must make a number of assumptions, which leads to large uncertainties in the derived mass outflow rate. The key elements that need to be constrained to estimate the mass outflow rate are: the geometry and the volume of the outflow; the physical properties of the gas (e.g. its electron density, temperature, and associated metallicity); the velocity of the outflowing gas and its radial evolution. For the scope of this paper, we focus on reporting the range of mass outflow rates derived from our measured quantities by exploring the possible assumptions for key physical quantities, described as follows.

Regarding geometry, both theoretical models (e.g. Ishibashi et al. 2019) and results from imaging and IFU observations of low redshift galaxies (e.g. Crenshaw et al. 2010; Müller-Sánchez et al. 2011; Liu et al. 2013; Wylezalek et al. 2016; Venturi et al. 2018) provide evidence of the presence of spherical and conical outflow morphology. Recent observations of low redshift AGN host galaxies have also suggested the possibility of expanding shell-like shocks driving an outflow across the host galaxy (e.g. Husemann et al. 2019). A shell-like expanding outflow model can take into account variable physical quantities such as the outflow luminosity, velocity, density, and temperature as a function of distance from the AGN. The total outflow rate throughout the galaxy is thus the sum over all the shells. As current observations with SINFONI do not allow us to infer the outflow morphology, we calculated mass outflow rates assuming both biconical and thin-shell geometry. A spherical geometry would change the overall normalisation factor in the mass outflow rate formula, but not the general trends, when compared with the biconical outflow geometry. We considered the volume up to a distance of 2 kpc from the AGN (see further discussion on radius and velocity later in the section). For the thin-shell model, we assumed a shell width of ∼500 pc, which is the at upper end of the limits inferred in low-redshift AGN host galaxies (e.g. Husemann et al. 2019).

In terms of luminosity, for an emission line modelled using multiple Gaussian components, the broad Gaussian component is often associated with the outflowing gas. However, as stated in earlier sections, the detection of an additional broad component depends on the S/N of the data and the model function used. Therefore, we continued to use the non-parametric approach and define the outflow luminosity as the luminosity calculated from the flux of the [O III] λ5007 emission line channels at |v| > 300 km s⁻¹. This definition is consistent with our outflow’s classification, that is, w₈₀ > 600 km s⁻¹, presented in Sect. 5.1. The outflow luminosity as defined above is reported in Table 4. We note that the average ratio between the broad component luminosity (Table 1) and the outflow luminosity ranges from 0.56–4.07 with a mean ratio of 1.80 ± 0.89. Therefore, on average the luminosity of the broad Gaussian component is higher than that estimated from the non-parametric procedure described above. We finally note that the luminosity calculated in this paper is not corrected for extinction and therefore represents a lower limit on the outflow luminosity.

Regarding radius and velocity, several spatially resolved observations of ionised gas at low redshift clearly show a radial dependence of outflow properties (e.g. Kakkad et al. 2018; Venturi et al. 2018; Husemann et al. 2019). It is therefore necessary to take into account the variations in the velocity across the outflowing medium while calculating the mass outflow rates. Due to the relatively low spatial resolution of high-redshift targets compared to the low-redshift AGN, deriving spatial variations in the physical properties of the outflow is not always feasible. Observations at low redshift report outflow radius values up to few kiloparsecs (e.g. Liu et al. 2013; Husemann et al. 2016). For the purpose of this paper, we calculated the mass outflow rate up to a fixed radius of ∼2 kpc, which is also the maximum radius for unresolved targets (see PSF values in Table 2). The outflow velocity at 2 kpc will be derived from the w₈₀ maps in case of resolved targets (Figs. 7 and 8) and the integrated w₈₀ value in case of unresolved targets. Other definitions for outflow velocity have also been used in the literature, such as broad Gaussian widths and v_max, which can introduce a systematic uncertainty up to a factor of ∼2 in the calculation of mass outflow rate (e.g. Kakkad et al. 2016; Fiore et al. 2017).

The electron density represents a major source of uncertainty in the calculation of mass outflow rates (e.g. Harrison et al. 2018), since it is very challenging to measure it, and in fact there is a wide range of values reported in the literature for the gas density in the NLR (Perna et al. 2017; Kakkad et al. 2018; Baron & Netzer 2019; Davies et al. 2020b). Similarly to the outflow’s velocity, variations in the electron density across the outflowing medium make the calculation of outflow rates even more challenging. Moreover, key diagnostic emission lines to estimate electron density such as [S II] λ6716, 6731 are too faint or remain undetected for the SUPER targets. Therefore, instead of assuming a single value for the gas density, we considered two values that cover the wide range observed in the literature, that is, from 500 cm⁻³–10 000 cm⁻³, which would return maximum and minimum possible outflow rates given the rest of the model assumptions described above.

Since Hβ is relatively faint, we calculated the mass of the outflowing ionised gas from the [O III] λ5007 line as follows (see e.g. Carniani et al. 2015; Kakkad et al. 2016):

$$\begin{array}{cc}& M_{[\text{O}\text{III}]}^{\mathrm{out}}=0.8\times {10}^8{M}_{\odot }\left(\frac{1}{{10}^{[\mathrm{O}/\mathrm{H}]-{[\mathrm{O}/\mathrm{H}]}_{\odot }}}\right)\left(\frac{L_{[\text{O}\text{III}]}}{{10}^{44}{\mathrm{erg}\mathrm{s}}^{-1}}\right)\hfill \\ \hfill & {\left(\frac{<n_e>}{500{\mathrm{cm}}^{-3}}\right)}^{-1},\hfill \end{array}$$(1)

where all the oxygen is assumed to be ionised to O²⁺, the electron temperature is 10 000 K, and [O/H], the metallicity of the outflowing material, is assumed solar. For a uniformly filled biconical outflow, the mass outflow rate takes the form (Fiore et al. 2017):

$$\begin{array}{c}\hfill {\dot{M}}_{[\text{O}\text{III}],\mathrm{cone}}^{\mathrm{out}}=3·v_{\mathrm{out}}\frac{M_{\mathrm{out}}}{R_{\mathrm{out}}},\end{array}$$(2)

where v_out is the outflow velocity, M_out is obtained from Eq. (1), and R_out is the distance at which the outflow rate is computed. We are making the assumption here that the biconical outflow region is uniformly filled with ionised gas with density ⟨n_e⟩. Equation (2) gives the instantaneous outflow rate of the ionised gas crossing the spherical sector at distance R_out from the AGN. On the other hand, if the outflow is assumed to propagate in thin shells of thickness ΔR (assumed 500 pc here), which cause a shock against the ISM, the outflow rate is simply (Husemann et al. 2019):

$$\begin{array}{c}\hfill {\dot{M}}_{[\text{O}\text{III}],\mathrm{thin}\text{-}\mathrm{shell}}^{\mathrm{out}}=v_{\mathrm{out}}\frac{M_{\mathrm{out}}}{\mathrm{\Delta }R},\end{array}$$(3)

where v_out is the outflow velocity at the location of the shell (assumed to be at 2 kpc here). Equation (3) provides the outflow rate within the shell of thickness ΔR instead of an integrated value.

We computed the mass outflow rates using both formulas described above and for two values of the electron density: 500 cm⁻³ and 10⁴ cm⁻³ (e.g. Baron & Netzer 2019; Davies et al. 2020b). The values obtained are reported in Table 4. We note that most of the SUPER targets have mass outflow rate values below 10 M_⊙ yr⁻¹ for the density range 500 cm⁻³–10⁴ cm⁻³. At the high-luminosity end, however, the mass outflow rate values exceed 10 M_⊙/yr even for a conservative electron density value of 10⁴ cm⁻³. In Fig. 11, we show the range of mass outflow rates versus the bolometric luminosities covered by our sample, from literature ionised gas outflow rates compiled in Fiore et al. (2017), and from the low-redshift X-ray AGN studied in Davies et al. (2020b). The mass outflow rates values from Fiore et al. (2017) and Davies et al. (2020b) were updated to match the electron density range used in this study. In addition, these literature mass outflow rates from Fiore et al. (2017) and Davies et al. (2020b) were also re-scaled to match the methodology used in this paper. Fiore et al. (2017) used the luminosity of the broad Gaussian component as a proxy for the outflow luminosity, while, as described above, we only considered channels with |v| > 300 km s⁻¹. For our objects, we verified that the mean ratio between the luminosity calculated from the broad Gaussian and the luminosity from [O III] channels at |v| > 300 km s⁻¹ is ≈2 and re-scaled the values from Fiore et al. (2017) by this factor. A similar scaling was used by Davies et al. (2020b), who used the [O III] luminosity to calculate the outflow mass. The Ṁ_out − L_bol relation in Fiore et al. (2017) was also used, when available, with either Hα or Hβ to calculate the total mass outflow rates for the ionised gas. Outflow rates for sources with only [O III] detection were scaled by a factor of three, which was consistent with the mean ratio between the Hα and [O III] luminosities in their sample. We therefore applied the same factor to their relation to make it consistent with us using the [O III] luminosities. After taking these factors into account, the Fiore et al. (2017) relation roughly matched the mass outflow rates observed for the Type 1 sample in this paper. Figure 11 clearly illustrates that the different assumptions on the outflow geometry (i.e. biconical outflow or thin spherical shell) have by far a smaller impact on the derived value for the mass outflow rate than the assumed value of the gas density. The latter is responsible for more than an order of magnitude difference in the range on n_e explored. Finally, combining the SUPER Type 1 sample presented in this paper with literature data, we confirm a correlation of the mass outflow rate with the AGN bolometric luminosity across five orders of magnitude. This relation could potentially be used to put constraints on the different AGN outflow models, but we caution that the real scatter and the slope of this relation is uncertain as long as we do not have better constraints on the gas density object by objects (see also Davies et al. 2020b).

Fig. 11. Ionised gas [O III] mass outflow rate versus bolometric luminosity of AGN for the SUPER Type 1 sample presented in this paper and literature data of low- and high-redshift AGN. The red shaded area and the black hatched area show the mass outflow rates for the SUPER targets assuming a biconical outflow model and a thin-shell model, as described in Sect. 5.5. The green shaded area shows the outflow rates for ionised gas from literature data compiled in Fiore et al. (2017) (after rescaling the relation with the same assumption as the Type 1 SUPER targets, see Sect. 5.5) and the blue shaded region shows the outflow rates for low redshift X-ray AGN sample from Davies et al. (2020b). The shaded region in all the studies correspond to mass outflow rates assuming an electron density from 500 cm−3–10 000 cm−3.

Finally, as argued in Baron & Netzer (2019), the ionised gas is most likely associated with a significant amount of neutral dusty gas that occurs in the same outflow. Along with the fact that outflows can exist in multiple gas phases (e.g. Perna et al. 2019; Serafinelli et al. 2019), this suggests that the actual mass outflow rates can be substantially larger than the ones deduced here.

6. Scaling relations

We now focus on the relation between the properties of the detected outflows (e.g. velocity and escape fraction) and the properties of the central super massive black hole (e.g. bolometric luminosity and black hole mass). Theoretical models of AGN outflows predict fast winds to originate from the accretion disc that later impact the ISM, resulting in a forward shock that expands within the host galaxy (e.g. Zubovas & King 2012; Zubovas & Nayakshin 2014; Faucher-Giguère & Quataert 2012; Wagner et al. 2013). This would naturally predict positive correlations between outflow properties (e.g. velocity, mass outflow rate) and AGN properties (e.g. King & Pounds 2015). The slope of such relation will in addition depend on the density profile of the medium (e.g. Faucher-Giguère & Quataert 2012). While there are already a few studies exploring such scaling relations for literature samples of AGN across a wide redshift range (e.g. Fiore et al. 2017; Fluetsch et al. 2019), we perform a similar analysis on the SUPER sample for which we have a homogeneous and detailed characterisation of both the outflows properties as well as those of the central black hole.

We first investigated the scaling relations for w₈₀ (Fig. 12) versus the physical properties of the central SMBH. In the top-left panel of Fig. 12, we plot w₈₀ against the AGN X-ray luminosity L_X(2 − 10 keV). A clear positive trend is present, that is, higher X-ray luminosity targets show larger [O III] width. A similar positive correlation is observed if we consider the AGN bolometric luminosity (top right panel in Fig. 12) instead of the X-ray luminosity. To cover a wider range in luminosity, we added the AGN sample observed in the KASHz survey matched in redshift (z ∼ 2–2.5) and the high luminosity targets from the WISSH survey (e.g. Bischetti et al. 2017; Martocchia et al. 2017; Zappacosta et al. 2020). We derived the Pearson coefficient and the p-value for each of these relations and the results are reported in Table 5. For the SUPER data alone, the Pearson coefficient for w₈₀ versus L_X relation is 0.51, while for the redshift matched KASHz sample the coefficient is 0.37 with a null hypothesis probability (for a non-correlation) of ∼2.4% and 2%, respectively. For the combined data sets of the SUPER, KASHz, and WISSH samples, the Pearson coefficient and the null hypothesis probability is 0.59 and < 0.001%, respectively, suggesting a strong correlation between w₈₀ and L_X.

Fig. 12. Top left panel: non-parametric velocity width, w80 derived from integrated spectrum versus X-ray luminosity measured in the 2–10 keV band. The red data points represent the SUPER sample, the grey data points show the luminosity and redshift-matched sample from the KASHz survey, and the green data points show the high-luminosity sample from the WISSH survey. The dashed red line shows the best fit relation for SUPER galaxies, the dashed grey line shows the best fit relation for the matched KASHz sample, and the blue line shows the overall relation for the SUPER, KASHz, and WISSH samples (Bischetti et al. 2017). Top right panel: w80 versus bolometric luminosity. Bottom left panel: w80 versus black hole mass, MBH. Bottom right panel: w80 versus Eddington ratio (Lbol/LEdd). In all the plots, the blue line shows the best fit to all the data points within the plot. The values plotted for the SUPER sample are obtained from Tables 1 and 3.

The statistical tests also confirm strong positive correlation for w₈₀ versus L_bol (Fig. 12, null-hypothesis probability < 0.1%), which may be interpreted as an increase in the turbulence of the gas within the AGN host galaxies with the overall radiation power of the central black hole. We note that there is a positive correlation also for v₁₀ versus L_X and L_bol, although they are weaker, according to the Pearson test, compared to those presented above for w₈₀.

Further exploring the correlation of the gas kinematics with the AGN properties, w₈₀ has a relatively weak correlation with the black hole mass (bottom left panel in Fig. 12) with a null hypothesis of ∼2%. Similar correlation is observed between w₈₀ and the Eddington ratio (λ_Edd = L_bol/L_Edd, where L_Edd = 1.3 × 10³⁸(M_BH/M_⊙ erg s⁻¹)) with a null-hypothesis probability of ∼0.4%. These values suggest that the outflow velocity correlates weakly with the black hole mass when compared with the bolometric luminosity correlation.

Apart from w₈₀, we tested the presence of a correlation in the L_bol versus maximum velocity of the [O III] line (Fig. 13), defined in the literature as v_max = Δλ_broad + 2σ_broad, that is, the sum of the difference between the centroid velocity of the broad and the narrow component and twice the velocity dispersion of the broad component (e.g. Rupke & Veilleux 2013). The Pearson coefficient for the relation between v_max and L_bol is 0.68 with a null hypothesis probability of 0.1%, also confirming in this case a positive correlation between these two quantities. We note that despite the similar proportionality, the normalisation of the curve in the SUPER sample and that in the Fiore et al. (2017) sample are different. We expect that the velocity of the outflowing gas depends on the properties of the medium it is moving through, which would result in some intrinsic scatter in any correlation between such velocity and the AGN bolometric luminosity. The difference in the normalisation can also stem from differences in the gas mass of the host galaxies, as explained below.

Fig. 13. Maximum velocity, vmax, defined in Sect. 6, versus AGN bolometric luminosity. The red circles represent the SUPER targets, the dashed blue line shows the best fit relation for literature data compiled in Fiore et al. (2017). The black curves are the Menci et al. (2019) predictions for the shock velocity (approximated here by the measured vmax) versus AGN luminosity for galaxies. The three curves correspond to molecular gas masses of 2.7 × 109 M⊙, 6 × 109 M⊙ and 2 × 1011 M⊙, respectively. More details are given in Sect. 6.

To test the latter hypothesis, we considered the Menci et al. (2019) model, which provides a 2D description of the expansion of an AGN-driven shock in a disc with an exponential gas density profile. Their predictions are expressed in terms of the total molecular gas mass and circular velocity of the host galaxy. In particular, from Eq. (14) in Menci et al. (2019), we can express the velocity of the shock, V_S, as a function of the AGN bolometric luminosity, L_bol and host galaxy gas mass, M_gas, in the following way:

$$\begin{array}{cc}& log(\frac{V_S}{\mathrm{km}{\mathrm{s}}^{-1}})=2.65+0.35log(\frac{L_{\mathrm{bol}}}{{10}^{45}{\mathrm{erg}\mathrm{s}}^{-1}})\hfill \\ \hfill & -\frac{1}{3}log(\frac{M_{\mathrm{gas}}}{{10}^{10}{M}_{\odot }})+log(\frac{V_c}{200{\mathrm{km}\mathrm{s}}^{-1}})\hfill \\ \hfill & -0.65log(\frac{R_S}{1\mathrm{kpc}}).\hfill \end{array}$$(4)

We used the above formula to predict the scaling of the outflow velocity versus the AGN bolometric luminosity for a range of host galaxies’ properties that encompass, to the best of our knowledge, those of our sample. We assumed that the circular velocity (V_c) of the galaxy is ≈300 km s⁻¹, which is a reasonable assumption from the Tully-Fischer relation for galaxies with M_* ≈ 10⁹ − 10¹¹ M_⊙ (e.g. Cresci et al. 2009). We fixed the radius of the shock (R_S) to ≈1 kpc, which is where the bulk of the outflowing gas is located following our spectroastrometry analysis in Sect. 5.4. We used our CO(3-2) ALMA follow-up observations (Circosta et al., in prep.) to estimate the range of molecular gas masses for our AGN hosts: the mean gas mass for the sample is M_gas ≈ 6 × 10⁹ M_⊙, with 90% of the sample above M_gas ≈ 2.7 × 10⁹ M_⊙ and maximum measured gas mass is M_gas ≈ 2 × 10¹¹ M_⊙. As can be seen in Fig. 13, the range of model tracks (shown by the black relations) is able to reproduce most of our observed values (shown as red circles), where we have assumed that v_max is a good approximation of V_S. Therefore, the Menci et al. (2019) description of an AGN-driven shock moving through a galaxy disc is able to reproduce our observations once the observed range of molecular gas masses for the host galaxies is used. The data points lying outside the model predictions from Menci et al. (2019) can be explained if galaxies with lower or higher molecular gas mass are taken into account.

Active galactic nuclei-driven outflows are suggested as a physical mechanism capable of effectively removing gas from their host galaxies, which is a very important process to reproduce the observed mass build-up of galaxies in the local universe and the chemical composition of the ISM (e.g. Muratov et al. 2015). We can use the [O III] velocity profile to estimate the fraction of the ionised gas that has the ability to escape the gravitational potential of the AGN host galaxy. We first calculated the escape velocity of a cloud, which is at a distance r from the AGN location (see Rupke et al. 2002):

$$\begin{array}{c}\hfill V_{\mathrm{esc}}=\sqrt{2}{V}_c{[1+ln\left(\frac{r_{max}}{r}\right)]}^{1/2},\end{array}$$(5)

where r_max is the extent of the dark matter halo (assumed ∼100 kpc here) and V_c is the circular velocity of the galaxy, assumed to be 300 km s⁻¹ as in Eq. (4) from the Tully-Fischer relation for the range of stellar mass in SUPER galaxies (see Cresci et al. 2009). At 2 kpc, which is roughly the spatial resolution of the H-band SINFONI observations, the escape velocity calculated using Eq. (5) is on average ∼950 km s⁻¹. We then estimated the escape fraction of ionised gas as the ratio between the flux of the [O III] λ5007 channels in the integrated spectra with |v| > V_esc and the total flux of the [O III] λ5007 line. The escape fraction of the ionised gas calculated for the Type 1 SUPER sample presented in this paper is reported in Table 4. For most galaxies, ≲10% of the outflowing ionised gas has the ability to escape the host galaxy, while the rest of the gas is expected to re-accrete back onto the host galaxy. The highest escape fractions of ∼43% and ∼97% are observed in J1333+1649 and J1441+0454, which are also the targets with high bolometric luminosity. We note that in J1441+0454, the exceptionally high escape fraction is a result of the extremely blueshifted [O III] profile, which is also blended with the iron emission as explained in Sect. 5.1. Figure 14 shows the escape fraction as a function of the bolometric luminosity, which suggests a positive correlation between the two quantities, but with a large scatter. The Pearson correlation coefficient between the escape fraction and the bolometric luminosity is 0.64 with a null hypothesis probability of ∼0.2%. Although the correlation is relatively weak compared to the w₈₀ − L_bol relation in Fig. 12, the data might suggest that outflows hosted in high luminosity AGN have a higher fraction of escaping gas.

Fig. 14. Escape fraction as a function of the bolometric luminosity of the AGN for the Type 1 sample presented in this paper. The escape velocity is calculated for a distance of 2 kpc from Eq. (5). The blue curve shows the best-fit relation.

We also note that SINFONI data only trace the ionised phase of the outflow, while a large fraction of the outflow might be present in the molecular gas phase (e.g. Cicone et al. 2018; Fluetsch et al. 2019). Therefore, follow-up studies in other gas phases such as the warm and cold molecular phase and neutral gas phase are required to complete the picture of the effects of radiation pressure on the ISM of the galaxies.

7. Summary and conclusions

We present near-infrared IFU spectroscopy for 21 X-ray selected Type 1 AGN from the SUPER survey (half of the survey size). We traced and characterised the velocity and extension of the ionised gas and its outflow component in the NLR using the [O III] λ5007 transition. The main conclusions from the work presented in this paper are summarised below.

1. Using a cut on the non-parametric velocity width, w₈₀ >  600 km s⁻¹, in the integrated spectra, we find that all the Type 1 AGN in the SUPER survey show the presence of ionised outflows. We consider the selected threshold in w₈₀ a conservative choice, based on the very diverse distributions of low-z and high-z mass matched sample of star forming galaxies compared to the w₈₀ distribution of AGN host galaxies (Sect. 5.1, Fig. 3). We also confirm a strong linear correlation between the [O III] luminosity and the X-ray luminosity, as previously observed in the literature (Sect. 5.1, Fig. 4).

2. Using three different methods (COG analysis, half-light radii, and the PSF-subtraction method), we find evidence of kiloparsec-scale extended ionised gas emission in seven out of the 11 targets for which [O III] is detected at an S/N >  5 in the integrated spectrum (Sect. 5.2, Figs. 5 and 6). Flux and velocity maps of these resolved targets reveal outflows extended to ∼6 kiloparsec (quantified by D₆₀₀ value, Table 4 and Figs. 7 and 8), with indications of redshifted outflows in three objects (cid_346, J1333+1649, and S82X1905).

3. We used the spectroastrometry method to determine the distance of the bulk of the gas moving at a given velocity. We find that the high-velocity outflowing gas (|v| > 600 km s⁻¹) is contained in the central ∼1 kpc for ∼80% of the targets for which the spectroastrometry analysis was performed (Sect. 5.4). For two objects (cid_346 and S82X1905) the bulk motion of high-velocity gas is extended to ∼3 kpc.

4. We explored a range of plausible assumptions on the physical properties of the outflow (its geometry, velocity, and radius) and of the outflowing gas (i.e. its electron density) and report the range of derived mass outflow rates for each target (Sect. 5.5). The mass outflow rates of the Type 1 sample are in the range of ∼0.01–1000 M_⊙ yr⁻¹ (Table 4). After factoring in the systematic uncertainties in the outflow models, these outflow rates seem to correlate with the bolometric luminosity of the AGN.

5. The non-parametric velocity, w₈₀, strongly correlates with the X-ray luminosity, L_X(2 − 10 keV), and the bolometric luminosity of the central super massive black hole (Fig. 12). The correlation is relatively weak for w₈₀ versus M_BH (or L_Edd). The maximum velocity versus bolometric luminosity plot for the SUPER sample agrees with model predictions for an AGN-driven shock driving an outflow through a galaxy disc (Fig. 13).

6. For most galaxies, < 10% of the outflows in the ionised gas have the potenial to escape the gravitational potential of the host galaxy. The escape fraction also increases with the bolometric luminosity of the AGN (Fig. 14).

While this paper focused on presenting the ionised gas kinematics of the Type 1 sample, an upcoming publication will present the SINFONI and ALMA results for the overall (Type 1 and Type 2) AGN SUPER sample. The different types of morphology of the ionised gas and the associated outflows presented in this paper clearly show the advantages of performing high spatial resolution observations. Even better spatial and spectral resolution and higher sensitivity is expected with future facilities such as ELT/HARMONI. Such observations will enable us to resolve the extended redshifted and blueshifted outflowing gas, which might reveal sub-structures within the outflow, and trace the radial evolution of the outflowing gas, which is certainly not possible at z ∼ 2 with today’s instrumentation. Finally, we are currently limited to studying ionised outflows from the ground up to z ≈ 4, observing the [O III] line in the K-band. With the IFU capabilities of NIRSpec on-board the JWST, a study similar to the one presented in this paper could be performed with comparable resolution at z >  4 and will allow us to further constrain the importance of AGN feedback at earlier comic epochs.

Acknowledgments

We thank the referee for the useful and constructive comments. We thank Michele Cirasuolo, Alice Concas and Dominika Wylezalek for providing the comparison data for the w₈₀ distribution. M. P. is supported by the Programa Atracción de Talento de la Comunidad de Madrid via grant 2018-T2/TIC-11715. Based on observations collected at the European organisation for Astronomical Research in the Southern Hemisphere under ESO program 196.A-0377.

References

Appendix A: Integrated H-band spectra and [O III] channel maps of the Type 1 SUPER sample

Fig. A.1. Integrated H-band spectrum of SUPER targets X_N_160_22, X_N_81_44, X_N_66_23, X_N_35_20, X_N_12_26, and X_N_44_64. The grey curve shows the observed spectrum, the red curve shows the reproduced overall emission line model, the blue dashed curve shows the iron emission, and the dashed green curves show the continuum emission and the individual Gaussian components (narrow, broad, and BLR) used to reproduce the profiles of various emission lines. The blue vertical lines indicate the locations of Hβ, [O III] λ4959, and [O III] λ5007. The vertical grey regions mark the channels with strong skylines, which were masked during the fitting procedure. The X-axis shows the rest frame wavelength after correcting for the redshift of the target and the Y-axis shows the observed flux.

Fig. A.2. Same as Fig. A.1 for X_N_4_48, X_N_102_35, X_N_115_23, cid_166, cid_1605, and cid_346.

Fig. A.3. Same as Fig. A.1 for cid_1205, cid_467, J1333+1649, J1441+0454, J1549+1245, and S82X1905.

Fig. A.4. Same as Fig. A.1 for S82X1940 and S82X2058.

Fig. A.5. 1.5″ × 1.5″ [O III] λ5007 channel maps of SUPER target X_N_160_22 at different velocity slices, after subtracting the Hβ, [O III] λ4959, and iron models from the raw cube. Each velocity slice is 200 km s−1 wide, and the displayed value is the centre velocity of the respective channel. The black cross marks the location of the H-band continuum peak, used as a proxy for the AGN position. North is up and east is left.

Fig. A.6. Same as in Fig. A.5 for X_N_81_44, X_N_66_23 and X_N_12_26.

Fig. A.7. Same as in Fig. A.5 for X_N_115_23 and cid_166.

Fig. A.8. Same as in Fig. A.5 for cid_346 and J1333+1649.

Fig. A.9. Same as in Fig. A.5 for J1549+1245, S82X1905, and S82X1940.

All Tables

All Figures

Fig. 1. Past and ongoing IFU surveys targeting the [O III] λ5007 emission line to study the properties of outflows in AGN host galaxies in the Lbol versus redshift plane. The grey data points represent surveys at z <  1 (Husemann et al. 2013, 2014, 2017; Liu et al. 2013; Harrison et al. 2014; Karouzos et al. 2016; Bae et al. 2017; Rupke et al. 2017), while the solid black symbols are collected from z >  1 studies (Nesvadba et al. 2006, 2017; Alexander et al. 2010; Harrison et al. 2012; Cresci et al. 2015; Carniani et al. 2015; Brusa et al. 2016; Kakkad et al. 2016; Vayner et al. 2017; Vietri et al. 2018). The blue shaded area covers the parameter space of the KASHz survey (Harrison et al. 2016) and the red shaded area that of the overall SUPER survey (Circosta et al. 2018). This paper presents an analysis on the Type 1 sub-sample of the SUPER survey, the data points of which are shown as cyan star symbols. Adapted from Circosta et al. (2018), to which we refer the reader for further details.

In the text

Fig. 2. Integrated H-band spectrum of X_N_115_23 shown here as an example. The grey curve shows the observed spectrum, the red curve shows the reproduced overall emission line model, the blue dashed curve shows the iron emission, and the dashed green curves show the continuum emission and the individual Gaussian components (narrow, broad, and BLR) used to reproduce the profiles of various emission lines. The blue vertical lines indicate the location of Hβ, [O III] λ4959 and [O III] λ5007. The vertical grey regions mark the channels with strong sky lines that were masked during the fitting procedure. The X-axis shows the rest frame wavelength after correcting for the redshift of the target, and the Y-axis shows the observed flux density. The integrated spectra of the rest of the targets are shown in Appendix A.

In the text

Fig. 3. Top panel: distribution of the non-parametric velocity, w80 measured for [O III] λ5007 line, for AGN and star forming samples from different surveys. The AGN sample consists of targets from the SUPER survey in red (this paper) and targets from KASHz survey matched in redshift in hatched black (Harrison et al. 2016). The Type 1 sample used in this paper has a median hard X-ray luminosity of 1044.8 erg s−1 and the redshift-matched KASHz targets have a median X-ray luminosity of 1043.9 erg s−1. The blue histogram shows the w80 distribution of mass-matched low redshift star forming sample from Wylezalek et al. (2020), while the ocean-blue vertical bar shows the w80 value from the redshift and mass-matched star forming sample from the KLEVER survey (334 ± 35 km s−1, Curti et al. 2020). The dashed black-line at 600 km s−1 corresponds to the w80 value used in this work to define if a target hosts an AGN-driven outflow. Lower panel: inverse cumulative w80 distribution with the same colour coding as in the top panel. It is clear from both the panels that targets from the SUPER survey occupy the higher end of the w80 distribution. Based on the w80 criteria, all the targets in the SUPER survey show the presence of outflows, and ∼52% of the redshift matched targets in the KASHz survey show outflows.

In the text

Fig. 4. Total [O III] luminosity versus hard X-ray luminosity (LX(2–10 keV) for the SUPER (red circles) and KASHz targets (blue data points), the latter from a similar survey targeting X-ray selected AGN at high redshift. The combined data set of KASHz and SUPER show a positive correlation between the two quantities, shown by the black line. The shaded region represents the 1σ scatter on the relation. KASHz targets occupy the lower X-ray luminosity range compared to SUPER targets, which may explain the difference in the distribution of w80 in the two surveys in Fig. 3.

In the text

Fig. 5. Curves of growth for the Type 1 sample as described in Sect. 5. The solid black line shows the curve from the image obtained by collapsing the PSF-observation cube at the location of [O III] λ5007 channels. The Hβ BLR component is shown as a dashed black curve with the grey area showing the 1σ uncertainly. The BLR component is unresolved and serves as a proxy to determine the PSF during the observations which is also compared with the dedicated PSF star observation. The red curve shows the curve of growth for the total [O III] emission, meaning narrow and broad components. All the uncertainties are at 1σ levels. The BLR component is not plotted for X_N_66_23 as it remains undetected. For the rest of the Type 1 sample in the SUPER survey, the S/N is not enough to perform such a curve-of-growth analysis.

In the text

Fig. 6. Residuals obtained using “PSF-subtraction” method across the SINFONI FoV as described in Sect. 5 for the Type 1 sample. Left panel: residual image obtained by collapsing the channels at the location of [O III] λ5007 after subtracting the emission from the AGN. The colour map is on a linear scale with the yellow regions showing excess emission (i.e. > 3σ). North is up and east is to left. Right panel: spectrum extracted at the location of excess emission, with an aperture of 0.1″ × 0.1″ shown by the red square in the left panel. The vertical grey lines in the residual spectrum shows the location of the [O III] λ4959 and [O III] λ5007, respectively, while the channels used to create the residual map on the left are indicated with a magenta box. Non-zero residuals in the images show that [O III] emission is resolved, while a noisy pattern is indicative of unresolved emission, which are, respectively, marked in the left panels. The AGN X_N_66_23 is not included in this analysis due to non-detection of the BLR component.

In the text

Fig. 7. Flux and velocity maps of targets showing resolved/extended emission inferred from the COG and PSF-subtraction method (from top to bottom: X_N_81_44, X_N_66_23, X_N_115_23, and cid_346). Column 1: background image shows the flux map of total [O III] emission (narrow+broad). Column 2: [O III] line centroid map. Column 3: v10 map. Column 4: w80 map. In all the maps, the black star denotes the centre position of the H-band continuum, as a proxy for AGN location. North is up and east is left.

In the text

	Fig. 8. Same as Fig. 7 for J1333+1649, J1549+1245, and S82X1905.
In the text

Fig. 9. Results of spectroastrometry applied on the Type 1 SUPER targets. The X-axis represents the distance between the continuum centre (vertical dashed line) and the position of the photo-centre of the [O III] line at the corresponding velocity channel given in the Y-axis. The error bars represent 1σ uncertainty on the values. The red arrow marks the cloud with velocity greater than 600 km s−1 (300 km s−1 in the [O III] line profile) whose centroid is at the maximum distance from the AGN.

In the text

Fig. 10. 1.5″  × 1.5″ [O III] λ5007 channel maps of X_N_115_23 at different velocity slices, after subtracting the Hβ, [O III] λ4959 and iron models from the raw cube. Each velocity slice is 200 km s−1 wide and the displayed value is the centre velocity of the respective channel. Black cross marks the location of the H-band continuum peak, used as a proxy for the AGN position, and the cyan star marks the location of the [O III] centroid determined from a 2D Gaussian fit. North is up and east is to the left.

In the text

Fig. 11. Ionised gas [O III] mass outflow rate versus bolometric luminosity of AGN for the SUPER Type 1 sample presented in this paper and literature data of low- and high-redshift AGN. The red shaded area and the black hatched area show the mass outflow rates for the SUPER targets assuming a biconical outflow model and a thin-shell model, as described in Sect. 5.5. The green shaded area shows the outflow rates for ionised gas from literature data compiled in Fiore et al. (2017) (after rescaling the relation with the same assumption as the Type 1 SUPER targets, see Sect. 5.5) and the blue shaded region shows the outflow rates for low redshift X-ray AGN sample from Davies et al. (2020b). The shaded region in all the studies correspond to mass outflow rates assuming an electron density from 500 cm−3–10 000 cm−3.

In the text

Fig. 12. Top left panel: non-parametric velocity width, w80 derived from integrated spectrum versus X-ray luminosity measured in the 2–10 keV band. The red data points represent the SUPER sample, the grey data points show the luminosity and redshift-matched sample from the KASHz survey, and the green data points show the high-luminosity sample from the WISSH survey. The dashed red line shows the best fit relation for SUPER galaxies, the dashed grey line shows the best fit relation for the matched KASHz sample, and the blue line shows the overall relation for the SUPER, KASHz, and WISSH samples (Bischetti et al. 2017). Top right panel: w80 versus bolometric luminosity. Bottom left panel: w80 versus black hole mass, MBH. Bottom right panel: w80 versus Eddington ratio (Lbol/LEdd). In all the plots, the blue line shows the best fit to all the data points within the plot. The values plotted for the SUPER sample are obtained from Tables 1 and 3.

In the text

Fig. 13. Maximum velocity, vmax, defined in Sect. 6, versus AGN bolometric luminosity. The red circles represent the SUPER targets, the dashed blue line shows the best fit relation for literature data compiled in Fiore et al. (2017). The black curves are the Menci et al. (2019) predictions for the shock velocity (approximated here by the measured vmax) versus AGN luminosity for galaxies. The three curves correspond to molecular gas masses of 2.7 × 109 M⊙, 6 × 109 M⊙ and 2 × 1011 M⊙, respectively. More details are given in Sect. 6.

In the text

	Fig. 14. Escape fraction as a function of the bolometric luminosity of the AGN for the Type 1 sample presented in this paper. The escape velocity is calculated for a distance of 2 kpc from Eq. (5). The blue curve shows the best-fit relation.
In the text

Fig. A.1. Integrated H-band spectrum of SUPER targets X_N_160_22, X_N_81_44, X_N_66_23, X_N_35_20, X_N_12_26, and X_N_44_64. The grey curve shows the observed spectrum, the red curve shows the reproduced overall emission line model, the blue dashed curve shows the iron emission, and the dashed green curves show the continuum emission and the individual Gaussian components (narrow, broad, and BLR) used to reproduce the profiles of various emission lines. The blue vertical lines indicate the locations of Hβ, [O III] λ4959, and [O III] λ5007. The vertical grey regions mark the channels with strong skylines, which were masked during the fitting procedure. The X-axis shows the rest frame wavelength after correcting for the redshift of the target and the Y-axis shows the observed flux.

In the text

	Fig. A.2. Same as Fig. A.1 for X_N_4_48, X_N_102_35, X_N_115_23, cid_166, cid_1605, and cid_346.
In the text

	Fig. A.3. Same as Fig. A.1 for cid_1205, cid_467, J1333+1649, J1441+0454, J1549+1245, and S82X1905.
In the text

	Fig. A.4. Same as Fig. A.1 for S82X1940 and S82X2058.
In the text

Fig. A.5. 1.5″ × 1.5″ [O III] λ5007 channel maps of SUPER target X_N_160_22 at different velocity slices, after subtracting the Hβ, [O III] λ4959, and iron models from the raw cube. Each velocity slice is 200 km s−1 wide, and the displayed value is the centre velocity of the respective channel. The black cross marks the location of the H-band continuum peak, used as a proxy for the AGN position. North is up and east is left.

In the text

	Fig. A.6. Same as in Fig. A.5 for X_N_81_44, X_N_66_23 and X_N_12_26.
In the text

	Fig. A.7. Same as in Fig. A.5 for X_N_115_23 and cid_166.
In the text

	Fig. A.8. Same as in Fig. A.5 for cid_346 and J1333+1649.
In the text

	Fig. A.9. Same as in Fig. A.5 for J1549+1245, S82X1905, and S82X1940.
In the text