© ESO 2021

1. Introduction

Feedback mechanisms from active galactic nuclei (AGN) are widely considered to play a key role in galaxy formation and evolution (Silk & Rees 1998; King 2010a; Fabian 2012; King & Pounds 2015). AGN feedback is indeed included in all theoretical, semianalytic, and numerical studies of galaxy formation and evolution (e.g., Granato et al. 2004; Di Matteo et al. 2005; Ciotti et al. 2010; Guo et al. 2016; Nelson et al. 2019) as it allows us to reconcile theoretical predictions with observed galaxy properties. Such AGN activity is considered the main factor responsible for the quenching of star formation in more massive galaxies (so-called −negative feedback). The energy output of a supermassive black hole (SMBH) accreting close to the Eddington limit is large enough to drive massive, wide-angle outflows on large scales (e.g., Zubovas & King 2012, 2014; Costa et al. 2014; King & Pounds 2015; Pontzen et al. 2017) that are capable of either sweeping the gas out of the host galaxy, thus reducing the host galaxy gas reservoir for star formation, or heating the intergalactic medium of the host galaxy through the injection of thermal energy, thus preventing the gas from cooling and collapsing to form stars (‘ejective’ versus ‘preventive’ feedback; e.g., see Woo et al. 2017; Cresci & Maiolino 2018). Both processes are expected to halt the accretion onto the central black hole and, consequently, to give rise to the SMBH mass values observed to correlate with the physical properties of the host galaxy bulge (i.e., its mass, velocity dispersion and luminosity; e.g., Ferrarese & Merritt 2000; Kormendy & Ho 2013). Additional observational evidence supporting the mutual influence of the central SMBH and the host galaxy comes from the observed similarity between the star formation (SF) and BH accretion histories across cosmic time. In fact, both activity histories are seen to peak at z ∼ 2 (Madau et al. 1996; Madau & Dickinson 2014), meaning that the bulk of both SF and BH accretion occurred within z ∼ 1−3 (e.g., Marconi et al. 2004; Aird et al. 2010, 2015). The epoch z ∼ 1−3 (also referred to as ‘cosmic noon’) is hence crucial to the study of such phenomena and their effects in action since this is the time when AGN feedback is expected to be more effective.

Thanks to advanced integral-field spectroscopic (IFS) facilities, AGN-driven outflows have been exhaustively observed from optical to IR and mm bands, in both local (e.g., Feruglio et al. 2010, 2015; Cicone et al. 2012, 2014) and high-redshift galaxies (e.g., Cano-Díaz et al. 2012; Maiolino et al. 2012; Carniani et al. 2015; Cresci et al. 2015). It is worth mentioning that while theoretical predictions usually refer to the whole outflowing gas, observations usually trace the emission produced by a single gas phase of the outflow. Therefore, in order to properly compare model predictions with observational results, it is fundamental to obtain a complete, multi-phase description of the outflow (e.g., Cicone et al. 2018; Harrison et al. 2018).

Even though the existence of AGN-driven outflows has been widely confirmed thanks to observations, there are a number of relevant open questions that remain unanswered, considering, for instance, how the energy released by the accreting BH is coupled with the interstellar medium (ISM), thus driving large-scale outflows, and how efficient the coupling is between the nuclear and galaxy-scale outflows.

Theoretical models (e.g., King 2003, 2005; King & Pounds 2015) predict fast (v ∼ 0.1c), highly ionised winds, accelerated on sub-pc scales by the AGN radiative force as the origin of the strong galaxy-scale feedback. As the nuclear wind impacts the ISM of the host galaxy, it produces an inner reverse shock that slows down the wind, along with an outer forward shock accelerating the galactic ISM. Depending on the efficiency of cooling processes (typically radiative) in removing energy from the hot shocked inner gas, there are two main wind driving-modes (e.g., King 2010b; Costa et al. 2014; King & Pounds 2015). If the cooling occurs on a timescale that is shorter than the wind flow time, most of the inner wind kinetic energy is lost (usually via inverse Compton scattering) and, therefore, the wind momentum is the only conserved physical quantity (‘momentum-driven’ regime). Vice versa, if the cooling is negligible, the postshock gas retains all the mechanical energy and expands adiabatically (‘energy-driven’ regime), sweeping up a significant amount of the host galaxy gas. According to a widely accepted picture, the observed scaling relations are the result of the effect of AGN feedback acting on the host galaxy through two distinct, subsequent phases (e.g., Zubovas & King 2012; King & Pounds 2015): an initial momentum-driven regime lasting so long as the BH mass has not yet reached the M_BH − σ relation and the outflow is confined within ∼1 kpc from the central BH, followed by a later energy-driven phase after the BH mass has settled on the relation, during which the outflow can propagate beyond 1 kpc scales.

From an observational point of view, the most promising candidates for acting as the ‘engine’ of the large-scale feedback are the so-called ultra-fast outflows (UFOs), which are highly-ionised, accretion disk winds with mildly relativistic velocities that originate at sub-pc scales. They are usually detected in AGN X-ray spectra (Chartas et al. 2002, 2016; Gofford et al. 2013, 2015; Nardini et al. 2015) via the presence of strongly blueshifted absorption lines of highly-ionised metals (typically iron, e.g., Fe XXV and Fe XXVI). UFOs are found in at least 40% of the local sources (Tombesi et al. 2010, 2011, 2012, 2013), with typical mass outflow rates of ∼0.01−1 M_⊙ yr⁻¹ and kinetic powers of log Ė_K ∼ 42 − 45 erg s⁻¹ (Tombesi et al. 2012).

Moving to high redshifts (z > 1), the UFO detection is hampered by the resolution or sensitivity limits of the current observational facilities. In fact, the number of AGN hosting UFOs at z > 1 that have been discovered thus far drastically decreases to 14 objects (see Dadina et al. 2018 on the list of published objects, and Chartas et al. 2021 for the latest updates); of these, twelve are gravitationally lensed systems (Hasinger et al. 2002; Chartas et al. 2003, 2007, 2009, 2016, 2021; Dadina et al. 2018). Strong gravitational lensing is indeed a well-known, powerful tool to investigate the physical properties of distant quasars (QSOs). The magnified view delivered by gravitational lenses allows us to separate the active nuclei from their hosts, enabling new measurements and spatially resolved studies, which otherwise would not have been possible beyond the local Universe (e.g., Peng et al. 2006; Ross et al. 2009; Bayliss et al. 2017; Spingola et al. 2020; Stacey et al. 2021).

To test theoretical predictions and to shed light on the acceleration and propagation mechanisms of large-scale outflows, we need to compare the energetics of the X-ray nuclear UFO with that of the large-scale outflow. In this work, we focus on the ionised phase of large-scale outflows traced by the optical emission of the [O III]λλ4959,5007 line doublet. Given that it is a forbidden transition, it preferentially traces the emission originating from the kpc-scale, typical of the AGN Narrow Line Region (NLR), since it cannot be produced on the high-density (n ∼ 10¹⁰ cm⁻³), sub-pc scale of the AGN broad line region (BLR). In the presence of outflows, the [O III] line profile is highly asymmetric with a broad, blueshifted wing corresponding to high speeds along the line of sight (v ≳ 1000 km s⁻¹ see e.g., Carniani et al. 2015; Cresci et al. 2015; Perna et al. 2015; Brusa et al. 2016; Marasco et al. 2020).

This work is aimed at studying the connection between nuclear, X-ray UFOs, and the ionised phase of large-scale outflows, for the first time, in two QSOs close to the peak of AGN activity (z ∼ 2). We use the [O III]λ5007 emission line to trace the ionised outflow, along with results on X-ray UFOs from the literature to properly compare the wind energetics on different scales. In doing so, we aim to highlight (and take advantage of) the crucial role of gravitational lensing as powerful tool for overcoming the current observational limits and going on to investigate the physical properties of distant QSOs.

This paper is organised as follows. In Sect. 2, we present the selected targets and describe the observation and data reduction procedure. In Sect. 3, we show our data analysis and spectral fitting. The inferred results are then presented in Sect. 4. In Sect. 5, we discuss the wind acceleration mechanism in our two QSOs and compare our results with findings from the literature. Finally, in Sect. 6, we outline our conclusions. We adopt a ΛCDM flat cosmology with Ω_m, 0 = 0.27, Ω_Λ, 0 = 0.73 and H₀ = 70 km s⁻¹ Mpc⁻¹ throughout this work.

2. Description of the observed QSOs

2.1. Selection of targets

Our sample consists of two z ∼ 1.5 multiple lensed QSOs, HS 0810+2554 and SDSS J1353+1138, observed with the Spectrograph for INtegral Field Observations in the Near Infrared (SINFONI, Eisenhauer et al. 2003) at the ESO Very Large Telescope (VLT) Unit Telescope 3 (UT3) within the framework of the program 0102.B-0377(A) (PI: G. Cresci). These objects were specifically selected as they are known to host UFOs (Chartas et al. 2016, 2021) and to be at redshifts (z ∼ 1.5 and z ∼ 1.6 for HS 0810+2554 and SDSS J1353+1138, respectively), such that the optical rest-frame emission (tracing ionised outflows) falls in the range of wavelengths observed by SINFONI, namely, in the near-IR J-band (λ ∼ 1.1−1.4 μm). Consequently, this selection in redshifts corresponds to study objects at epochs close to the peak of AGN activity (z ∼ 2). In total, there are fourteen high redshift (z > 1) QSOs with UFO detection, of which seven are found in the literature (among them, HS 0810+2554; see Dadina et al. 2018 for an updated list), while the rest have not yet been published (including SDSS J1353+1138; Chartas et al. 2021). These are among the brightest – in terms of 2–10 keV luminosity (L_2−10 keV > 10⁴⁵ erg s⁻¹, except for PID352 L_2−10 keV ∼ 10⁴⁴ erg s⁻¹) – QSOs at high redshift; this is either because they are intrinsically luminous (of the total 14-QSO sample, only HS 1700+6416 and PID352; Lanzuisi et al. 2012; Vignali et al. 2015) or because they are subject to gravitational lens magnification (including APM 08279+5255, PG1115+080, H1413+117, HS 0810+2554 and MG J0414+0534¹; Hasinger et al. 2002; Chartas et al. 2003, 2007, 2009, 2016; Dadina et al. 2018). Thus, these objects deliver high quality X-ray spectra which clearly exhibit UFO absorption features, in spite of their high redshift (z > 1). Amongst this z > 1 sample, APM 08279+5255 (z ∼ 3.9; Hasinger et al. 2002) has been the only z > 1 QSO known to host a large-scale (molecular) outflow (Feruglio et al. 2017), hence, this is where a connection between UFO and galaxy-scale outflow has been put forward thus far. Therefore, within the field of large-scale ionised outflows in HS 0810+2554 and in SDSS J1353+1138, this work will extend the number of z > 1 QSOs for which the connection between nuclear and large scales has been accessed. In the following, we provide a short description of HS 0810+2554 and SDSS J1353+1138, with their main properties listed in Table 1.

2.1.1. HS 0810+2554

HS 0810+2554 is a radio-quiet, narrow absorption line (NAL; FWHM ≲ 500 km s⁻¹) QSO at z ∼ 1.5, which was discovered by Reimers et al. (2002). It consists of four lensed images in a typical fold lens configuration with the two southern, brightest images in a merging pair configuration (A+B), as shown in the HST image in Fig. 1 (left panel). The lens galaxy (labelled with G) is detected in the HST image, and its redshift is estimated to be z_l ∼ 0.89 from the separation and the redshift distribution of existing lenses (Mosquera & Kochanek 2011). Quadruply lensed QSOs occur in strong gravitational lensing regimes (e.g., Narayan & Bartelmann 1996), when the compact and bright UV accretion disk and X-ray corona emission regions overlap the lens caustics. This leads to high magnification factors, whose values strongly depend on the image and lens positions. As a consequence, a small change in the input parameters to the lens models (image and lens positions) can lead to a significant change in the image magnifications. For HS 0810+2554, estimates of the magnification factor μ in different spectral bands are found in the literature, in particular for the X-ray (μ ∼ 103; Chartas et al. 2016), optical (μ ∼ 120; Nierenberg et al. 2020), and radio emission (μ ∼ 25; Jackson et al. 2015).

Fig. 1. Lensed images of HS 0810+1154 (left) and SDSS J1353+1138 (right). Left: HST ACS F555W image of HS 0810+2554 showing the four magnified images of the background quasar in fold lens configuration: the C and D images are spatially resolved, while the pair A+B is blended together. At the centre, we can see the emission from the foreground lens galaxy. Right: V and H-band images of SDSS J1353+1138 taken at the UH88 telescope (upper panels) and corresponding images after the subtraction of A and B components (lower panels), clearly showing the lens galaxy (component G). Image from Inada et al. (2006).

HS 0810+2554 was singled out as an exceptionally X-ray bright lensed object during an X-ray survey of NAL-AGN with outflows of UV absorbing material (Chartas et al. 2009). More recent Chandra and XMM-Newton observations (Chartas et al. 2016) provided definitive proofs for the presence of a highly ionised, relativistic wind in the source nuclear region. The strongly blueshifted absorptions of highly-ionised metals (i.e., Fe XXV, Si XIV) indicate that the outflow velocity components are within the range of 0.1−0.4c. The VLT/UVES spectrum of HS 0810+2554 also shows blueshifted absorptions of the C IV and N V doublets, indicating the existence of UV absorbing material moving with an outflowing speed of v_C IV = 19 400 km s⁻¹ (Chartas et al. 2014, 2016). Even though is classified as radio-quiet object, VLA observations at 8.4 GHz (Jackson et al. 2015) indicate that HS 0810+2554 hosts a radio core, thus demonstrating that it is not radio-silent.

HS 0810+2554 was also recently observed with ALMA in the mm-band (Chartas et al. 2020). The analysis of ALMA data has shown the tentative detection of high-velocity clumps of CO(J = 3 → 2) emission, suggesting the presence of a massive molecular outflow on kpc-scales. With our characterisation of the ionised outflow in HS 0810+2554, we now have, for the first time ever, a three-phase description of an AGN-driven wind at high redshift, from the nuclear to the galaxy scale: the highly-ionised (on nuclear scales), ionised, and neutral molecular phases (on galaxy scales) thanks to a broadband spectral coverage ranging from the X-ray to the optical and mm bands.

2.1.2. SDSS J1353+1138

Unlike HS 0810+2554, which has been widely observed in several spectral bands, SDSS J1353+1138 has been less intensively studied, as its discovery is more recent (Inada et al. 2006). This object was selected from the Sloan Digital Sky Survey (SDSS) as candidate double lensed QSO at z ∼ 1.6. Inada et al. (2006) during the University of Hawai’i 88-inch Telescope (UH88) follow-up observations of SDSS J1353+1138, obtaining V, R, I, and H-band images of the source. The two lensed images are well-distinguishable (see the right panel in Fig. 1), with an angular separation of $\mathrm{\Delta }\sim 1\stackrel{″}{.}40$ (Inada et al. 2006).

More recently, on 2016 January 13, SDSS J1353+1138 was observed with XMM-Newton. The analysis of the X-ray spectrum (Chartas et al. 2021) revealed a significant absorption at ∼6.8 keV (consistent with Fe XXV), indicating the presence of a ∼0.31c UFO.

2.2. SINFONI observations and data reduction

SINFONI observations of HS 0810+2554 and SDSS J1353+1138 were carried out on two different nights in February and March 2019, respectively, in the near-IR J-band (λ ∼ 1.1−1.4 μm) and with a spectral resolution R = 2000. The observations were performed in seeing-limited mode², using the $0\stackrel{″}{.}250\times 0\stackrel{″}{.}125$ pixel scale which provides a total field of view (FOV) of 8″ × 8″, which is essential for mapping the gas dynamics on galaxy scales. The airmasses are different for each target, spanning the ranges of ∼1.7−1.9 and ∼1.2−1.3 during the observations of HS 0810+2554 and SDSS J1353+1138, respectively. The data were obtained in eight and sixteen integrations of 300 s each, for a total of 40 min for HS 0810+2554, and 80 min for SDSS J1353+1138. During each observing block, an ABBA pattern was followed: the target was put alternatively in two different positions of the FOV about $4\stackrel{″}{.}3$ apart, to perform the sky subtraction through a nodding technique. A dedicated star observation to measure the point-spread-function (PSF) was not available in either case but the estimated angular resolution is $\sim 0\stackrel{″}{.}7$ ($\sim 0\stackrel{″}{.}8$) for HS 0810+2554 (SDSS J1353+1138), based on the measured extent of the spatially unresolved BLR emission (see Sect. 3.3). Finally, a standard B-type star for telluric correction and flux calibration was observed shortly before or after the on-source exposures.

We used the ESO-SINFONI pipeline (v. 3.2.3) to reduce the SINFONI data. Before flux calibration and co-addition of single exposure frames, we corrected for atmospheric dispersion effects consisting in a significant change of the AGN continuum emission across the FOV of both sources. This is a consequence of the differential atmospheric dispersion at different wavelengths, which makes it so the measured spectra not ‘straight’ as expected. In practical terms, as the wavelength increases, an increasingly larger fraction of the emission gets deposited into adjacent pseudo-slits, producing a coherent spatial shift of the target as a function of λ. Additionally, possible flexures in the instrument may contribute to producing the observed shift. In order to limit the impact of these optical distortions, we spatially aligned the emission centroid channel-by-channel in each single-exposure, sky-subtracted cube by adopting the following procedure.

For each cube, we first determined, in every spectral channel, the average position of the emission centroid on the FOV through a 2D-Gaussian fit. Then we calculated the shift of the centroid mean position with respect to the centroid position in the first spectral channel, assumed as the reference channel. In both spatial directions on the FOV we found an increasing trend of the shift at increasing wavelengths, which we modelled with a two-degree polynomial to neglect the presence of some spikes that are due to noisier channels. The spatial shift, totally observed from the bluest to the reddest spectral channel, spanned the range of ∼0.5−1 pixel among the various single-exposure cubes of both targets. As the spatial shifts were fractional in units of pixels, we adopted the Drizzle algorithm (Gonzaga et al. 2012; Fruchter & Hook 2019) to perform the optimised alignment of every spectral channel in each raw single-exposure data cube.

Finally, we performed the flux calibration and the co-addition of the single-exposure cubes. The final sky-subtracted, flux-calibrated data consist of 100 × 72 × 2234 data cubes, hence, each one including more than 7000 spectra. Each spectrum is sampled by 2234 channels with a 1.25 Å channel width and covers the spectral range ∼1.1−1.4 μm, corresponding to about 4400−5600 Å rest-frame wavelengths.

2.3. Lens models for the two quasars

As both objects are gravitationally lensed QSOs, lens models are required in order to infer the intrinsic (i.e., unlensed) physical properties of the outflow, such as the intrinsic radius and unlensed flux, which are key ingredients for calculating the outflow energetics.

Both QSOs are lensed by a foreground elliptical galaxy and detailed lens models for these two objects can be found in the literature. In particular, for HS 0810+2554, there are several lens models reported in literature obtained from observations performed in different spectral bands (e.g., from VLA-radio data in Jackson et al. 2015, from ALMA-mm data in Chartas et al. 2020). In this work, we adopted, for HS 0810+2554, the most recent model by Nierenberg et al. (2020), inferred from HST-WFC3 IR observations: images and lens galaxy positions have been measured from direct F140W wide imaging (central wavelength ∼1392 nm), while slitless dispersed spectra have been provided by the grism G141 (useful range: 1075−1700 nm). Assuming z_l ∼ 0.89 for the lens galaxy (Mosquera & Kochanek 2011), Nierenberg et al. (2020) modelled the deflector mass distributions with a singular isothermal ellipsoid (SIE), plus an external shear to account for tidal perturbations from nearby objects.

Detailed lens models for SDSS J1353+1138 are presented in Inada et al. (2006) and Rusu et al. (2016), based on imaging observations in the i-band with the Magellan Instant Camera (MagIC) on the Clay 6.5 m Telescope and in the K-band with the Subaru Telescope adaptive optics system, respectively. Inada et al. (2006) modelled the lens mass distribution using either a SIE model, or a singular isothermal sphere (SIS) model plus a shear component (γ), and estimated the lens redshift to be z_l ∼ 0.3 based on the Faber–Jackson relation (Faber & Jackson 1976). The resulting total magnification factors μ are 3.81 and 3.75 from the SIS+γ and SIE model, respectively. Also assuming z_l ∼ 0.3, Rusu et al. (2016) found slightly lower values for total magnification: μ ∼ 3.47 (SIS+γ), μ ∼ 3.42 (SIE) and μ ∼ 3.53 (SIE+γ). We used all these magnification values from the literature to determine the unlensed flux carried by the ionised outflow in SDSS J1353+1138. Instead, for HS 0810+2554, we were able to provide an estimate of the ionised outflow magnification (μ_out ∼ 2) starting from our data. Such values will be discussed further in Sect. 4.2 and Appendix A.

3. Data analysis

3.1. Fitting procedure

For the spectral analysis of SINFONI data, we adopted the fitting code used in Marasco et al. (2020) to analyse MUSE data of two local QSOs. Here, we implemented the code to also handle the SINFONI data, introducing adjustments and new functionalities depending on the specific necessities of our data. In the following, we illustrate the basics of our spectral-fitting method, highlighting the required changes for the analysis of our SINFONI data (see Marasco et al. 2020 for the detailed description of the fitting code).

The entire fitting procedure aims at performing the kinematical analysis of the diffuse ionised gas, with a primary focus on the [O III]λ5007 emission line (hereafter [O III]) that is the optimal tracer of ionised outflows, as previously mentioned. Our strategy consists of the following three key-steps.

Phase I. We built a template model for the bright BLR emission using an integrated, high signal-to-noise ratio (S/N) spectrum.

Phase II. The BLR template built in phase I was used to map spaxel-by-spaxel the contribution of the BLR to the emission across the entire FOV. The resulting BLR model cube was then subtracted from the data cube.

Phase III. Finally, we performed, spaxel-by-spaxel, a finer modelling of the faint emission lines produced by the diffuse ionised gas, originating on galactic scales.

Hereafter, we refer to the diffuse gas emission lines as ‘narrow’ in order to distinguish them from the typical ‘broad’ emission lines (FWHM > 1000 km s⁻¹; e.g., Osterbrock 1981) originating from within the dense and highly turbulent BLRs.

A single noise value has been associated to each channel in our SINFONI data cubes, computed as the root mean square (rms) of the fluxes extracted spaxel-by-spaxel in a region with no significant emission from the target. The details of the spectral analysis of SINFONI data of both QSOs are provided below.

3.1.1. Modelling the BLR emission

The fitting code starts with modelling the bright BLR emission in a spectrum extracted from the nuclear region, while also fitting the other spectral components. The spectral components to be fitted are: the AGN continuum, BLR emission lines, and narrow emission lines from the diffuse gas. In principle, we should also have the stellar continuum emission, but in our data the AGN continuum is entirely dominant. The BLR model is built by the code as the sum of two independent components: the broad Balmer hydrogen emission lines (Hβ in HS 08010+2554, Hβ and Hγ in SDSS J1353+1138) and several Fe II broad emission lines, which are the two main BLR components in the rest-frame range observed by SINFONI (λ ∼ 4200−5600 Å). The diffuse emission instead consists of the [O III] emission doublet and the narrow components of the Balmer hydrogen lines.

For HS 0810+2554, we extracted a high-S/N spectrum from an aperture of $0\stackrel{″}{.}3$ radius, centred on the observed blended emission of A+B images (see Fig. 1), and fitted all the previously mentioned components simultaneously. We modelled the AGN continuum through a 1st-degree polynomial. The Fe II emission lines were modelled using the semi-analytic templates of Kovačević et al. (2010), while the BLR component of Hβ was fitted by two broad Gaussian components. The narrow emission lines were fitted through two Gaussian components. Given the complexity of the BLR–Hβ line profile in HS 0810+2554, we additionally associated spatially unresolved residuals from the fit to this component, following the approach detailed in Marasco et al. (2020).

In case of SDSS J1353+1138, where the two lensed images (A and B in Fig. 1) are well-distinguishable and spatially resolved, matters were complicated by the fact that we observed a significant change between nuclear spectra extracted from the two distinct images. As a consequence, this prevented us from considering a single BLR template. The procedure followed for SDSS J1353+1138 is described separately in Sect. 3.2.

3.1.2. Mapping the unresolved BLR emission across the FOV

As the BLR emission goes unresolved in our data, we expect that its observed spatial variation follows the PSF of our observations. Therefore, we allowed the BLR template obtained in phase I to change only in amplitude across the FOV and we proceeded to fit the whole data cubes with the software pPXF (Cappellari 2017). For the modelling of the narrow emission lines, we used multiple Gaussian components and adopted a statistical approach (a Kolgomorov–Smirnov test) to select, spaxel-by-spaxel, the minimal, optimal number of Gaussian components to aptly reproduce the emission line profiles (see Marasco et al. 2020 for details). For both HS 0810+2554 and SDSS J1353+1138, we considered a number of Gaussian components ranging from 1 to 3, finding the latter optimal for reproducing the most complex line profiles. Then we subtracted spaxel-by-spaxel the BLR and AGN continuum emissions and, for each galaxy, we created a cube containing only the residual emission lines due to the diffuse gas. In the following, we refer to this cube as the ‘subtracted cube’.

3.1.3. Modelling the narrow emission lines

In phase III, we focused on the finer modelling of the narrow emission lines that remained after the subtraction of the BLR and AGN continuum emission components. The only significant residual emission in our data is the [O III] emission doublet, while the narrow components of Balmer hydrogen lines are very weak and marginally resolved. Therefore, we modelled the residual narrow emission lines through multiple Gaussian components, adopting some reasonable constraints: the two emission lines of the [O III] doublet were fitted by imposing the same central velocity and velocity dispersion, with the intensity ratio I(5007)/I(4959) fixed at 3 according to the theoretical expectations of the atomic theory; whereas, for the weak narrow components of Balmer hydrogen lines, we assumed the same [O III]λ5007 line profile shape and left the flux as a free parameter. This is a reasonable assumption as we expect that the narrow Balmer hydrogen and [O III] emission lines come from regions with the same gas kinematics.

Similarly to the fit in phase II, we considered a number of Gaussian components ranging from 1 to 3, using, as before, a statistical approach to determine spaxel-by-spaxel the optimal number of components required. In both QSOs, most of the line profiles are well reproduced by two Gaussian components: a narrow, bright component close to the systemic velocity, plus a broad blueshifted component to reproduce the [O III] blue wing observed in most of the FOV, which we identify with approaching outflow emission. In HS 0810+2554, the 3-Gaussian fit was selected in some spaxels to reproduce properly the faint but still visible red wing in the [O III] line profile, tracing the fainter outflow component receding from us. On the contrary, in SDSS J1353+1138, two Gaussian components are sufficient to reproduce the most complex line profiles, as we do not detect the [O III] red wing anywhere.

In order to study the physical and dynamical properties of the outflow emission, which is the main focus of this work, we had to properly identify the [O III] emission due to the high-speed outflowing gas by disentangling its contribution to the [O III] emission line due to the gas bulk motion within the host galaxy. To do so, we adopted the same selection criterium used by Marasco et al. (2020). For each Gaussian component used to reproduce the [O III] line profile in a given spaxel, we focused on the fraction of the total line flux contained in the line wings with a velocity shift |v − v_peak| larger than a certain threshold width w_th, where v_peak is the peak velocity of the line in each spaxel. If the fraction of total flux in the line wings was higher than a given threshold τ, the Gaussian component has been classified as a possible ‘outflow’ component, to be confirmed by the following kinematical analysis (Sect. 4.1); otherwise, it has been classified as a ‘narrow’ component, due to systemic gas motions in the host galaxy. We verified the decomposition in several representative spaxels to select the optimal threshold values. In HS 0810+2554, we used τ = 0.5 and w_th = 300 km s⁻¹: the Gaussian component reproducing the narrowest, brightest emission near the systemic velocity has been typically classified as narrow; while any additional Gaussian component used to model either the blue or the red wing in [O III] profile has been identified as outflow component. In SDSS J1353+1138, in order to get the expected classification we slightly relaxed the width threshold (w_th = 250 km s⁻¹).

Figure 2 summarises our strategy. The top panels show J-band spectra of HS 0810+2554 and SDSS J1353+1138, extracted from SINFONI data cubes with an aperture of $\sim 0\stackrel{″}{.}44$-radius, centred on the peak of the observed emission (located on image A in SDSS J1353+1138). The best-fit models shown were obtained in phase II. Since no distinctions between a possibly broad outflow and narrow systemic components were made at this stage of the procedure, the diffuse gas model is plotted as sum of multiple components. We notice that in HS 0810+2554, there is a faint, broad emission line at ∼4700−4750 Å (rest-frame) present in the Fe II templates, which is not present in our data. However, this does not affect the overall fitting procedure, as the observed Fe II emission is reproduced well by the templates at all the other wavelengths. The middle panels show the spectra extracted from the subtracted cubes using the same aperture as above, along with the results from the finer multi-Gaussian fit of the diffuse gas emission lines (phase III). In both QSO-subtracted spectra, the [O III] line profile exhibits a prominent, asymmetric blue wing that is already visible in the full spectra shown above. This strongly suggests the presence of high-speed outflowing material moving towards the observer, which we discuss in greater detail in Sect. 4.1.

Fig. 2. Representative scheme of our fitting procedure. Top panels: SINFONI J-band spectra of HS 0810+2554 (left) and SDSS J1353+1138 (right), zoomed in the spectral region of [O III] and Balmer hydrogen emission lines. Both spectra were extracted using an aperture of $\sim 0\stackrel{″}{.}44$-radius, centred on the peak of the AGN continuum emission (located on image A in SDSS J1353+1138). Black dashed lines show the data, while red solid lines show the best-fit models obtained in phase II of our analysis. The latter is partitioned into the contributions of various components: AGN continuum (yellow), BLR emission in Balmer hydrogen emission lines (green) and Fe II (purple), and narrow line emission from the diffuse gas (light blue) represented as sum of multiple Gaussian components. Middle panels: J-band spectra extracted from the subtracted cubes, with the same aperture used in the top panels. Subtracted data (black lines) are compared to the best-fit models (red lines) resulting from the finer emission lines (EL) modelling implemented in phase III. The green lines highlight the outflow component alone. Bottom panels: residuals obtained by subtracting the full EL model from the subtracted spectra.

3.2. Modelling a ‘double’ BLR in SDSS J1353+1138

As noted in Sect. 3.1.1, for SDSS J1353+1138, we found a significant change in the spectral shape within the wavelength range including the Hβ and [O III] lines, while comparing the nuclear spectrum of the brighter image A (spectrum A) with that of the fainter image B (spectrum B), shown in the upper and lower panels of Fig. 3, respectively. In particular, while in the former it is possible to easily identify the [O III] emission lines, we did not detect any counterpart in the latter. Moreover, the Hβ line profile in spectrum B is broader, with an evident brighter blue wing. Both effects are likely due to an overall increase in the Fe II emission in image B, as the Hβ line width is not expected to intrinsically vary between different lensed images. The anti-correlation between Fe II and [O III] emissions in AGN spectra reflects a well-known effect known as Eigenvector-1 (Boroson & Green 1992) and represents one of the most frequent differences among AGN properties. Even though it has been the subject of many studies, a clear physical understanding of its origin is still lacking. Boroson & Green (1992) suggest that high column densities in the BLR enhance Fe II, while reducing the ionising radiation able to reach the NLR. In a spectral analysis of AGN principal components in SDSS, Ludwig et al. (2009) argue instead that the covering factor of the NLR is the likely cause of the range in [O III] strength, while Ferland et al. (2009) suggest that the higher column densities required for the infall in more luminous AGNs would additionally account for the observed correlation of Fe II strength with L/L_Edd.

Fig. 3. Best-fit models of nuclear spectra of SDSS J1353+1138 extracted from an aperture of $0\stackrel{″}{.}3$ radius, centred on image A (upper panel) and on image B (lower panel). The various spectral components, the total model and data are represented with different colours (see the plot legend). In spectrum A, a broken power law distribution is perfectly suited to reproduce the BLR–Hβ profile, while in spectrum B, an additional broad Gaussian component was required to adequately reproduce the broad peak in the BLR–Hβ line profile, which is entirely dominant over the barely detected Hβ narrow component (solid light-blue line). The two spectra clearly differ from each other mostly for the lack of [O III] detection and the presence of a prominent blue wing in the Hβ line profile in spectrum B.

In spite of its still unclear origin, there are two main possible explanations for the observed significant variation in the Fe II emission between the two lensed images. The first one is based on the typical short time scales (i.e., days, weeks; e.g., Kaspi et al. 2000) on which the BLR is seen to vary. Because of the different path followed by the light from the background QSO, the two lensed images are produced with a time delay of about 16 days (Inada et al. 2006). This temporal shift is comparable with the typical BLR variation timescale, therefore, it could be sufficient to have a significant change in the BLR emission explaining the effect we observed. Given the short time scale probed here, this could carry interesting implications for the accretion variations in the AGN and in the consequent response of the BLR gas. Alternatively, the observed variation could be the consequence of microlensing effects (e.g., see Nierenberg et al. 2020) produced by either single stars or low-mass dark matter halos intervening along the line of sight. Microlensing effects typically affect only the emission originated on small scales, while the emission from the NLR is insensitive. Of the two possibilities, the latter seems to be less likely, as we do not observe any significant counterpart variation in the strength of the BLR–Hβ component, in addition to that observed in the Fe II strength. However, a remarkable simultaneous variation in both Hβ and Fe II strength is what we would expect in the case when the two emissions are strictly co-spatial, whereas we know that the BLR is stratified and that microlensing magnifies the emission from the most compact regions more strongly. Therefore, we do not exclude the microlensing hypothesis. A detailed analysis of the different BLR spectra from the two images is beyond the scope of this work and will be presented in a forthcoming paper. Therefore, in this work, we focus on how we accounted for this effect during the spectral analysis.

Consequently, in SDSS J1353+1138, we extracted two distinct nuclear spectra, namely, spectrum A and B that are shown in Fig. 3, using a 0$\stackrel{″}{.}$3 radius aperture centred on the emission peak of each lensed image. We proceeded to fit them separately. In both spectra we used a first-degree polynomial to model the AGN continuum that is still dominant over the stellar continuum. Unlike the BLR modelling of HS 0810+2554, the multiple Gaussian fit was not sufficient to reproduce the broader and complex profile of the broad Balmer emission lines, especially the Hβ line in spectrum B. In fact, even though both Hβ prominent wings are likely due to the Fe II emission, as discussed above, the Fe II templates employed by the code were not able to reproduce such observed emission. Therefore, we modelled both wings as part of the broad Hβ line profile without any focus on their physical interpretation, as we were simply interested in identifying the overall BLR spectrum in order to remove it. For the modelling of the broad Balmer emission lines observed in SDSS J1353+1138 (i.e., Hβ and Hγ), we used a broken power law distribution convolved with a Gaussian profile (Nagao et al. 2006):

$$\begin{array}{c}\hfill F(\lambda )=\{\begin{array}{cc}{F}_0\times {\left(\frac{\lambda }{{\lambda }_0}\right)}^{+\alpha },\hfill & \text{for}\lambda <{\lambda }_0\hfill \\ F_0\times {\left(\frac{\lambda }{{\lambda }_0}\right)}^{-\beta },\hfill & \text{for}\lambda >{\lambda }_0\hfill \end{array},\end{array}$$(1)

where the free parameters of the fit, for each line, are the central wavelength λ₀, the two power-law indices α and β, the normalisation F₀, and the width σ of the Gaussian kernel. In spectrum A (upper panel), the Hβ and Hγ were modelled separately through the line profile described in Eq. (1). Spectrum B (lower panel) required instead an additional broad Gaussian component to suitably reproduce both the more extended red wing and the broad peak (σ ∼ 800 km s⁻¹) in the Hβ profile. Moreover, we had to constrain the Hγ line profile through that of Hβ. Different Fe II templates have been selected in the BLR best-models for the two nuclear spectra. To model the narrow emission line profiles, we used three Gaussian components in spectrum A, while a single Gaussian component in spectrum B, since we did not actually observe any narrow component.

Then we proceeded to fit the whole data cube following the procedure described in Sect. 3.1.2, with the main difference that in each spatial pixel, pPXF considered a linear combination of the two BLR models, weighing their relative contribution and providing their most suitable combination as the BLR model in that specific spaxel. In general, in those spaxels close to one of the lensed images, we basically got the BLR model directly obtained in the modelling of the respective nuclear spectrum; while a combination of the two BLR-models in those spaxels ended up roughly between the two images, as we indeed expected.

3.3. Testing the spatially resolved emission of the ionised outflows

Before analysing the outflow kinematics across the FOV, we tested whether the emission of the detected ionised outflows was spatially resolved. This is crucial for the calculation of the outflow energetics. As our observed data have missed a dedicated PSF star, we compared the spatial extent of the [O III] outflow emission with that of the BLR emission (both obtained from the previous spectral modelling). In fact, given that the latter is unresolved in our data, it is suitable for reproducing the instrumental response. We fitted a 2D-Gaussian profile to the BLR flux map, obtained by integrating in wavelength our BLR model, in order to estimate the angular resolution of our seeing-limited observations. The resulting best-fit Gaussian profiles are not circularly symmetric, especially in the case of HS 0810+2554, whose profile is elongated in the NW-SE direction. Such an elongation is mostly due to lens stretching effects and to the blending of A and B images, rather than to a possible intrinsic asymmetry of the PSF. Therefore, for both sources we assumed the angular size of the minor axis of the best-fit Gaussian profile as representative for the true PSF extent, as this is heading in the direction where the lens stretching effects are expected to be minimised. For HS 0810+2554 and SDSS J1353+1138, we estimated for the PSF a FWHM (θ_res) of 0$\stackrel{″}{.}$7 and 0$\stackrel{″}{.}$8, respectively. These NIR values are slightly smaller than the optical seeing measurements obtained with the differential image motion monitor (DIMM) during the observations, namely, $0\stackrel{″}{.}9$ and $1\stackrel{″}{.}0$, respectively (Oya et al. 2016).

To test whether the detected ionised outflows are spatially resolved, we adopted the following procedure. First, we created the flux maps for the [O III] outflow and BLR components. Then we calculated spaxel-by-spaxel the ratio of the [O III] outflow ([O III]_out) flux over the BLR flux, and produced the maps reporting the flux-ratio values across the FOV.

The ratio maps obtained for both QSOs are shown in the left panels of Fig. 4. The increasing trend in [O III]_out-to-BLR ratios with the distance from the emission peak indicates that the [O III] outflows are spatially resolved in both QSOs. Moreover, the ratio map of HS 0810+2554 highlights the existence of a preferred NE-SW direction along which the highest ratio values are found. Such a direction is perpendicular to the blending direction of the images A+B. Unlike HS 0810+2554, the two lensed images of SDSS J1353+1138 are spatially well-resolved and not affected by significant lens stretching effects. As a consequence, the [O III]_out-to-BLR ratio maps shows an isotropic pattern of increasing ratio values moving outwards from the centre of image A (we recall that we detected the [O III] emission only from this image, as previously discussed in Sect. 3.2).

Fig. 4. Spatially resolved ionised outflows in HS 0810+2554 (top) and in SDSS J1353+1138 (bottom). Left panels: maps of the ratio between the [O III] outflow flux and the BLR flux (both from best-fit models). Coloured pixels refer to a S/N ≳ 2 (S/N ≳ 3) on the full [O III] emission line (i.e., narrow + outflow components) for HS 0810+2554 (SDSS J1353+1138). The positions of the continuum emission peak of the lensed images are marked with white ‘+’, while the dotted white lines indicate the contour levels of the BLR emission at 75%, 50%, and 25% of its peak; for SDSS J1353+1138, it is represented also the level at 90%. Right panels: normalised intensity profiles along the pseudo-slit (black dotted-dashed lines in the ratio map) and in circular annuli of increasing radius for HS 0810+2554 (top) and SDSS J1353+1138 (bottom), respectively: BLR model (red lines), [O III] outflow model (dotted blue lines) and [O III] blue wing from data (cyan lines). Blue points represent ratio values of the [O III] outflow flux over the BLR flux and they are referred to the right-hand logarithmic scale. The dashed black line in the plot of SDSS J1353+1138 corresponds to the radial distance of the centre of image B.

What we discuss above, based on the use of 2D-ratio maps, can be better appreciated using spatial profiles. We determined the spatial profiles of the [O III] outflow and BLR emissions, as well as of their ratio values and studied their variability with increasing distance from the peak of the overall emission. In HS 0810+2554, since the [O III] emission is preferentially located along the NE-SW direction, we defined a pseudo-slit in such a direction (θ ∼ 130°) with a width of five spaxels, along which we calculated the emission spatial profiles. On the contrary, given the isotropic pattern of the whole emission in SDSS J1353+1138, we determined the spatial profiles in circular annuli of increasing radial distance from the centre of image A (lower panel). All the spatial profiles thus obtained are shown in the right panels of Fig. 4. Each one has been normalised to its own 0″-value, that is, to the value at the peak position of the overall emission; in the case of the [O III] outflow and BLR profiles, the 0″-value corresponds also to their own peak value. This is helpful in further confirming our previous conclusion that the [O III] outflows are spatially resolved in both galaxies since the [O III] outflow profiles are broader than the respective BLR profiles. A unique exception occurs in SDSS J1353+1138, in a correspondence of about 1$\stackrel{″}{.}$5 from the centre where we can observe a clear bump in the BLR emission profile: this is due to the flux contribution from image B. Therefore, this does not affect our previous conclusion. As a further test, in addition from that obtained from our fitting procedure, we determined the [O III] outflow profile by collapsing the spectral channels in the subtracted-data cube, including the blue wing of the [O III] line profile (4976−5000 Å and 4970−4996 Å for HS 0810+2554 and SDSS J1353+1138, respectively). The two [O III] outflow profiles, from the fit (blue dotted line) and from the spectrally-integrated (cyan solid line) subtracted-data, agree very well.

In order to estimate the angular extent of ionised outflows on the image plane³, we focused on the ratio values of the [O III] outflow flux over the BLR flux. These are plotted in logarithmic scale (on the right-hand side of the plots), after having been rescaled to 1 in the central pixel. In this way, we can easily identify values higher than 1 as regions producing a significant [O III] outflow emission and, hence, we can determine the spatial extent of resolved ionised outflows. The associated errorbars were computed by propagating the uncertainties on the [O III] and BLR fluxes in the spatial pixels involved. These were computed by propagating the error (mostly due to the noise) associated to the spectral channels, over which we integrated to get the total flux contained in that spatial pixel. We took the maximum distance, including solely ratio values not consistent with 1, as both radius within which to spatially integrate the flux of the [O III] outflow component, and outflow observed radius (R_obs) to be still corrected for the SINFONI-PSF and lensing stretching effects. To correct for the PSF-smearing, we applied the correction $R_{\mathrm{PSF}}=\sqrt{R_{\mathrm{obs}}^2-{({\theta }_{\mathrm{res}}/2)}^2}$, where R_PSF is the PSF-corrected radius, R_obs is the radius observed in the image plane and θ_res is our seeing estimate (0$\stackrel{″}{.}$7 and 0$\stackrel{″}{.}$8 for HS 0810+2554 and SDSS J1353+1138, respectively). In HS 0810+2554, we spatially integrated the [O III] outflow flux up to $R_{\mathrm{obs}}=1\stackrel{″}{.}25$, finding a total observed flux of (3.73 ± 0.05) × 10⁻¹⁵ erg s⁻¹ cm⁻² and $R_{\mathrm{PSF}}\sim 1\stackrel{″}{.}2$. In SDSS J1353+1138, we took $R_{\mathrm{obs}}=1\stackrel{″}{.}13$ and assessed a total observed flux of (8.6 ± 0.6) × 10⁻¹⁶ erg s⁻¹ cm⁻² and $R_{\mathrm{PSF}}\sim 1\stackrel{″}{.}06$⁴ for the [O III] outflow. In Sect. 4.2, we account for the lens effects and estimate the intrinsic extent (by correcting for stretching effects) and unlensed flux (by correcting for magnification effects) of the ionised outflows in HS 0810+2554 and in SDSS J1353+1138.

4. Results

4.1. Distribution and kinematics of the ionised gas

The main purpose of this work is to map the kinematics of the [O III] emission, with a primary focus on the outflow component. Figure 5 shows a global overview of the distribution and the kinematics of the ionised gas resulting from the modelling of the [O III] emission line. The moment-0 (intensity field), moment-1 (velocity field), and moment-2 (dispersion field) maps for the narrow and the outflow components are shown separately in order to better trace their distinct spatial and velocity distributions. All maps have been produced reporting only spatial pixels with a S/N equal or higher than 2 for HS 0810+2554 and higher than 3 for SDSS J1353+1138.

Fig. 5. Moment-0 (intensity field), moment-1 (velocity field), and moment-2 (dispersion field) maps of the [O III] line emission in HS 0810+2554 (left) and in SDSS J1353+1138 (right). The maps for the total, narrow and outflow components are shown separately, reporting only spatial pixels with a S/N equal or higher than 2 for HS 0810+2554, and than 3 for SDSS J1353+1138. The black ‘+’ indicates the emission centroid in each QSO, while the dotted lines represent the contour levels of the BLR emission at 75%, 50%, and 25% of its peak.

The candidate [O III]-outflow component is extended up to large distance from the galaxy centre of both QSOs, and it stands out for its strongly-blueshifted velocities and high velocity dispersion values (|v|≳500 km s⁻¹ and σ ≳ 600 km s⁻¹, respectively; see moment-1 moment-2 maps in Fig. 5 relative to the outflow component). Such velocity dispersions are well above the values measured in typical star-forming systems at these redshifts (σ ∼ 100 km s⁻¹, e.g., Cresci et al. 2009; Law et al. 2009) and, along with the overall blue-shifted motion, they provide clear evidence for large-scale outflows in these galaxies. Moreover, while in HS 0810+2554 the outflow and the narrow components have almost the same intensity across the FOV, we note that in SDSS J1353+1138, the [O III] outflow emission is brighter than the [O III] narrow emission produced by the bulk of the gas of the host galaxy.

For the narrow component, which is expected to trace mostly systemic galactic motions, we obtained low-velocity and low velocity-dispersion values (|v|≲50 km s⁻¹, σ ≲ 200 km s⁻¹ for HS 0810+2554, and |v|≲100 km s⁻¹, σ ≲ 300 km s⁻¹ for SDSS J1353+1138; see moment-1 moment-2 maps in Fig. 5 relative to the narrow component) in the central region where also the outflow emission is detected, further supporting our decomposition in the spectral analysis of both QSOs. In the more external region of HS 08101+2554, we observe slightly higher values in the velocity dispersion (σ ≲ 300 km s⁻¹), as the line profile is modelled with a unique Gaussian component, given that the [O III] emission line is fainter and the S/N is lower. This could indicate that the [O III] outflow component is still present, but cannot be isolated from the [O III] narrow component because of its faintness and the worse S/N.

Similarly to Marasco et al. (2020) and given the complexity of the [O III] line profile across the FOV, we preferred adopting the following definitions of velocity and width for the outflow characterisation (e.g., see also Harrison et al. 2014; Zakamska & Greene 2014; Cresci et al. 2015; Carniani et al. 2015; Brusa et al. 2016), rather than the moment-1 and moment-2 values. The latter are indeed more affected by geometrical projection and dust absorption effects. In each spatial pixel, we determined the 10th and 90th velocity percentiles (v₁₀ and v₉₀) of the overall emission line profile (i.e., narrow + outflow components if present), as representative velocities of the approaching and receding outflow components, respectively. The null velocity value corresponds to the systemic velocity peak of the narrow component in the central spectrum. From v₁₀ and v₉₀, we computed the line width W₈₀ defined as v₉₀ − v₁₀. The W₈₀ width is approximately equal to the full width at half maximum (FWHM) for a Gaussian profile. Maps of v₁₀, v₉₀ and W₈₀ are shown in Fig. 6.

Fig. 6. v10, v90, and W80 maps of the [O III] emission line in HS 0810+2554 (upper panels) and in SDSS J1353+1138 (lower panels). We applied the same cut in S/N as in the moment maps of Fig. 5, that is, with a S/N equal or higher than 2 and for HS 0810+2554 and higher than 3 for SDSS J1353+1138.

The maps of v₁₀ show highly blueshfited velocities in most of the field of HS 0810+2554 and SDSS J1353+1138. In the former, a slightly steeper velocity gradient is present in the west direction from the centre, where we observe velocities as high as about −2170 km s⁻¹; in the latter, the outflow region is preferentially elongated in the NE-SW direction with highest velocity values (up to −2410 km s⁻¹) at the SW end of the strongly blueshifted region. In HS 0810+2554, we clearly detect also the redshifted component of the outflow in the two reddest regions in the v₉₀ map, where the outflow is seen receding from us at velocities up to about 1730 km s⁻¹ along the line of sight. Looking at the W₈₀ maps, we observe extremely large values (1100 km s⁻¹ ≲ v ≲ 3500 km s⁻¹) in the outflow regions, which is consistent with other z ∼ 2 QSO outflows found in the literature (Carniani et al. 2015). Furthermore, by comparing the [O III]-outflow moment-1 map with the v₁₀ map for each QSO, we note that the shape of the [O III]-outflow moment-1 map reflects the bluest region in the v₁₀ map, suggesting that any additional Gaussian component added to model the wings in the [O III] profile has been correctly classified as outflow component (compare also [O III]-outflow moment-2 maps and the respective W₈₀ maps).

We rule out the possibility of alternative scenarios, such as galactic inflows or a galaxy merger event. In fact, in the few reported cases of their detection, galactic inflows have been observed mostly in absorption and with quite small bulk velocities (∼200 km s⁻¹) and velocity dispersions (e.g., Bouché et al. 2013). Moreover, for the inflowing gas theoretical modelling predicts a small covering factor (e.g., Steidel et al. 2010), making its direct observation rare especially at high redshift (e.g., Cresci et al. 2010). Finally, we exclude also the galaxy-merger scenario since the deep optical images of both HS 0810+2554 and SDSS J1353+1138 (see Fig. 1) do not show any continuum emission counterpart in correspondence of the outflow region, which could support such a scenario.

Finally, we stress that the obtained maps are relative to the lens plane and, thus, they do not account for gravitational lensing effects. While these are expected not to significantly affect the observed gas kinematics (hence the outflow velocity), they strongly alter the observed gas spatial distribution and the observed surface brightness: fluxes are magnified and spatial dimensions are stretched. Therefore, the obtained maps could not be used to infer directly the outflow intrinsic radius and its total flux, which are key ingredients, along with the outflow velocity, in the computation of the outflow energetics. We first need to quantify the lensing effects and then we can derive the unlensed physical properties of the outflow. This aspect is discussed in Sect. 4.2 (see also Appendix A).

4.2. Inferring unlensed size and flux of the outflow

As we discuss in Sect. 4.1, we performed both the spectral analysis and the kinematical study in the lens plane. For this reason, we had to correct for the lens effects to determine the actual extent and flux of the outflow.

There are several adaptive-mesh fitting codes which, given a mass distribution for the lens and a surface brightness profile for the background source, fit the resulting forward lensed image to the observed data and use a statistics test (e.g., the minimum χ² method) to establish the best-fit models for both the lens and the source. These algorithms usually require the knowledge of the instrumental PSF to allow a correct comparison with the observed data. The output of these fitting-codes is a 2D or 3D reconstruction (depending on the code used) of the unlensed source. In order to achieve an accurate reconstruction, it is required that the lensed images are all detected and spatially resolved⁵, as their position depends on the first derivative of the gravitational potential of the lens, while their flux on the second derivative (e.g., Jackson et al. 2015; Nierenberg et al. 2020). In other words, the knowledge of the position of the multiple lensed images and of their flux provides strong constraints on the lens and background source models.

Unfortunately, we could not use such fitting codes to fully reconstruct the unlensed outflow in the source plane for HS 0810+2554 nor SDSS J1353+1138 since our data did not satisfy the necessary requirements. In fact, in the case of HS 0810+2554 the spatial resolution of the SINFONI data was too low to resolve the various lensed images and, thus, we were not able to achieve an accurate full reconstruction of the background source. However, we were able to recover a partial reconstruction of the background outflow using the lens-fitting code from Rizzo et al. (2018) and adopting an approximated procedure (see Appendix A). In this way, we estimated the outflow magnification factor and intrinsic radius to be, respectively, μ_out = 2.0 ± 0.2 and R_out = (8.7 ± 1.7) kpc, with z = 1.508 ± 0.002 as the redshift measured from the nuclear spectrum extracted for the BLR-modelling (described in Sect. 3.1.1) and, here, adopted to convert the outflow angular size into kpc units. Our z measurement is consistent with the centroids of the ALMA CO(J = 3 → 2) and CO(J = 2 → 1) emission lines of HS 0810+2554 reported in Chartas et al. (2020). By correcting for μ_out the observed [O III] outflow flux (determined in Sect. 3.3), we found the unlensed outflow flux to be F_out = (1.9 ± 0.2) × 10⁻¹⁵(2.0/μ_out) erg s⁻¹ cm⁻². Our μ_out ∼ 2 estimate is close to what Chartas et al. (2020) found for a high-velocity CO-clump at similar distance in ALMA data of HS 0810+2554. On the contrary, it remarkably differs (up to two orders of magnitude) from the values from the literature, presented in Sect. 2.1.1. This follows from the fact that the latter are estimates of the magnification of the emission from the more compact (UV disk and X-ray corona) region, while our estimate is relative to a large-scale (∼8 kpc) emission. Moreover, the reconstructed unlensed outflow does not intercept the lens caustics (see Appendix A) and, hence, it misses the magnification contribution from those regions, where the lens magnification is drastically larger.

In SDSS J1353+1138, the complete lack of [O III] detection in image B prevented us from attempting any background source reconstruction starting from our data, as no constraints could be put on this image. Therefore, we had to make simplified assumptions, referring to the lens models for SDSS J1353+1138 by Inada et al. (2006) and Rusu et al. (2016) (discussed in Sect. 2.3), based on AGN plus host galaxy emission in the i and K-band, respectively. Our assumptions rely on the fact that in SDSS J1353+1138 gravitational lensing effects are supposed to be smaller than in HS 0810+1154 (e.g., Narayan & Bartelmann 1996). As a consequence, as compared to HS 0810+1154, for this object we expect: (1) a milder and almost isotropic stretching of physical dimensions, as we indeed observed; (2) smaller and less variable-in-space values of differential magnification. On the basis of the first argument, we neglected the stretching lens effect and approximated the unlensed outflow angular size to $R_{\mathrm{PSF}}=1\stackrel{″}{.}06\pm 0\stackrel{″}{.}13$ (determined in Sect. 3.3). Considering instead the lens magnification, we expect total magnification factors of a few units that is weakly dependent on the geometrical details of the flux distribution for background emissions with comparable spatial extent. Consequently, we used the average between the i-band (i.e., μ = 3.81 and μ = 3.75; Inada et al. 2006) and K-band total magnification factors (i.e., μ = 3.47, μ = 3.42 and μ = 3.53; Rusu et al. 2016) as a proxy for the total outflow magnification μ_out, under the assumption of comparable unlensed physical sizes. Given the unknown real unlensed flux distribution of the J-band outflow, we conservatively assumed an uncertainty of 10% on our adopted μ_out value, thus obtaining μ_out = 3.6 ± 0.4. Correcting, finally, for the lens magnification and converting to kpc-units, we found the unlensed flux for the outflow to be F_out = (2.4 ± 0.3) × 10⁻¹⁶(3.6/μ_out) erg s⁻¹ cm⁻² and its intrinsic radius to be R_out = (9.2 ± 1.1) kpc, using z = 1.632 ± 0.002 as measured from the nuclear spectra extracted during the BLR-modelling (described in Sect. 3.2).

In Table 2, we summarise the main outflow properties for HS 0810+2554 and SDSS J1353+1138, obtained up to this point. The first columns show the maximum velocity values observed in the v₁₀, v₉₀ and W₈₀ maps (described in Sect. 4.1), respectively, referred to as $v_{10}^{\max }$, $v_{90}^{\max }$ and $W_{80}^{\max }$. Then we report our lens-corrected estimates of R_out and F_out (inferred as discussed above), the latter corrected for the outflow magnification factors, μ_out, shown in the last column. Most of these physical quantities are also employed in the computation of the outflow energetics in Sect. 4.3.

4.3. Outflow energetics

We derived the physical properties of the large-scale ionised outflows in HS 0810+2554 and SDSS J1353+1138 from the observed [O III]λ5007 emission, following the prescriptions described in Cano-Díaz et al. (2012) as done also in Marasco et al. (2020). The [O III] line luminosity is given by:

$$\begin{array}{c}\hfill L_{[\mathrm{O}\mathrm{III}]}={\displaystyle {\int }_{\mathrm{V}}{ϵ}_{[\mathrm{O}\mathrm{III}]}f\mathrm{d}V,}\end{array}$$(2)

where V is the volume occupied by the ionised outflow, f is the filling factor of the [O III] emitting clouds in the outflow, and ϵ_[O III] is the [O III] line emissivity which, at the temperatures typical of the NLR (∼10⁴ K), is weakly dependent on the temperature (∝T^0.1) and can be written as:

$$\begin{array}{c}\hfill {ϵ}_{[\mathrm{O}\mathrm{III}]}=1.11\times {10}^{-9}{E}_{[\mathrm{O}\mathrm{III}]}{n}_{{\mathrm{O}}^{2+}}{n}_{\mathrm{e}}{\mathrm{erg}\mathrm{s}}^{-1}{\mathrm{cm}}^{-3},\end{array}$$(3)

with E_[O III] as the energy of the [O III] photons, n_O²⁺ and n_e, the volume densities of the O²⁺ ions, and of the electrons, respectively. Then assuming that most of the oxygen in the ionised outflow is in the form of O²⁺, it follows that:

$$\begin{array}{c}\hfill {ϵ}_{[\mathrm{O}\mathrm{III}]}\approx 5\times {10}^{-13}{E}_{[\mathrm{O}\mathrm{III}]}{n}_{\mathrm{e}}^2{10}^{[\mathrm{O}/\mathrm{H}]}{\mathrm{erg}\mathrm{s}}^{-1}{\mathrm{cm}}^{-3},\end{array}$$(4)

where [O/H] gives the oxygen abundance in solar units. The mass of the outflowing ionised gas can be derived from the following expression:

$$\begin{array}{c}\hfill M_{\mathrm{out}}={\displaystyle {\int }_{V}1.27m_{\mathrm{H}}{n}_{\mathrm{e}}f\mathrm{d}V,}\end{array}$$(5)

where m_H is the mass of the hydrogen atom and the factor of 1.27 follows from including the mass contribution of helium. By combining Eqs. (2) and (5), we get:

$$\begin{array}{c}\hfill M_{\mathrm{out}}=5.33\times {10}^7(\frac{L_{[\mathrm{O}\mathrm{III}]}}{{10}^{44}{\mathrm{erg}\mathrm{s}}^{-1}})(\frac{⟨n_{\mathrm{e}}⟩}{{10}^3{\mathrm{cm}}^{-3}}{)}^{-1}\frac{C}{{10}^{[\mathrm{O}/\mathrm{H}]}}{M}_{\odot },\end{array}$$(6)

where ⟨n_e⟩ is the electron density averaged over the ionised outflow volume (i.e., $⟨n_{\mathrm{e}}⟩={\int }_V{n}_{\mathrm{e}}f\mathrm{d}V/{\int }_{V}f\mathrm{d}V$) and $C={\langle n_{\text{e}}\rangle }^2/\langle n_{\text{e}}^2\rangle$ is the so-called ‘condensation factor’. Under the simplifying hypothesis that all ionising gas clouds have the same density, we get C = 1 and eliminate the outflow mass dependance on the filling factor of the emitting clouds.

In order to compute the energetics of the ionised outflow, we have to make further simplifying assumptions for the outflow geometry and structure: we assume that the outflow has a (bi-)conical geometry with an opening angle Ω and a radial extent, R_out, and that it consists of a collection of ionised gas clouds, uniformly distributed within the cone and outflowing with a speed v_out. The mass outflow rate is given by:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{out}}=⟨\rho ⟩v_{\mathrm{out}}\mathrm{\Omega }{R}^2,\end{array}$$(7)

where ⟨ρ⟩ is the average mass density computed over the total volume V occupied by the conical outflow⁶. By substituting ⟨ρ⟩ in Eq. (7) with its definition in terms of Ṁ_out (using Eq. (6)) and V, we obtain:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{out}}=164(\frac{L_{[\mathrm{O}\mathrm{III}]}}{{10}^{44}{\mathrm{erg}\mathrm{s}}^{-1}})(\frac{⟨n_{\mathrm{e}}⟩}{{10}^3{\mathrm{cm}}^{-3}}{)}^{-1}(\frac{v_{\mathrm{out}}}{{10}^3{\mathrm{km}\mathrm{s}}^{-1}})(\frac{R_{\mathrm{out}}}{\text{kpc}}{)}^{-1}\frac{M_{\odot }{\mathrm{yr}}^{-1}}{{10}^{[\mathrm{O}/\mathrm{H}]}}\end{array}$$(8)

where we have assumed C = 1. The mass outflow rate thus inferred does not depend on either the opening angle Ω of the outflow cone or the filling factor f of the emitting clouds.

Finally, we calculate the kinetic energy (E_kin), kinetic luminosity (L_kin) and momentum rate (ṗ_out) of the outflow by means of the following expressions:

$$\begin{array}{cc}& E_{\mathrm{kin}}=9.94\times {10}^{42}(\frac{M_{\mathrm{out}}}{M_{\odot }})(\frac{v_{\mathrm{out}}}{{\mathrm{km}\mathrm{s}}^{-1}}{)}^2\mathrm{erg},\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill & L_{\mathrm{kin}}=3.16\times {10}^{35}(\frac{{\dot{M}}_{\mathrm{out}}}{M_{\odot }{\mathrm{yr}}^{-1}})(\frac{v_{\mathrm{out}}}{{\mathrm{km}\mathrm{s}}^{-1}}{)}^2{\mathrm{erg}\mathrm{s}}^{-1},\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill & {\dot{p}}_{\mathrm{out}}=6.32\times {10}^{30}(\frac{{\dot{M}}_{\mathrm{out}}}{M_{\odot }{\mathrm{yr}}^{-1}})(\frac{v_{\mathrm{out}}}{{\mathrm{km}\mathrm{s}}^{-1}})\text{dyne}.\hfill \end{array}$$(11)

Equations (6)–(11) require the knowledge of different physical properties of the outflow, some of which we were able to derive, while others had to be assumed. These are: the oxygen abundance, which we fixed to the solar abundance, and the electron density, which we assumed to be n_e ∼ 1000 cm⁻³, in agreement with the values measured in similar studies at high redshift (see e.g., Perna et al. 2017; Förster Schreiber et al. 2019). The second assumption, in particular, affects the derived outflow energetics (Davies et al. 2020; Kakkad et al. 2020) but it is necessary, nonetheless, since we cannot measure n_e directly from our data.

We now focus on the physical quantities we were able to calculate. In Sect. 4.2, we provided the values of the intrinsic radius R_out and flux F_out of the ionised outflows, traced by the [O III] line emission. From F_out we calculated the intrinsic [O III] line luminosity used in the mass outflow expression (Eq. (6)). We found (μ_out-corrected) [O III] luminosity values of L_[O III] = (2.8 ± 0.3) × 10⁴³ erg s⁻¹ and L_[O III] = (4.4 ± 0.6) × 10⁴² erg s⁻¹ for HS 0810+2554 and SDSS J1353+1138, respectively.

In order to establish the velocity of the ionised outflows, we focused on the kinematical maps of v₁₀ and v₉₀, shown in Sect. 4.1. The spectral analysis and the study of the gas kinematics in HS 0810+2554 and in SDSS J1353+1138 have revealed the presence of an extended central region hosting outflowing gas directed towards us along the line of sight (Sect. 4.1), leading to very high values of v₁₀. Moreover, in HS 0810+2554, we detected also a red wing in the [O III] line profile in two smaller regions apart from the peak of the overall emission, corresponding to the high-velocity receding component of the ionised outflow. In this case, following Marasco et al. (2020), we defined the outflow velocity as:

$$\begin{array}{c}\hfill v_{\mathrm{out}}=\max (|v_{10}^{\max }-v_{\mathrm{sys}}|,|v_{90}^{\max }-v_{\mathrm{sys}}|),\end{array}$$(12)

where $v_{10}^{\max }$ and $v_{90}^{\max }$ are the maximum value, respectively, observed in the v₁₀ and v₉₀ maps (described in Sect. 4.1), and v_sys is the bulk (or systemic) velocity of the galaxy, inferred from the nuclear spectrum used for the BLR-fitting (in Sect. 3.1.1) and set to the value of 0 km s⁻¹, as previously described in Sect. 4.1. This definition is required by the unknown geometry and orientation of the outflow with respect to the line of sight: since we ignore the true angle of the outflow with respect to the line of sight, and given that the bulk of the outflow unlikely points towards the observer, we assume that the best representation of the outflow speed is provided by the velocity ‘tail’ of the line profile, that is, v₁₀ and v₉₀ in Eq. (12). These values are thought to be more suited to represent v_out than the mean (or median) velocity of the line, which strongly depends on projection effects and dust absorption (e.g., Cresci et al. 2015).

The same argument holds also for the determination of the outflow velocity in SDSS J1353+1138. But in this case, possibly because of the different orientation with respect to the observer and higher dust absorption, we did not observe any asymmetric red wing in the [O III] profile produced by the receding part of the outflow, as stressed in Sect. 4.1. Therefore, for SDSS J1353+1138 we focused only on the blue tail of the [O III] line associated to the outflow, and hence assumed $v_{\mathrm{out}}=v_{10}^{\max }$ as outflow velocity for SDSS J1353+1138.

To calculate the quantities in Eqs. (6)–(11) and their uncertainties, we used the error propagation considering the inferred errors on R_out, L_[O III] and v_out, and a typical uncertainty of 50% on n_e (e.g., Perna et al. 2017; Förster Schreiber et al. 2019). The physical properties of the ionised outflows detected in HS 0810+2554 and in SDSS J1353+1138 are reported in Table 3, including also our estimates of the kinetic efficiency and of the momentum-boost. The former is defined as the ratio between the outflow kinetic luminosity, L_kin, (defined in Eq. (10)) and the AGN bolometric luminosity, L_Bol, (corrected for the lens magnification), while the latter is defined as the ratio between the momentum rate of the outflow (ṗ_out) and the momentum initially provided by the AGN-radiation pressure (i.e., L_Bol/c), which is approximately identified also with the momentum rate of the X-ray UFO. The values of L_Bol adopted in this work are for (2.5 ± 0.9) × 10⁴⁵ erg s⁻¹ and (39 ± 2) × 10⁴⁵ erg s⁻¹ HS 0810+2554 and SDSS J1353+1138, respectively, and they will be discussed in Sect. 5.1.

We found mass outflow rates of ∼2 M_⊙ yr⁻¹ and ∼12 M_⊙ yr⁻¹, and kinetic efficiencies of ∼9 × 10⁻⁵ and ∼700 × 10⁻⁵ for SDSS J1353+1138 and HS 0810+2554, respectively. The values obtained for HS 0810+2554 are in good agreement with the predictions at L_Bol ∼ 2 × 10⁴⁵ erg s⁻¹ of the Ṁ_out − L_Bol and L_kin − L_Bol scaling relations (Carniani et al. 2015; Fiore et al. 2017) for the ionised outflow component. On the contrary, the inferred Ṁ_out value for SDSS J1353+1138 is small (by a factor of ∼100) compared to the predictions for the ionised outflow component at L_Bol ∼ 4 × 10⁴⁶ erg s⁻¹. Our findings are however relative only to the ionised phase of large-scale outflows, while a significant amount of outflowing gas may be in neutral molecular and/or atomic phase. Given the typical bolometric luminosities of our galaxies (L_Bol ∼ 10⁴⁵ − 10⁴⁶ erg s⁻¹), outflow mass rates predicted for the molecular component may indeed exceed our measurements for the ionised gas by a factor ∼100. Chartas et al. (2020) have recently claimed the tentative detection of a massive (M_{out,  mol} ∼ 4 × 10⁹ M_⊙) CO-molecular outflow in HS 0810+2554, with total mass rate and velocity of ∼400 M_⊙ yr⁻¹ and 1040 km s⁻¹, respectively. Compared to molecular outflows, which have been observed in z > 1 QSOs in spite of the uncertainty due to detection limits (e.g., Carniani et al. 2017; Feruglio et al. 2017; Brusa et al. 2018), neutral atomic outflows mainly reside in ultra-luminous infrared galaxies (ULIRGs), showing both intense SF and AGN activity (e.g., Rupke et al. 2005; Cazzoli et al. 2016; Perna et al. 2017; Fluetsch et al. 2019, 2020). This may suggest that neutral atomic outflows are mostly powered by SF rather than by AGN activity (e.g., Concas et al. 2017; Bae & Woo 2018) or that they can occur in obscured AGN hosting large quantities of cold gas to be channeled into galaxy-scale outflows (Perna et al. 2017). However, the contribution of any additional gas phase may significantly increase the overall outflow mass and energetics (e.g., Cicone et al. 2014; Carniani et al. 2016; Fiore et al. 2017). The implications of all these results will be discussed in details in Sect. 5.

5. Discussion

5.1. Connection with the nuclear X-ray UFOs

The main purpose of this work is to shed light on the acceleration mechanism of ionised outflows on large scales. In this regard, we compared the energetics of the galaxy-scale ionised outflows to the UFOs present at nuclear scales, in order to test whether they are causally connected (i.e., whether they are subsequent phases of the same AGN-accretion burst). In the case of a causal connection, we could go on to investigate the nature of the acceleration mechanism distinguishing between momentum-driven and energy-driven winds.

We show the energetics measurements of the large scale ionised outflows in HS 0810+2554 and SDSS J1353+1138 in Sect. 4.3. In this section, we present the X-ray measurements of the hosted UFOs from the literature that are used to determine the sub-pc scale wind energetics. For both QSOs, we refer to Chartas et al. (2021), who present the first X-ray spectral analysis of SDSS J1353+1138, and new UFOs measurements for HS 0810+2554 obtained from a new Chandra observation acquired in 2016 using the updated version of the photoionisation code XSTAR (Bautista & Kallman 2001).

As done in Marasco et al. (2020), we followed Nardini & Zubovas (2018) in order to re-compute the UFOs energetics in a consistent way based on the same assumption we made in the calculation of the large-scale outflow energetics. We assumed that the UFO is launched from the escape radius $r_{\text{esc}}\equiv 2G{M}_{\text{BH}}/v_{\text{UFO}}^2$ of the BH and we derived the mass outflow rate for the nuclear wind as:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{UFO}}\simeq 0.3(\frac{\mathrm{\Omega }}{4\pi })(\frac{N_{\mathrm{H}}}{{10}^{24}{\mathrm{cm}}^{-2}})(\frac{M_{\mathrm{BH}}}{{10}^8{M}_{\odot }{\mathrm{yr}}^{-1}})(\frac{v_{\mathrm{UFO}}}{c}{)}^{-1}{M}_{\odot }{\mathrm{yr}}^{-1},\end{array}$$(13)

where Ω is the solid angle subtended by the UFO, M_BH the black hole mass, N_H the hydrogen gas column density, and v_UFO the wind speed. We took the values of v_UFO and N_H inferred by Chartas et al. (2021). For HS 0810+2554, these values (reported in Table 4) slightly differ from those previously published (i.e., $v_{\mathrm{UFO},1}=0.{12}_{-0.01}^{+0.02}c$, $N_{\mathrm{H},1}=3.4_{-2.0}^{+1.9}\times {10}^{23}{\mathrm{cm}}^{-3}$ and $v_{\mathrm{UFO},2}=0.{41}_{-0.04}^{+0.07}c$, $N_{\mathrm{H},2}=2.9_{-1.6}^{+2.0}\times {10}^{23}{\mathrm{cm}}^{-3}$ for the two distinct UFO components; Chartas et al. 2016), but are still in agreement. The adopted M_BH values are virial estimates based on Hβ for HS 0810+2554 (Assef et al. 2011) and on C IV for SDSS J1353+1138 (Chartas et al. 2021). The C IV-based measurement of M_BH for SDSS J1353+1138 has been corrected following the prescription published in Coatman et al. (2017) for C IV-based virial black hole mass estimates, which are known to be systematically biased compared to masses derived from the Balmer hydrogen lines. Finally, we assumed for both UFOs a covering factor of $f=\frac{\mathrm{\Omega }}{4\pi }=0.4_{-0.2}^{+0.2}$, on the basis that about 40% of local AGNs have been observed to host UFOs (Tombesi et al. 2010; Gofford et al. 2013). With all these ingredients, we calculated the mass rate Ṁ_UFO (with Eq. (13)), the momentum rate ṗ_UFO (i.e., ṗ_UFO = Ṁ_UFOv_UFO) and the momentum-boost of the nuclear winds defined as ṗ_UFO/(L_Bol/c). As uncertainty on all quantities, the minimum–maximum range of possible values was considered, considering the values of v_UFO, N_H and M_BH as inferred in Chartas et al. (2021). For L_Bol we took the average between the μ-corrected values obtained in Chartas et al. (2021) with two different methods: the 2–10 keV bolometric correction (Lusso et al. 2012) and the estimate from the continuum luminosity at 1450 Å (Assef et al. 2011; Runnoe et al. 2012).

Table 4 summarises the physical properties of the UFOs in HS 0810+2554 and SDSS J1353+1138, derived following the prescriptions in Nardini & Zubovas (2018), along with the physical quantities taken from Chartas et al. (2021). For HS 0810+2554 the values of Ṁ_UFO and ṗ_UFO/(L_Bol/c) refer to the whole hosted UFO, as sum of the two detected UFO components at different speeds (for which we separately report N_H and v_UFO in Table 4), originally discovered in Chartas et al. (2016) and confirmed in Chartas et al. (2021).

Figure 7 shows the momentum-boost as a function of the outflow speed for HS 0810+2554 and SDSS J1353+1138. Given the two X-ray UFOs points (represented as circles), the predictions for a momentum-driven and an energy-driven scenario are represented by the dashed horizontal and diagonal lines, respectively. Errorbars on UFOs points indicate the minimum-maximum range of possible values, as mentioned previously. In comparing with ionised outflow points (represented as squares), we observe that while the [O III] outflow of HS 0810+2554 is consistent (within the uncertainty) with the expectations for the momentum-driven mode, that of SDSS J1353+1138 is smaller by a factor of ∼100 than the predictions for a momentum-driven propagation, therefore, it is far from any tentative connection with the UFOs on nuclear scale. We reiterate that our measurement is a lower limit, as we did not detect [O III] emission from the second lensed image and, moreover, we approximated the outflow intrinsic radius to the one observed (just PSF-corrected). Nevertheless, even accounting for all these issues and approximations we made, such a discrepancy between observations and theoretical predictions could hardly be explained, as it amounts to about two orders of magnitude. Even associating half of the [O III] flux observed in image A to image B, on the basis of the flux ratio measurement of the two H-band images (Inada et al. 2006), we would obtain a momentum-boost larger than the previous one by only a factor ∼1.5. Similarly, it is unlikely that we are underestimating the stretching effects so much as to overestimate R_out by a factor of ∼100. Therefore, we conclude that such a significant discrepancy must have a different origin. The more plausible hypotheses are: (1) either the likely presence of a massive molecular outflow in this galaxy that our work is not accounting for; or (2) the possibility that the observed UFO is caused by an extraordinary burst episode (see e.g., Zubovas & Nardini 2020) in the BH accretion activity of SDSS J1353+1138, while the large-scale outflow must be considered as the resultant effect of the AGN activity averaged over longer time-scales (Woo et al. 2017). For this object, Chartas et al. (2021) estimate a photon index Γ ∼ 2.2, which is typical of Narrow-line Seyfert 1 galaxies (Leighly 1999; Vaughan et al. 1999) and a typical signature of highly accreting systems (Huang et al. 2019). The Eddington ratio λ_Edd (defined as L_Bol/L_Edd) estimated for SDSS J1353+1138 is λ_Edd = 0.20 ± 0.02 (Chartas et al. 2021), which is larger by a factor ∼3 than that inferred for HS 0810+2554 (λ_Edd = 0.07 ± 0.03; Chartas et al. 2021). Such results could support the recent post-burst scenario. Certainly, each hypothesis does not necessarily exclude the other and the real situation can be a combination of the two (effects).

Fig. 7. Momentum-boost versus wind velocity diagram for both UFOs (represented as circles) and galactic outflow components (represented as squares) of HS 0810+2554 and SDSS J1353+1138. The dashed lines represent the predictions for a momentum-conserving regime (horizontal) and energy-conserving mode (diagonal). For HS 0810+2554, also the point relative to the ionised+molecular large-scale outflow is shown, identified by the star symbol.

The small values of momentum-boost (∼0.02−2) and of kinetic efficiency (∼9−700 × 10⁻⁵), inferred for the ionised outflows in SDSS J1353+1138 and HS 0810+2554, could be explained by an overall scarcity of [O III] (not only for the outflow component) in highly accreting AGNs, observed in other local galaxies with similar properties (e.g., Sulentic et al. 2000). In our two QSOs, the poor outflow energetics is indeed due mainly to small values of outflow mass, and not to particularly low velocities. Such a scenario is related to the Eigenvector-1 effect: while the sources accreting at high rates (close to the Eddington limit) are actually the most promising candidates for hosting an active UFO (e.g., Nardini et al. 2019), they usually present a very bright Fe II emission and a faint, outshined [O III] emission. As a consequence, the [O III] may be not ideal to trace the ionised phase of the outflow in AGNs accreting at high rates since the bulk of the ionised gas could be in the form of different chemical species.

It is also possible we underestimated the uncertainty of the [O III] outflow in HS 0810+2554, mostly due to our approximated procedure in the unlensed reconstruction. Moreover, given the recent tentative detection of a more massive CO-outflow on large scales claimed in Chartas et al. (2020), such a small value of the momentum-boost of the ionised outflow is not as unexpected given that it accounts only for the ionised gas traced by the [O III] emission. Hence, we also determined the momentum-boost relative to the ionised plus molecular large-scale outflow (p_{out,  tot}/(L_Bol/c) ∼ 38), assuming its velocity to be equal to the mass-weighted average between the ionised and the CO-molecular outflow velocities (v_out, tot ∼ 1044 km s⁻¹). Given the two orders of magnitude of difference between ionised and molecular outflow masses, the mass-weighted velocity is essentially the molecular outflow velocity (v_{out,  tot} = 1040 km s⁻¹; Chartas et al. 2020).

In Fig. 7, we report the CO+[O III] point with its uncertainty. Once the contribution of the molecular component is included, the energetics of the overall large-scale outflow in HS 0810+2554 is compatible with an energy-driven scenario of wind propagation, within the (large) UFO uncertainty. However, deeper observations are required to confirm the CO-outflow detection and to constrain its energetics.

5.2. Comparison with other QSOs hosting UFOs

Our study has revealed that the energetics of the galaxy-scale ionised outflow in HS 0810+2554 is consistent with the expectations of the momentum-driven mechanism, whereas in SDSS J1353+1138, the ionised outflow does not seem to be related with the detected UFO event on sub-pc scales. However, the exiguity of our two-QSOs sample prevents us from testing the predictions of theoretical models, as well as making general considerations on the nature of the mechanism powering outflows on large scales. Therefore, we considered our results along with those of a sample of well-studied QSOs, hosting both UFOs and galactic outflows, recently collected in Marasco et al. (2020). The QSO-sample consists of the two local objects, MR225−178 and PG 1126−041, analysed in Marasco et al. (2020) to trace the ionised phase of the large-scale outflows, similarly to this work, but at low redshift; in addition, the other local QSOs having reliable UFOs and molecular or atomic outflow measurements. In terms of redshift, the only exception is the lensed QSO APM 08279+5255 at z ∼ 3.9 (Hasinger et al. 2002), which has been found to host a molecular outflow in energy-driven regime (Feruglio et al. 2017). In addition to HS 0810+2554 and SDSS J1353+1138, this is the only other case at high redshift, for which it has been possible to study the connection between the nuclear and the galaxy-scale winds. This highlights once again the importance of gravitational lensing as a powerful tool to overcome the limits imposed by the current observations.

Marasco et al. (2020) re-computed the UFO mass rate (and the wind energetics, consequently) for each QSO of their gathered sample, starting from the known estimates for M_BH, N_H, and v_UFO. All measurements for molecular outflows have been re-scaled to the same luminosity-to-mass conversion factor α_CO = 0.8 (K km s⁻¹ pc²)⁻¹ M_⊙, typical of QSOs, starburst and submillimeter galaxies (Downes & Solomon 1998; Bolatto et al. 2013; Carilli & Walter 2013). The main AGN, X-ray-wind, and large-scale outflow properties of the known QSO-sample are listed in Table B.1 in Marasco et al. (2020).

Figure 8 represents the updated version of Fig. 9 in Marasco et al. (2020), including the measurements relative to HS 0810+2554 and SDSS J1353+1138: it shows the ratios between the galaxy-scale and sub-pc scale outflow momentum rates (ṗ_out and ṗ_UFO, respectively), ordered according to the increasing L_Bol of the host AGN. This is an alternative way (with regard to Fig. 7) to compare observational results with theoretical predictions for each QSO and to distinguish between the two main regimes. The energetics of the large-scale outflow in 10 out of 12 QSOs results to be consistent (or nearly consistent) with either a momentum-driven or an energy-driven regime: in seven (six) and three (four) objects, respectively, if we exclude (include) the contribution of the tentatively detected molecular outflow in HS 0810+2554. This globally confirms the conclusion drawn by Marasco et al. (2020) that models of either momentum- or energy-driven outflows describe the mechanism of wind propagation on galaxy-scales very well. The only two exceptions are SDSS J1353+1138 and IRAS 17020+4544, representative of the two extreme opposite situations in which the energetics on large and nuclear scales seem to be completely unrelated. As discussed in Sect. 5.1, the low value of ṗ_out/ṗ_UFO in SDSS J1353+1138 could be due to the presence of a massive molecular outflow, which our work does not account for. Alternatively, it could indicate that SDSS J1353+1138 has recently undergone a burst episode of its AGN activity. In terms of high AGN variability, the high value of ṗ_out/ṗ_UFO inferred for IRAS 17020+4554 could be explained by invoking a past higher BH accretion activity compared to present day, but observations reveal that IRAS 17020+2554 is now accreting at a substantial fraction of its Eddington rate (λ_Edd ∼ 0.7; Longinotti et al. 2015). Finally, we do not observe any remarkable trend in ṗ_out/ṗ_UFO values with L_Bol, nor any evident dependence of wind acceleration mechanism on galaxy redshift upon separately inspecting the results obtained for the high-redshift (our two objects plus APM 08279+5255) and low-redshift QSOs (the remaining ones) in the sample.

Fig. 8. Ratio between the galaxy-scale and sub-pc scale outflow momentum rates for different QSOs hosting UFOs. Measurements for individual objects are shown in blue with the respective errorbars, using different markers according to the gas phase of the observed large scale outflow. The galaxy points are ordered by increasing LBol. The horizontal dashed line shows the prediction for a momentum-driven wind (ṗout/ṗUFO = 1), while the orange rectangles indicate the individual predictions for energy-driven winds. Filled and empty squares represent ionised outflow measurements based on Hα and [O III] emission, respectively. For HS 0810+2554, our [O III]-based measurement is shown both alone and combined with the CO-measurement of the molecular outflow (Chartas et al. 2020), with the respective symbol and combination of symbols (see the plot legend).

6. Conclusions

Galaxy-wide outflows powered by AGN activity are thought to play a fundamental role in shaping the evolution of galaxies, as they allow us to reconcile theoretical models to observations. However, even though observations have widely confirmed their presence in both local and high-redshift galaxies, a clear understanding of the mechanism which accelerates these powerful, galaxy-scale winds is still lacking. To test the predictions of the current theoretical models, we need to compare, in a given object, the energetics of the sub-pc wind with that of the galaxy-wide outflow. The optimal sources in attempting the make a connection between different scales are the powerful QSOs near the peak of AGN activity (z ∼ 2), where the AGN feedback is expected to be more effective.

Given such a perspective, this work focuses on two z ∼ 1.5 multiple lensed QSOs, specifically selected to host UFOs (HS 0810+2554 and SDSS J1353+1138) and observed with the near-IR integral field spectrometer SINFONI. Thanks to the strong lens magnification and to the spatially resolved SINFONI data, which trace the dynamics of the ionised gas phase through rest-frame optical emission lines, it has become possible, for the first time, to attempt to make the connection between the sub-pc winds and the large-scale ionised outflows in two QSOs near the peak of AGN activity. The only other well-studied case at high redshift is APM 08279+5255 (z ∼ 3.9), which has been suggested to host a molecular outflow in energy-conserving regime (Feruglio et al. 2017). Moreover, the recent analysis of ALMA data of HS 0810+2554 has revealed the tentative detection of CO-molecular outflow (Chartas et al. 2020). Therefore, the characterisation of the ionised phase of the outflow in HS 0810+2554 provides the first three-phase description of an AGN-driven wind at high redshift, from the X-ray to the optical and mm bands, corresponding to highly ionised, ionised, and molecular gas phases, respectively.

In the following, we summarise the main results obtained in this work:

1. We studied the gas kinematics to identify the presence of outflowing gas by tracing the emission of the forbidden line doublet [O III]λλ4959,5007, whose line profile is highly asymmetric in presence of outflows, with a typical broad, blueshifted wing corresponding to high speeds along the line of sight. In both QSO spectra, we detected the presence of extended (∼8 kpc) ionised outflows moving up to v ∼ 2000 km s⁻¹ in the image lens plane.

2. After correcting for the gravitational lensing effects, we found that the ionised outflow in HS 0810+2554 is consistent within the uncertainty with the predictions for a momentum-driven regime and with an energy-driven propagation once the contribution of the molecular outflow is included. On the contrary, the ionised outflow in SDSS J1353+1138 appears to be unrelated to the nuclear scale energetics, likely requiring either the presence of a massive molecular outflow or a high variability among the QSO activity.

3. By comparing our inferred results with those of the small sample of known QSOs from the literature, each hosting both sub-pc scale UFOs and neutral or ionised winds on galaxy scales, we found that the momentum- and energy-driven frameworks describe all the observed targets very well, with the exception of SDSS J1353+1138 and IRAS 17020+4544. Therefore, these driving mechanisms appear to explain how the energy released by the AGN activity is coupled with the galactic ISM, thus driving the wind propagation on a large scale.

Altogether, the observations presented in this work provide important pieces of information on the long sought-after ‘engine’ of large scale outflows and feedback, for the first time, at the crucial epoch for AGN feedback in galaxies, highlighting once again the power of integral field spectroscopy in this type of study. Follow-up CO observations of these sources will be necessary to confirm the molecular outflow detection and to constrain its energetics in HS 0810+2554 and to test whether a massive molecular outflow is responsible for the mismatch between the wind energetics at different scales observed in SDSS J1353+1138.

Acknowledgments

We thank the anonymous referee for comments and suggestions, which have improved the paper. We acknowledge support from the Italian Ministry for University and Research (MUR) for the BLACKOUT project funded through grant PRIN 2017PH3WAT.003. MP is supported by the Programa Atracción de Talento de la Comunidad de Madrid via grant 2018-T2/TIC-11715. MP acknowledges support from the Spanish Ministerio de Economía y Competitividad through the grant ESP2017-83197-P, and PID2019-106280GB-I00. Some data shown in this work were obtained from the Mikulski Archive for Space Telescopes (MAST).

References

Appendix A: Approximated reconstruction of the unlensed outflow in HS 0810+2554

As mentioned in Sect. 4.2, for HS 0810+2554 we were able to obtain a partial 2D-reconstruction of the outflow emission in the source plane by using the gravitational lens fitting-code presented in Rizzo et al. (2018), which adopts the lensing operator described in Vegetti & Koopmans (2009). This depends, in particular, on lensing operators describing the lens as a power law plus a shear component (see Sect. 2.2 in Rizzo et al. 2018 for details).

We began by simulating possible intrinsic geometries of the outflow. We considered multiple conical configurations, differing in radius and aperture angle, while the position angle of the cone was fixed to θ ∼ 130°, namely, the direction in which we observed the major [O III] outflow emission in the lens plane (see Sect. 3.3). In order to reduce further the number of free parameters characterising the outflow geometry, we fixed the origin of all simulated cones to the position of the emission centroid of the unlensed outflow, obtained through a first tentative full 2D-reconstruction with the lens fitting-code by Rizzo et al. (2018) and a lens model fixed to that found in Nierenberg et al. (2020). In fact, even though it could not be used to establish the intrinsic extent of the [O III] outflow because of the presence of residual effects from the PSF-deconvolution, it still provided reliable constraints of the centre position of the background outflow emission. Each cone was modelled as a uniform distribution of homogeneous point-like sources, and then each point-source was individually lensed forward (i.e., singularly mapped and magnified in the lens plane) through the reconstruction-algorithm of Rizzo et al. (2018), keeping the lens model fixed to Nierenberg et al. (2020). As a result, we obtained the whole forward lensed image of each starting background cone, differing in radius and aperture angle.

In order to establish the intrinsic size of the outflow, we selected those cones with a radius compatible with the detected [O III] outflow emission by visually comparing the extent of the forward lensed emission produced by a given radius background cone, with the maximum distance at which we observe the [O III] outflow, once corrected for the SINFONI-PSF, that is, R_max = 9.5 pixels⁷. In Fig. A.1, we show an example of our forward lensing method: starting, for instance, from a background homogeneous cone with radius of 8 pixels and aperture of 60° (left) and its forward lensed image (right). The solid red circumference has a radius equal to R_max = 9.5 pixels, thus indicating the maximal extent reached by the observed (PSF-corrected) [O III] outflow emission (see caption of Fig. A.1 for a detailed description). In this way, we found that the cones with a radius ranging from 6.5 to 9.5 pixels were compatible with R_max. Taking the average of these more plausible radii and the maximum deviation from the mean as error, and converting into kpc-units, we estimated the intrinsic radius of the outflow to be R_out = (8.7 ± 1.7) kpc, with z = 1.508 ± 0.002 as the redshift we measured from the nuclear spectrum extracted during the BLR-fitting (described in Sect. 3.1.1).

Fig. A.1. Example of forward lensing (Rizzo et al. 2018) of a background homogeneous cone with radius of 8 pixels and aperture of 60° (left) into the respective lensed image (right). In both images, the orange stars, red and black ‘+’ indicate, respectively, the position of the cone origin, of the emission peak observed by SINFONI, and of the lens centre. The dashed and solid red circumferences in the right panel have a radius equal to the intrinsic size of the cone in the source plane (i.e., 8 pixel in the example shown here) and to the maximum extent reached by the [O III] emission observed by SINFONI, i.e., Rmax = 9.5 pixels; both circumferences are centred in the observed emission peak (red ‘+’). The dashed and dotted black lines identify the position angle (∼130°) and the aperture (here 60°) of the intrinsic cone. We note that the conical source does not intercept the centre of the lens, thus subject to a low magnification.

For the estimate of the magnification factor with its error, we considered multiple aperture angles: 30°, 45°, 60°, and 90°, for each defined-radius cone, and calculated the magnification factor as the ratio between the total flux in the lens plane and the total flux in the source plane. All simulated conical configurations provided low magnification factors that are weakly dependent on the assumed geometry of the cone (they differ by a factor ∼1.4 at most). This follows from the fact that no simulated cones intercept the lens caustics and extend to large scales, where the lens magnification is reasonably lower. Therefore, to determine the outflow total mean magnification, we focused only on the magnification values obtained in correspondence to the established range of more plausible outflow radii, that is, 6.5−9.5 pixels, and averaged over them, finding μ_out = 2.0 ± 0.2, where the maximum deviation from the mean has been taken as error. We used this value to correct the observed [O III] outflow flux (determined in Sect. 3.3), thus obtaining F_out = (1.9 ± 0.2) × 10⁻¹⁵(2.0/μ_out) erg s⁻¹ cm⁻².

All Tables

All Figures

Fig. 1. Lensed images of HS 0810+1154 (left) and SDSS J1353+1138 (right). Left: HST ACS F555W image of HS 0810+2554 showing the four magnified images of the background quasar in fold lens configuration: the C and D images are spatially resolved, while the pair A+B is blended together. At the centre, we can see the emission from the foreground lens galaxy. Right: V and H-band images of SDSS J1353+1138 taken at the UH88 telescope (upper panels) and corresponding images after the subtraction of A and B components (lower panels), clearly showing the lens galaxy (component G). Image from Inada et al. (2006).

In the text

Fig. 2. Representative scheme of our fitting procedure. Top panels: SINFONI J-band spectra of HS 0810+2554 (left) and SDSS J1353+1138 (right), zoomed in the spectral region of [O III] and Balmer hydrogen emission lines. Both spectra were extracted using an aperture of $\sim 0\stackrel{″}{.}44$-radius, centred on the peak of the AGN continuum emission (located on image A in SDSS J1353+1138). Black dashed lines show the data, while red solid lines show the best-fit models obtained in phase II of our analysis. The latter is partitioned into the contributions of various components: AGN continuum (yellow), BLR emission in Balmer hydrogen emission lines (green) and Fe II (purple), and narrow line emission from the diffuse gas (light blue) represented as sum of multiple Gaussian components. Middle panels: J-band spectra extracted from the subtracted cubes, with the same aperture used in the top panels. Subtracted data (black lines) are compared to the best-fit models (red lines) resulting from the finer emission lines (EL) modelling implemented in phase III. The green lines highlight the outflow component alone. Bottom panels: residuals obtained by subtracting the full EL model from the subtracted spectra.

In the text

Fig. 3. Best-fit models of nuclear spectra of SDSS J1353+1138 extracted from an aperture of $0\stackrel{″}{.}3$ radius, centred on image A (upper panel) and on image B (lower panel). The various spectral components, the total model and data are represented with different colours (see the plot legend). In spectrum A, a broken power law distribution is perfectly suited to reproduce the BLR–Hβ profile, while in spectrum B, an additional broad Gaussian component was required to adequately reproduce the broad peak in the BLR–Hβ line profile, which is entirely dominant over the barely detected Hβ narrow component (solid light-blue line). The two spectra clearly differ from each other mostly for the lack of [O III] detection and the presence of a prominent blue wing in the Hβ line profile in spectrum B.

In the text

Fig. 4. Spatially resolved ionised outflows in HS 0810+2554 (top) and in SDSS J1353+1138 (bottom). Left panels: maps of the ratio between the [O III] outflow flux and the BLR flux (both from best-fit models). Coloured pixels refer to a S/N ≳ 2 (S/N ≳ 3) on the full [O III] emission line (i.e., narrow + outflow components) for HS 0810+2554 (SDSS J1353+1138). The positions of the continuum emission peak of the lensed images are marked with white ‘+’, while the dotted white lines indicate the contour levels of the BLR emission at 75%, 50%, and 25% of its peak; for SDSS J1353+1138, it is represented also the level at 90%. Right panels: normalised intensity profiles along the pseudo-slit (black dotted-dashed lines in the ratio map) and in circular annuli of increasing radius for HS 0810+2554 (top) and SDSS J1353+1138 (bottom), respectively: BLR model (red lines), [O III] outflow model (dotted blue lines) and [O III] blue wing from data (cyan lines). Blue points represent ratio values of the [O III] outflow flux over the BLR flux and they are referred to the right-hand logarithmic scale. The dashed black line in the plot of SDSS J1353+1138 corresponds to the radial distance of the centre of image B.

In the text

Fig. 5. Moment-0 (intensity field), moment-1 (velocity field), and moment-2 (dispersion field) maps of the [O III] line emission in HS 0810+2554 (left) and in SDSS J1353+1138 (right). The maps for the total, narrow and outflow components are shown separately, reporting only spatial pixels with a S/N equal or higher than 2 for HS 0810+2554, and than 3 for SDSS J1353+1138. The black ‘+’ indicates the emission centroid in each QSO, while the dotted lines represent the contour levels of the BLR emission at 75%, 50%, and 25% of its peak.

In the text

	Fig. 6. v10, v90, and W80 maps of the [O III] emission line in HS 0810+2554 (upper panels) and in SDSS J1353+1138 (lower panels). We applied the same cut in S/N as in the moment maps of Fig. 5, that is, with a S/N equal or higher than 2 and for HS 0810+2554 and higher than 3 for SDSS J1353+1138.
In the text

Fig. 7. Momentum-boost versus wind velocity diagram for both UFOs (represented as circles) and galactic outflow components (represented as squares) of HS 0810+2554 and SDSS J1353+1138. The dashed lines represent the predictions for a momentum-conserving regime (horizontal) and energy-conserving mode (diagonal). For HS 0810+2554, also the point relative to the ionised+molecular large-scale outflow is shown, identified by the star symbol.

In the text

Fig. 8. Ratio between the galaxy-scale and sub-pc scale outflow momentum rates for different QSOs hosting UFOs. Measurements for individual objects are shown in blue with the respective errorbars, using different markers according to the gas phase of the observed large scale outflow. The galaxy points are ordered by increasing LBol. The horizontal dashed line shows the prediction for a momentum-driven wind (ṗout/ṗUFO = 1), while the orange rectangles indicate the individual predictions for energy-driven winds. Filled and empty squares represent ionised outflow measurements based on Hα and [O III] emission, respectively. For HS 0810+2554, our [O III]-based measurement is shown both alone and combined with the CO-measurement of the molecular outflow (Chartas et al. 2020), with the respective symbol and combination of symbols (see the plot legend).

In the text

Fig. A.1. Example of forward lensing (Rizzo et al. 2018) of a background homogeneous cone with radius of 8 pixels and aperture of 60° (left) into the respective lensed image (right). In both images, the orange stars, red and black ‘+’ indicate, respectively, the position of the cone origin, of the emission peak observed by SINFONI, and of the lens centre. The dashed and solid red circumferences in the right panel have a radius equal to the intrinsic size of the cone in the source plane (i.e., 8 pixel in the example shown here) and to the maximum extent reached by the [O III] emission observed by SINFONI, i.e., Rmax = 9.5 pixels; both circumferences are centred in the observed emission peak (red ‘+’). The dashed and dotted black lines identify the position angle (∼130°) and the aperture (here 60°) of the intrinsic cone. We note that the conical source does not intercept the centre of the lens, thus subject to a low magnification.

In the text