© The Authors 2024

1. Introduction

Since the discovery of supermassive black holes (SMBHs, 10⁶–10¹⁰ M_⊙) in the centres of galaxies (Lynden-Bell 1969; Rees 1984), it has become clear that these SMBHs are closely connected to the formation and evolution of galaxies. There is a correlation between the properties of the host galaxy and the properties of the black hole for example the black hole mass and stellar velocity dispersion (Tremaine et al. 2002; Kormendy & Ho 2013; McConnell & Ma 2013; DeGraf et al. 2015). It is believed that a mechanism connects the innermost regions, where the black hole’s gravitational field is dominant to the larger scales. One proposed mechanism is the presence of outflows in the form of energetic jets or winds (Alexander & Hickox 2012; Harrison 2017; Harrison & Ramos Almeida 2024). These outflows are believed to play a crucial role in galaxy evolution by regulating the accretion onto and ejection of material from the SMBH (King & Pounds 2015). During this period called the ‘feedback phase’, the radiation pressure accelerates in the form of winds from the accretion disk in rapidly accreting sources into the interstellar medium (ISM), driving outflows and providing an efficient feedback mechanism (Silk & Rees 1998; Hopkins & Elvis 2009).

The outflows in active galactic nucleus (AGN) host galaxies have a complex nature, consisting of multiple phases (ionised, molecular and neutral; Cicone et al. 2018; Harrison & Ramos Almeida 2024). The ionised phase has been observed to extend over kilo-parsec scales in both local (Kakkad et al. 2022; Speranza et al. 2022; Toba et al. 2022; Shen et al. 2023; Venturi et al. 2023; Zanchettin et al. 2023) and high redshift (z) AGN (Carniani et al. 2015; Scholtz et al. 2020; Cresci et al. 2023) using spatially resolved data. While the presence of ionised outflows can be detected by broad or shifted wings in the [O III]λ5007 Å emission line in integrated spectra (Mullaney et al. 2013; Zakamska et al. 2016; Toba et al. 2017), high-resolution and high signal-to-noise ratio observations with integral field units (such as MUSE, KMOS and SINFONI) have played a crucial role in determining their extent (radius) and kinematics (velocity and geometry) (e.g. Scholtz et al. 2020; Kakkad et al. 2022, 2023; Cresci et al. 2023). These parameters are essential for calculating the mass outflow rates and energetics of the outflows.

Theorists suggest that a short, luminous and dust-enshrouded phase plays a significant role in the accretion of SMBHs when their activity is at its peak (high Eddington ratio; Hopkins et al. 2005; Di Matteo et al. 2005; Blecha et al. 2018). The reddest quasars are known to have these characteristic high Eddington ratios (Urrutia et al. 2012). The dusty environments produce moderate obscuration and reddening that makes it difficult to identify these objects. The most effective way to identify these sources and gather information about nuclear obscuration and ongoing AGN activity is through mid-infrared (mid-IR) and large-area hard X-ray surveys, coupled with multi-wavelength data such as optical surveys. Previous studies have utilised this approach to select these rare sources to provide insights into AGN activity (e.g. Perna et al. 2015; Cresci et al. 2015; Brusa et al. 2016; Musiimenta et al. 2023) and others to explore their radio properties (e.g. Klindt et al. 2019; Fawcett et al. 2020).

To obtain significant samples of obscured quasars, it is essential to have an all-sky X-ray survey. This need is being addressed by the eROSITA (extended ROentgen Survey with an Imaging Telescope Array) telescope (Predehl et al. 2021), thanks to its large field of view, spectral resolution, and angular resolution. The first data release for the all-sky survey is now out (eRASS1, Merloni et al. 2024).

To test the capabilities of the eROSITA survey, a mini-survey called eFEDS (eROSITA Final Equatorial-Depth Survey; Brunner et al. 2022; Salvato et al. 2022) was conducted during the performance and verification phase, covering ∼140 deg² of the sky in the Subaru Hyper-Supreme Cam SSP field and reaching the average depth of the planned all-sky survey. Musiimenta et al. (2023) exploited the large eFEDS catalogue (∼11 000 AGN at z = 0.5–3) by applying, for the first time simultaneously, a combination of 5 selection methods to isolate ∼1400 candidate AGN in the feedback phase. The five methods were based on optical and mid-IR colours and optical and X-ray spectral properties and they include;

• i − W3 > 4.6

• 

• 

• (r−W1> 4) ∧ (i − W4 > 7)

• (N_H > 10^21.5 cm⁻²) ∧ (λ_Edd> effective Eddington limit for different values of N_H).

Out of the selected candidates, only one source, eFEDS J091157.4+014327 (referred to as XID 439 in Brusa et al. 2022 and ID 608 in Salvato et al. 2022; Musiimenta et al. 2023 and hereafter) has been selected based on all five different criteria. It is therefore, referred to as a red quasar based on the combination of optical and mid-IR colours (Musiimenta et al. 2023) see also Urrutia et al. (2012), Brusa et al. (2022). The SDSS spectrum available for ID 608 (z = 0.603) displayed a broad redshifted [O III]λ5007Å line ([O III] hereafter) with a Full Width at Half Maximum (FWHM) of approximately 1650 km s⁻¹. The observed luminosity of the [O III] emission (L_O III) was measured to be around 2.6 × 10⁴² erg s⁻¹, which is one to two orders of magnitude larger than what has been observed in local Seyfert galaxies. However, it is comparable to the [O III] luminosities of 1.26 × 10⁴³ erg s⁻¹ (Reyes et al. 2008) and 2.4 × 10⁴² erg s⁻¹ (http://research.iac.es/galeria/cristina.ramos.almeida/qsofeed/) for a sample of optically selected QSO2s at 0.5 < z < 0.83 and a subsample of type 2 quasars at z ∼ 0.14, respectively

ID 608 has also been previously reported in Brusa et al. (2022) as an archetypal quasar in the feedback phase selected by eROSITA from the much smaller hard X-ray selected catalogue (Nandra et al. 2024). From the X-ray point of view, ID 608 exhibits type 2 characteristics, with a luminosity of approximately 6.5 × 10⁴⁴ erg s⁻¹ and a neutral hydrogen column density (N_H) ≈ 2.75 × 10²² cm⁻², when a simple power-law model is assumed. Waddell et al. (2024), based on a comprehensive spectral analysis of the eFEDS hard X-ray sample, reported that ID 608 is best fit with an ionised warm absorber, supporting a scenario in which the outflowing ionised gas seen in the optical band may imprint its presence also in the X-ray spectrum. The Eddington ratio of ID 608 (L_bol/L_Edd) was derived in Brusa et al. (2022) as ∼0.25 and revised later in Musiimenta et al. (2023) as 0.99 using the black hole mass (log(M_BH) = 8.4) from Shen et al. (2023). These properties and the fact that it is the only candidate quasar in the feedback phase selected simultaneously from multiple diagnostics make ID 608 an ideal target for further studies on the spatially resolved analysis of ionised winds in obscured quasars at intermediate redshifts and their interaction with the surrounding environment.

This study aims to analyse the observational data for ID 608 from the Multi Unit Spectroscopic Explorer (MUSE) instrument. Our focus is to investigate the kinematics of the outflows observed in this particular system and the environment around it and shed light on the underlying mechanisms responsible for these outflows. This paper is organised as follows. Section 2 describes the MUSE observations of our source and data reduction. Section 3 describes methods and results concerning companion galaxies, kinematics, and ionised outflow properties. Section 4 describes the interacting system discovered and mass outflow rates and energetics. Conclusions are presented in Sect. 5. Throughout the paper, we adopt the cosmological parameters H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3 and Ω_Λ = 0.7(Spergel et al. 2003).

2. MUSE observations and data reduction

The MUSE instrument located at the Very Large Telescope (VLT) of the European Southern Observatory (ESO) combines the power of a wide-field imager and an integral field spectrograph allowing scholars to simultaneously obtain detailed images and spectra of galaxies. It operates in various optical wavelengths, capturing light across a wide spectrum from ultraviolet to near-infrared. Its high-resolution capabilities allow for detailed mapping of structures within galaxies revealing intricate details of phenomena such as AGN outflows and star-forming regions among others.

We conducted observations of ID 608 at RA(J2000); 09:11:57.557 and Dec(J2000); +01:43:27.54 using the MUSE instrument (Bacon et al. 2010) as part of proposal ID:108.22M6.001. Our observations spanned five nights; December 28th 2021, December 29th 2021, January 3rd 2022, January 4th 2022, and January 5th 2022 in Wide Field Mode (WFM). These yielded an average seeing of 0.6″ for a total observation time of 6.1 hours. The field of view is 1 × 1 arcmin² with a pixel size of 0.2 arcsec. Our observations cover a spectral range of 4650–9300 Å sampled at intervals of 1.25 Å with a spectral resolution R ≈ 2800.

Fig. 1. Full data cube over the whole wavelength range covered in our MUSE observation and the continuum and [O III] emission for size-reduced cubes. Panel a (left) is the image of the full MUSE’s FOV (1′×1′≅402 kpc × 402 kpc) obtained after data reduction. ID 608 is marked by a purple circle. Panel b (centre) shows an inner zoom of the cube (30″ × 30″ ≅ 201 kpc × 201 kpc) with the continuum in blue, extracted from two regions of the spectrum (around [7825, 7900] Å and [8100, 8150] Å observed wavelengths) and the [O III] emission in red extracted from observed wavelengths 7890 Å to 8100 Å. Panel c (right) shows a further zoom (8″ × 7.5″ ≅ 53 kpc × 50 kpc) of the [O III] emission centred on the quasar position. The position of the nucleus of the quasar obtained from the maximum of the [OIII] emission is marked by the red star. The contours in both panels b and c represent 3σ, 5σ, 15σ, 30σ, and 60σ. North is up, East is left.

The ESO archive provides both raw and reduced data cubes that are ready for science analysis. However, to ensure the highest data quality, we opted to work with the raw data and perform the data reduction ourselves following the manual from the AIP website¹. This process uses EsoReflex (Freudling et al. 2013) and allows us to have full control over the reduction process instead of relying on solely automatically reduced MUSE data cubes generated by the ESO pipeline.

For the raw data pre-processing, we utilised the MUSE data reduction software (Version 2.8 of the pipeline Weilbacher et al. 2020). We used all the files provided by the ESO archive for our observation, except for the line spread function (LSF) profile and geometry and astrometry tables. To ensure accuracy we calculated the LSF profile ourselves and substituted the geometry and astrometry tables with more recent ones obtained in December 2021, which were closer in time to our observation compared to the ones provided by ESO from April 2021.

The pre-processing steps included bias subtraction, dark calibration, flat-fielding, wavelength calibration, LSF calibration, and twilight calibration. Since our source was observed for five nights, we combined the preprocessed data into a single datacube using muse_exp_align and muse_exp_combine recipe as explained in the manual. The full cube collapsed over the entire wavelength range is shown in panel a of Fig. 1, where the nucleus is marked with a red star.

Since we are interested in ionised outflows traced by the [O III] emission that only extends over a small portion of the FOV, we created a smaller cube around the position of ID 608. Figures 1b, c show the reduced datacubes of sizes of 30″ × 30″ (201 × 201 kpc) and 8″ × 7.5″ (54 kpc × 50 kpc), respectively. This [O III] emission line falls within the wavelength range of 7890–8100 Å at our source redshift z = 0.6031. Upon examining the intensity map in panel c of Fig. 1, we observe extended [O III] emission in the southwestern (SW) region of the source and, to a lesser extent, also in the northeastern (NE). The structure of this source’s [O III] emission as seen in panel c resembles the shape of the Goldfish. Therefore, we name this quasar a ‘Goldfish’ galaxy for future reference.

3. Methodology and results

3.1. Gas distribution and searching for companion galaxies

To understand how the [O III] emitting gas is distributed within our system, we binned the data into different velocity channels. The channels included velocities ranging from −650 to −250 km s⁻¹, −250 to −50 km s⁻¹ (the ‘blue’ channel), −50 to 50 km s⁻¹ (the ‘systemic’ channel), 50 to 250 km s⁻¹ (the ‘red’ channel), and 250 to 650 km s⁻¹.

We created flux maps or colour maps from these velocity channels, which are shown in Figs. 2a–f. The top panels (a–e) show all the velocity channels individually, while panel f shows a combined RGB image of the blue, systemic, and red channels. To determine the centroid of the quasar, we used a fixed wavelength (centroid) of 8026.5 Å (corresponding to the [OIII]5007 emission line at z = 0.6031). The centroid has been determined using QFitsView (Ott 2012) by subtracting the continuum from the quasar nuclear spectrum extracted within an aperture of 0.6″ × 0.6″ (4 × 4 kpc). The pixel with the maximum [O III] emission was considered the ‘quasar nucleus’ and is marked by a star in all the panels. Its wavelength with the maximum [O III] emission was used as the centroid of the quasar (quasar centroid hereafter).

Fig. 2. RGB colour image, velocity channels and velocity distribution for the [O III] emission obtained from data without fitting. The different panels are named as a, b, c, d, e, f, and g. Panels a–e: [O III] intensity maps extracted in different velocity channels, as labelled. The arrows indicate the position of companion one (C1). Panel f: RGB colour image constructed from three [O III] velocity bins ([−250, −50] km s−1, [−50, −50] km s−1 and [50, 250] km s−1). Panel g shows the velocity distribution map attained without fitting to show the different regions as will be referred to in this study. The regions are NE, NW, SW, MSW, and FSW. The star shows the position of the nucleus of the quasar. The figures in all the panels are at the same scale and extracted from 8″ × 7.5″ ≅ 53 kpc × 50 kpc FOV.

We generated a velocity map from the data without fitting, which involved extracting a spectrum for each pixel around the 7890–8100 Å wavelength range and determining the wavelength corresponding to the peak of the [O III] emission line. The velocity map (shown as ‘data velocity distribution’) in Fig. 2g, indicates the velocity distribution of the gas that is asymmetric with velocity gradients. It also shows the different regions as explained below.

In panel a, the gas distribution shows an extended structure towards the south-west (SW) direction and a small extension towards the north-east (NE) direction. This SW extension is more visible in panel b, where an additional extension in the NE appears as a second peak. This peak corresponds to higher velocities of up to ∼600 km s⁻¹ in the NE region of panel f. Likely, this second peak located at 1.4″ or 9.4 kpc away from the nucleus of the quasar represents a companion galaxy or a merger (C1 in the figure, see further discussion below). The colour maps of velocity channels [−650, −250] km s⁻¹ and [−250, −50] km s⁻¹ suggest the presence of a bubble-like structure in the middle SW position (MSW region in panel f and hereafter).

The systemic velocity, as shown in panel c, is mainly centralised, with only a small portion showing up in the SW position (SW region in panel f and hereafter). Additionally, we observed the emission in the red velocity channel, shown in panel d, with most of it centralised near the quasar nucleus and extra emission in the SW region and far SW position (FSW region in panel g and hereafter). In panel e, the velocity channel [250, 650] km s⁻¹ mainly exhibits centralised emission near the quasar nucleus. Complete velocity channels from −1200 km s⁻¹ to 1200 km s⁻¹ at intervals of 100 km s⁻¹ are shown in Appendix A with more visible extended structures.

To search for additional emitting sources, we analysed the spectra of the sources within the larger FOV closer to the quasar (201 × 201 kpc) and estimated their redshifts. Based on the presence of continuum emission and spectral lines (absorption or emission), in addition to C1 we discovered other two companion galaxies (C1 and C2) with redshifts consistent with that of the quasar.

The spectra of ID608 and the three companion galaxies are shown in Fig. 3. The redshifts of these sources were estimated from different emission and absorption lines observed, that is [Ne V] 3425, [O II] 3726, Ca K and H, Hδ 4101, Hβ 4861, [O III] 5007, and Mg 5175. [Ne V] 3425, Hβ and [O III] 5007 were only used if they appeared significant in the spectrum. All these lines are marked in Fig. 3. C1 is found at a redshift of z = 0.6010 ± 0.001, while ID 608 is at z = 0.6024 ± 0.0008. Both C1 and ID 608 exhibit strong emission lines in their spectra. C2 on the other hand, displayed weak emission lines and is continuum-dominated with a redshift of z = 0.6033 ± 0.0023. In the case of C3, absorption lines were observed and it is also continuum-dominated with a redshift of z = 0.6012 ± 0.0028 displaying a significant 4000 Å break.

Fig. 3. Spectra of the quasar and the companion galaxies C1, C2, and C3 with their estimated redshifts. The spectrum in red indicates the quasar (ID 608), companion galaxy one in black, companion galaxy two in yellow, and companion galaxy three in blue.

The three companions are marked as yellow stars in Fig. 4. C1 is at 1.4″ (9.4 projected kpc), C2 is at 1.3″ (8.7 projected kpc) and C3 is at 7 arcsec (47 projected kpc) away from the quasar nucleus. The extended structure in the SW direction that is visible in the [O III] emission stretches towards the position of C3. In the bottom left corner of Fig. 4, we show the zoom-in map of the continuum, [O III], and Hβ emission around the quasar with positions of C1 and C2. From the zoom-in map, the peak of C2 is visually significant in the continuum emission. No continuum emission is seen in the MSW, SW, and FSW regions. C1 is slightly observed in Hβ while C2 is not, this is also seen in the spectra of the quasar and the companion galaxies in Fig. 3.

Fig. 4. Colour map of the continuum emission overlapped with the [O III] emission in 30″ × 30″ ≅ 201 × 201 kpc FOV. The continuum emission is shown in blue with corresponding white contours and the [O III] in red (no contours). The companion galaxies are marked by yellow stars, while ID 608 nucleus is marked by a red star. In the bottom left of this figure, we show the zoom-in map of the continuum (blue contours), [O III] (red without contours), and Hβ (green contours) emission. The positions of C1 and C2 are indicated that do not appear to emit in Hβ. The contours in the zoom in represent 3σ, 5σ, 15σ, 30σ, and 60σ.

In summary, the analysis of the velocity distribution and the spectra suggests a complex system with various components. The presence of different companion galaxies at the same redshift as the quasar suggests the possibility of an ongoing merging of galaxies or galaxy interactions. Moreover, there is an extended structure (MSW, SW, and FSW) in the [O III] emission without continuum emission. This is likely to have resulted from an ongoing merger and/or interaction between C3 and the quasar.

3.2. Analysis of the [O III] emission from different regions

We extracted spectra from a 0.6″ × 0.6″ box from the regions (NE, NW, SW, MSW, FSW and the centre) where we observed [O III] emission with velocity gradients. We fit the spectra with multi-Gaussian components (n = 1–6) around [O III]λλ4949,5007 Å doublet using the lmfit Python package and the in-house script described in Speranza et al. (2024). In this method, more than one Gaussian is fitted when the χ² improves more than 10% with respect to the previous fit. It is worth noting that the parameters of both emission lines were tied together and the flux ratios between the [O III]λ4959 and [O III]λ5007 lines were fixed at a ratio of 1:3 (Osterbrock 1981). The velocity dispersion of each kinematic component is left to vary. The different spectra are shown in Fig. 5.

Fig. 5. Spectra relative from the six regions shown in Fig. 2. The spectra are extracted from a small area of a 0.6″ × 0.6″ box at the given position in the FOV similar to Fig. 2f. Within the six regions shown, only the nuclear region and NE are best fitted by more than two Gaussian functions. One or two Gaussian components can adequately describe all the remaining spectra. The legend for all the spectra is represented in the spectrum from the NE region where the black solid spectrum shows the data, the red dotted shows the best fit, the grey dotted spectrum shows the residual and all the other colours (blue, magenta, solid red, green) show the different components fitted. The three vertical lines highlight the centres of the Hβλ4861, [O III]λ4949 and [O III]λ5007 emission lines. The yellow stars indicate the positions of C1 and C2. In the bottom right of the map, the white arrow points towards the position of C3 that is outside this map size. The different regions that is NE, NW, MSW, SW, FSW and ID 608 are indicated in the respective spectrum.

The spectrum from the NE region where an extra emission peak (C1) was observed shows an [O III] emission peak at a velocity offset of ∼ − 235 km s⁻¹ from the quasar centroid. We measured the line ratios to test the possibilities of AGN presence or quasar ionisation in the different regions and/or companion galaxies, if necessary. In the quasar, the line ratio log([O III]/H_β) ∼ 1. Since C1 is detected in both [O III] and H_β emission, we measured the line ratio log([O III]/H_β) ∼ 0.9 which is compatible with the values that indicate the possible AGN ionisation. To check whether ID 608 is responsible for powering C1 or if there is a secondary AGN, we followed the method in Perna et al. (2023). By assessing whether the line emissions in C1 can be attributed to the ionisation originating from ID 608, we calculated the ionisation luminosity using equation one (Eq. (1)) in Perna et al. (2023). We obtained ionizing luminosity of 3 × 10⁴³ erg s⁻¹, which is lower than the incident ionizing luminosity (see Perna et al. 2023 for the formula) from ID 608 of 2 × 10⁴⁵ erg s⁻¹. These values indicate that ID 608 serves as the ionizing source for C1 and exclude the possibility of the presence of a secondary AGN in C1.

The spectrum from the NW region (C2) shows a broad [O III] emission which could be due to the contamination from the quasar, though the continuum map shows significantly the presence of another source. Its [O III] emission peak observed is at a velocity offset of ∼ − 93 km s⁻¹ from the peak of the quasar. C2 is detected in [O III] but not in H_β.

We do not derive line ratios in C3 since it is not detected in [O III] and H_β. Its spectrum discussed in Sect. 3.1 and Fig. 3 is similar to that of an elliptical galaxy.

The spectrum from the MSW region shows a blue shift in relation to the quasar centroid, while the spectrum from the SW region perfectly aligns with the quasar centroid. The spectrum from the FSW region exhibits a distinct red shift when compared to the quasar centroid. As seen in Fig. 4, the H_β detection in these regions is not significant because the Balmer line is faint (even in the centre). The [O III] extended emission in this extended structure is likely due to the interactions between the quasar and C3. The following section will explore the debate regarding whether these extended structures (MSW, SW and FSW) are indicative of outflows.

3.3. Spaxel-by-spaxel analysis of the [O III] emission

3.3.1. Is [O III] emission resolved?

Before performing a spaxel-by-spaxel [O III] fitting to produce kinematic maps, we followed the method in Speranza et al. (2024) to determine whether the [O III] emission is spatially resolved or not, which is based on Carniani et al. (2013) and also used in previous studies (e.g. Carniani et al. 2015; Kakkad et al. 2020, 2023). This method known as ‘point spread function (PSF) subtraction’ allows us to assess the contribution of the PSF to the emission from the AGN. If the [O III] emission is resolved, the spectrum obtained from pixels away from the central region should differ from the spectrum obtained from the centre. For a complete description of this method, we refer the reader to Sect. 4.2 in Speranza et al. (2024). To perform PSF subtraction, we extracted a spectrum from five central pixels and fitted it with a multi-Gaussian model (see Sect. 3.2) to obtain a ‘PSF model’. Then, we extracted the [O III] spectrum for each pixel in the data cube and subtracted the PSF model. Before subtraction, a normalisation factor (considering the entire wavelength in the residual spectrum in Fig. 6) is applied to scale the amplitude of the PSF model. This scaling factor minimises the residues between the PSF model and the spectra from other pixels, as in Speranza et al. (2024; see also Kakkad et al. 2023). This process yielded a set of residual spectra for each pixel, as shown in the left panel of Fig. 6. To visualise the residuals across the cube, we collapsed the residual spectra and created a residual map, displayed in the right panel of Fig. 6. We observed excess emission at > 3σ near the central region in both the residual map and spectrum. This indicates that the [O III] emission is resolved.

Fig. 6. Residual spectrum and map after PSF subtraction show that the [O III] emission in our quasar is resolved. The left panel shows the residual spectrum analysis of a randomly selected pixel. The two vertical lines serve as markers highlighting the wavelengths of [O III]λ5007 and [O III]λ4949 emission lines for the centre of ID 608. The right panel showcases a residual map resulting from collapsing the residual spectrum of each pixel. The red square on the right panel corresponds to the specific pixel from which the residue spectrum was extracted in the left panel and the centre is marked by a white star.

To properly measure the outflow energetics, it is important to account for the smearing effect caused by the spatial PSF, which can result in an overestimation of outflow measured values, particularly in objects where the nucleus outshines the extended emission (Husemann et al. 2016). This effect can easily be addressed in unobscured quasars (e.g using QDeblend3D presented in Husemann et al. 2013, where the nuclear spectrum is scaled using the PSF modelled from broad line region (BLR) emission lines like Hβ and subtracted from the spectrum of each off-nuclear spaxel).

We note that ID 608 is a type 2 (obscured) AGN. Modelling the emission from the (faint) BLR is therefore complex. In Brusa et al. (2022) and Musiimenta et al. (2023), using the integrated SDSS spectrum, the best-fit model for the Hβ complex included a narrow component, and an additional broad component with kinematics consistent with the ones observed in [O III] tracing the outflows, without any contribution from the unresolved BLR. The Hβ extracted from the central 0.6″ × 0.6″ box in the MUSE data cube is very faint (see the bottom left spectrum in Fig. 5) and we do not detect any BLR component. The PSF is therefore likely to go to zero in the spaxels very close to the centre, since it is constructed by scaling the nuclear spectrum to match the BLR wings in each spaxel. Therefore, no PSF deblending is required as BLR wings are undetected at radius, r > 4 pixels.

In addition, we note that although the method removes emissions that can be ascribed to unresolved BLR, it may also wrongly remove flux from the components we are interested in, causing a severe underestimation of the outflow parameters. This is particularly important because the physical scale of our PSF is ∼4 kpc, large enough to contain most of the unsettled and outflowing gas.

Given that we cannot detect clear emission from BLR and at the same time we do not want to disregard a significant portion of the outflowing flux just because it is unresolved, we will use the data cube without PSF subtraction for the subsequent analysis. At the end of Sect. 3.4, we comment on the consequences of this assumption.

3.3.2. Kinematics

After confirming that the [O III] emission is spatially resolved, we fit the [O III]λλ4949, 5007Å doublet for each spaxel using multi-Gaussian components as described in Sect. 3.2. We used the non-parametric approach where no physical meaning is ascribed to the single Gaussian components to avoid spaxel-by-spaxel parametric fitting that is too complex because of the many kinematic components that would be involved. These components are aimed at reproducing the line profile. By analysing the best-fit model obtained from this approach, we determined the kinematics of the gas around the quasar.

The kinematic maps are represented in Fig. 7. We define the five percentile velocity (V₀₅) to trace outflow velocities of the approaching gas and 95 percentile velocity (V₉₅) to trace outflow velocities of the receding gas. The line width that contains 80% of the line emission flux that is, W₈₀ = 90 percentile velocity −10 percentile velocity, is an approximation of the velocity dispersion. This velocity dispersion is equal to FWHM/2.35 for a Gaussian and tells us how turbulent the gas is. These parameters were obtained for each pixel resulting in the generation of velocity and flux maps in Fig. 7.

Fig. 7. Results from the spaxel-by-spaxel multi-Gaussian fitting of 8″ × 7.5″ ≅ 53 kpc × 50 kpc region around the quasar. Top left panel: the flux map. Top right panel: the velocity peak. In the middle panel, we present the V05 and V95 maps indicating the velocities at the 5th and 95th percentile, respectively. Bottom left panel: the line width (W80). Bottom right panel: the number of Gaussians fitted in each pixel.

The flux map in the top left panel shows an extended emission in the NE and SW directions. The positions of C1 and C2 are indicated. C3 is outside this FOV, we indicate the direction of its position from the centre with the arrow.

The velocity peak (V_p) map in the top right panel displays systemic velocities close to the central region, followed by a sharp increase reaching ∼ − 600 km s⁻¹ in the NE and MSW regions. In FSW, there is a rise in velocities to ∼200 km s⁻¹. This velocity peak map helps us visualise the kinematics of the ionised gas in our data: the pattern is asymmetric with a velocity gradient but no significant rotation is seen. The velocities in the extended structure in the SW direction are low and can be attributed to the ordinary movement of the gas than outflows.

It is evident from the middle left panel that the quasar exhibits ionised outflows moving towards the observer, indicated by velocities (V₀₅) reaching up to a maximum of −1200 km s⁻¹. Within the spatial extension of up to ∼2″ (13.4 kpc) away from the centre, these outflows have blue-shifted velocities greater than −500 km s⁻¹. On the other hand, the middle right panel shows receding gas (V₉₅) with velocities up to 1000 km s⁻¹, which is confined to a smaller region near the AGN, about 6.7 kpc in size for velocities greater than 500 km s⁻¹. These velocities greater than 500 km s⁻¹ are high enough to be associated with outflows. The red-shifted velocities observed in these kinematic maps were also identified from the SDSS integrated spectrum in Brusa et al. (2022), Musiimenta et al. (2023) with a velocity shift of 100 km s⁻¹ and 171 km s⁻¹.

In the bottom left panel, the line width is notably high ranging from 600–1800 km s⁻¹ in regions where there is a high receding and approaching velocities. The high turbulence in these regions is associated with ionised outflows. As one moves farther away from the centre of the AGN, the velocity dispersion decreases. The velocity dispersion values in the MSW region are significantly lower than 500 km s⁻¹. Lastly, the bottom right panel displays the number of Gaussian fits performed in each pixel, providing further detail about the distribution and complexity of the observed data.

We tested a three-Gaussian fitting for n = 1–3 with and without Voronoi binning. We found similar results (see Appendix B).

3.4. Ionised outflow properties

In this section, we measure the ionised outflow properties including the mass of the outflow, the mass outflow rate and the kinetic power. To measure the outflow properties in this study, we construct flux maps for two different regions from the multi-Gaussian best-fit model described in Sect. 3.3 that is, the region with blueshifted outflow velocities [−2200, −250] km s⁻¹ (left panel of Fig. 8, blue region hereafter) and the region with redshifted outflow velocities [250, 2000] km s⁻¹ (right panel of Fig. 8, red region hereafter). This allows us to isolate the flux that is only related to the outflow. The outflow properties are then measured for two separate regions and final values are later obtained by summing.

Fig. 8. Flux maps for the blue-shifted and red-shifted velocities from the multi-Gaussian fitting. We select the region [−2200, −250] km s−1 to extract the flux map for negative velocities (the left panel). In the right panel, the flux map for the positive velocities [250, 2000] km s−1. The light blue, yellow, and green contours indicate the 4, 3, and 2σ levels, respectively. These contours help to determine the significance level of the measured flux values. The blue dashed line represents the maximum distance to the 3σ fitted ellipse in white.

To measure the mass of the outflow, we employed the approach presented in Carniani et al. (2015) using the formula,

(1)

where is the outflow luminosity of the [O III] emission line in units of 10⁴⁴ erg s⁻¹, n_e is the outflowing gas electron density in units of 10³ cm⁻³, 10^[O/H] is the oxygen abundance in solar units and C is a factor that encodes the ‘condensation factor’ of the gas clouds (see Cano-Díaz et al. 2012, for the details). The ionising gas clouds are always assumed to have the same density thus the condensation factor C is approximated to one and the metallicity of the outflowing material is always assumed solar. To determine the luminosity of the outflow (), we follow the same approach outlined in Speranza et al. (2024). We fit an ellipse at 3σ of the flux (white ellipse in Fig. 8) and measure the total flux within the ellipse as the flux associated with the outflow, for both the red and blue regions. The flux is corrected for dust extinction using the E(B − V) = 0.63 value reported in Kim et al. (2015) and relations presented in Cardelli et al. (1989) and Calzetti et al. (2000). The extinction-corrected [O III] luminosity is 1.94 × 10⁴² erg s⁻¹ and 1.47 × 10⁴² erg s⁻¹ for red and blue regions, respectively, for a total of ∼3.4 × 10⁴² erg s⁻¹. This value is comparable with the values reported in Brusa et al. (2022) (2.6 × 10⁴² erg s⁻¹) and Musiimenta et al. (2023) (2.2 × 10⁴² erg s⁻¹). The electron density can be measured directly using various methods such as the [SII]λλ6716,30 line ratio, the auroral and transauroral line method and the ionisation parameter method (see Davies et al. 2020 for detailed discussion of the methods). Unfortunately, the necessary emission lines for estimating the electron density are not available in our data. As a result, we assume an electron density (n_e) of 500 cm⁻³, consistent with the assumption made in Brusa et al. (2022) also within the range of measured electron densities at z ∼ 0.5 − 1.5 (e.g. Perna et al. 2015, 2017; Cresci et al. 2023). We calculate an ionised outflow mass of 15.5 × 10⁶ M_⊙ and 11.7 × 10⁶ M_⊙ for red and blue regions, respectively. The total outflow mass is 27 × 10⁶ M_⊙.

The spherical mass outflow rate can be deduced from Eq. (1) as

(2)

where V_out is the outflow velocity, R_out (in kpc) is the radius at which the outflow is computed and is the mean density of an outflow from the fluid field continuity equation with Ωπ as the solid angle covered by the outflow. Following Speranza et al. (2024), to measure the radius, we consider the maximum radius from the centre of the quasar to the 3σ ellipse in Fig. 8. This procedure is carried out separately for the red and blue regions. To account for the seeing effect, we subtract half of the seeing FWHM (half width at half maximum) in quadrature from the measured radius following Kang & Woo (2018). It is important to note that this correction only caters for the extension and does not account for the full smearing effect (see end of Sect. 3.3.1). We find a radius of kpc and kpc for red and blue regions, respectively. This value is comparable to 10 kpc in Brusa et al. (2022) and Musiimenta et al. (2023) where they assumed outflow extension as a half-light radius corresponding to the entire galaxy scale.

Regarding the measurement of the velocity, many methods have been used to measure the velocity for spatially resolved data (see Hervella Seoane et al. 2023, for different velocity definitions and comparisons). We adopt as the velocity of the outflow (V_out) the mean value between V₀₅ and V₉₅. We find the mean value, V_out = 807.2 km s⁻¹ (red region) and 772.6 km s⁻¹ (blue region). The values of V₀₅ and V₉₅ used to calculate V_out are shown in Table 1. We measure mass outflow rates of ∼ 5.9 M_⊙ yr⁻¹ (red region) and ∼ 4.0 M_⊙ yr⁻¹ (blue region). The total ionised outflow mass rate is 10 M_⊙ yr⁻¹.

The corresponding kinetic power of the ionised outflow is estimated using

(3)

This gives Ė_kin ∼ 1 × 10⁴² erg s⁻¹ (red region), Ė_kin ∼ 1 × 10⁴² erg s⁻¹ (blue region) and total Ė_kin = 2 × 10⁴² erg s⁻¹.

The kinetic coupling efficiencies are defined as

(4)

where L_bol is the AGN bolometric luminosity. ξ was estimated based on the AGN bolometric luminosities obtained from Brusa et al. (2022) (7.8 × 10⁴⁵ erg s⁻¹, LB22 hereafter) and Musiimenta et al. (2023) (1.6 × 10⁴⁶ erg s⁻¹, LM23 hereafter). We obtain total kinetic coupling efficiencies of ξ_LM23 = 0.01% and ξ_LB22  =  0.03%.

The ionised outflow momentum rate is given by

(5)

The total value of Ṗ_out is determined to be 4.9 × 10³⁴ g cm s⁻². The rate of AGN momentum output is obtained through Ṗ_AGN = L_bol/c, resulting in values of 2.6×10³⁵ g cm s⁻² (LB22) and 5.3 × 10³⁵ g cm s⁻² (LM23). The ratio Ṗ_out/Ṗ_AGN is ∼0.2 for both luminosities (LB22 and LM23) and outflow regions (red and blue). The results showing the different values for both red and blue regions and the total (where applicable) are summarised in Table 1.

There are several assumptions on different parameters such as geometry, n_e, temperature, the velocity of the outflowing gas and radius involved in computing outflow properties. These have been explored in several studies (e.g, Husemann et al. 2016; Fiore et al. 2017; Davies et al. 2020; Musiimenta et al. 2023; Hervella Seoane et al. 2023; Speranza et al. 2024). The largest uncertainty is due to the assumptions made on n_e that can change results by a couple of orders of magnitude as discussed in detail in Davies et al. (2020). Given that we assumed n_e = 500 cm⁻³, our results on the outflow properties can be considered as upper limits. Similarly, as discussed in Hervella Seoane et al. (2023), using different methods to measure the outflow velocity gives different results. For example, they reported higher values of mass outflow rates from parametric analysis using maximum outflow velocity than those derived from non-parametric analysis. Overall the values of outflow energetics in this study may have up to two orders of magnitude uncertainties due to the assumptions applied.

Finally, a comparison of outflow properties derived with and without PSF deblending for the few sources in Husemann et al. (2016) that have similar characteristics as our quasar (very small or low Hβ, spatial resolution, L [OIII] luminosity and redshift) showed that the mass outflow rate obtained from the PSF subtracted datacube was generally only a factor of 2–5 lower than those obtained from the total flux datacube.

4. Discussions

4.1. Interactions or merging in ID 608

The colour maps discussed in Sect. 3.1 and the kinematic analysis indicate that we may be dealing with a complex system with multiple dynamics occurring simultaneously. In Sect. 3.1, we present three galaxies (C1, C2, and C3) clearly associated with our quasar within ∼50 kpc from its nucleus. While C1 and C2 are both within 10 projected kpc from the nucleus (9.4 and 8.7 projected kpc away, respectively), C3 is instead at a farther distance of 47 projected kpc. The structure seen in the MSW, SW and FSW region extends from the quasar towards the position of C3 and appears to connect these two galaxies.

Although we are limited by the available data to establish the physical connection between these companion galaxies and ID 608, we hypothesise that ID 608 is embedded in an interacting system with three other galaxies due to the following reasons: 1) The extended structure in the MSW, SW and FSW region has no corresponding continuum emission and it is likely a cloud of ionised gas traced by the [O III] emission; 2) based on the kinematic analysis (see Sect. 3.3.2), this extended structure has low-velocity dispersion (W₈₀  <  600 km s⁻¹) and corresponding low velocities (< 500 km s⁻¹) which does not align with the presence of outflows. Therefore, this extended structure could be an ‘ionised cloud’ of gas in the form of a tidal tail which formed as a result of merging or gravitational interactions (post and ongoing) between ID 608 and C3. Moreover, there is gas rotation going on in this structure as seen in the colour images and the velocity peak map (with velocity gradient) which confirms our hypothesis. It seems that C3 may be fated to merge with the quasar or we are seeing a post interaction. This is a common feature of luminous red quasars showing prominent outflow signatures (e.g. Wylezalek et al. 2022). The right panel of Fig. 9 visually depicts this ongoing interaction.

Fig. 9. Cartoon showing the quasar (ID 608) and the deductions we have made based on analysis of the data, velocity maps and flux maps. The companion galaxies are marked. The outflows are shown in the northeast (NE) and southwest (SW) directions, which we have denoted on the diagram. The extended structure is indicated as an ionised cloud in red and blue colours that indicate the receding and approaching movement of the cloud as seen in Fig. 2. The possible interaction between the quasar and C3 is denoted by the black arrows. The distances are not to scale. There is no meaning attached to the sizes of the galaxies and outflows.

This discovery points to the fact that mergers trigger the growth of SMBH (Hopkins et al. 2008). Moreover, the morphology of this system explains the observed properties of quasars with outflows which include optical and X-ray obscuration and red colours. The complexity discovered in this obscured red quasar also highlights the importance of X-ray selection methods to select these kinds of objects that are buried in interacting systems.

4.2. Comparison of outflow properties with theoretical predictions

The mass outflow rates of ID 608 have been previously studied using the SDSS integrated spectrum (Brusa et al. 2022; Musiimenta et al. 2023). Using spatially resolved analysis, we estimated the mass outflow rates and energetics of ID 608 as presented in Sect. 3.4. We obtained total outflow mass of 27 × 10⁶ M_⊙, ionised outflow mass rate of 10 M_⊙ yr⁻¹, kinetic coupling efficiencies of ∼0.01%–0.03%. The momentum boosts are ∼0.2.

The mass outflow rate values obtained in this study using our spatially resolved data are generally higher than the values reported in Musiimenta et al. (2023) (∼ 2.6 M_⊙ yr⁻¹), derived under similar assumptions on n_e, T and geometry. This is because the methodology used to compute outflow properties such as luminosity, radius and velocity differ from the ones assumed previously and can bring discrepancies in the final values obtained. Moreover, the spatially resolved information available in this study, allowed us to compute a more reliable mass outflow rate by summing the different regions of the outflow instead of assuming only the redshifted broad component of the integrated spectrum. This highlights the power of spatially resolved spectroscopy to better characterise the kinematics and extension of the outflowing gas.

Our results are comparable to other ionised outflow properties in previous studies of obscured AGN at z  ∼  0.1 that have been computed using spatially resolved data. For example in Speranza et al. (2024) where a similar method (see Sect. 3.4) was used to compute the outflow properties, for a sample of four local type 2 quasars with spatially resolved outflows, they obtained ionised mass outflow rates of 3.3–6.5 M_⊙ yr⁻¹ within 3.1–12.6 kpc which are less than the values obtained in this study but within the same extension radius. Our values are within the ranges of 0.1–23 M_⊙ yr⁻¹ as ionised outflows in Musiimenta et al. (2023) at redshift z ∼ 0.5–1 and with the same AGN bolometric luminosity range as our quasar.

Musiimenta et al. (2023) reported a strong correlation between the mass outflow rate and the AGN bolometric luminosity, holding also for sources that show high outflow velocity but relatively low bolometric luminosities. The revised values of outflow parameters for ID 608 presented in this work are still fully consistent with the proposed correlation. ID 608 has high outflow velocity with a corresponding low ionised gas mass, further supporting the hypothesis of the presence of a fast-moving, less massive wind characteristic of the initial stages of the blow-out phase.

Moreover, as described in Sect. 4.1, ID 608 is possibly going through early-stage merger activity and was selected based on its red colours, X-ray brightness and obscuration. This indicates that these selection methods can be used to effectively isolate objects in the initial stages of outflowing phases that is, objects buried in dust enshrouded environments.

Although the direct comparison of observations results to theoretical prediction is not yet reliable, we note that our low values of the kinetic coupling efficiencies are not consistent with theoretical predictions of both energy-conserving mechanisms (0.5%–5%; Costa et al. 2014) and radiation-pressure driven scenario (0.1%–1%; Costa et al. 2018). Overall, the kinetic coupling efficiencies inferred from our results are an indication that the outflow observed is not very relevant from the energetic point of view. This study only considers ionised outflow and it could be that more massive outflows with high coupling efficiencies are being driven in the molecular or other phases. For example in Fiore et al. (2017), the molecular phase gives mass outflow rates that are higher than the ionised phase. This implies that the outflows in this quasar could be more relevant in other gas phases.

The momentum boosts obtained in this study are too low compared to predictions for either energy (∼20) or momentum (∼1) conserving mechanisms. With the uncertainties that arise from different assumptions in measuring outflow properties such as outflow luminosity, electron density, radius, and velocity, we can not draw firm conclusions about the driving mechanisms in this system. Observation from other wavebands is necessary to explore outflows of this source in other phases and the driving mechanisms for example radio and deeper X-rays.

We measure ionised mass loading factor using, MLF. We used the SFR, 2 M_⊙ yr⁻¹ from Brusa et al. (2022). This gives MLF of five which implies that the outflow plays a vital role in AGN feedback and can potentially affect the gas in the ISM faster than star formation and also an indication that the outflows observed are likely AGN-driven than SF-driven. Values of MLF that indicate such scenarios have been reported in the literature (e.g. Venturi et al. 2023; Speranza et al. 2024) for obscured quasars at low redshifts.

The improved resolution and sensitivity of MUSE data allow for more accurate measurements of the relevant parameters such as radius, velocity and [O III] luminosity. Despite the several assumptions made, the spatially resolved data enables improved estimates of the mass outflow rate, energetics and insights into the environment around the quasar.

5. Conclusions

In this study, we performed a detailed analysis of the morphology and kinematics of the [O III] emission line in object ID 608. The data used in this analysis were obtained through MUSE observations for program 108.22M6.001 period 108. Our findings are summarised as follows.

• Our source is part of an interacting system featuring two close-by companion galaxies and a more distant galaxy situated within ∼50 projected kpc away from the nucleus of the quasar as identified for the first time in this study. These companions include C1 which shows signatures of ionisation from the quasar, C2, and C3 with a spectrum that resembles that of an elliptical galaxy. The observed extended emission emanating from our source in the SW can be attributed to a cloud of gas, which is highly probable to have originated from the ongoing gravitational interactions or merging transpiring within the system.

• We have observed ionised outflows in our quasar with significant velocities extending up to – kpc. These outflows exhibit both blue- and red-shifted velocities, with the blue-shifted component reaching speeds as high as −1200 km s⁻¹ and the red-shifted component reaching speeds of up to 1000 km s⁻¹. For the first time, we have discovered the presence of a blueshifted counterpart to this outflow.

• Mass outflow rates were computed using the [O III] luminosity, assuming an electron density of 500 cm⁻³, estimating the radius and velocity. The resulting ionised mass outflow rates were determined to be 10 M_⊙ yr⁻¹, with a corresponding kinetic power of 2 × 10⁴² erg s⁻¹. The mass loading factor was estimated to be five which is an indication that these outflows are more likely to be AGN-driven than star formation-driven.

• The kinetic coupling efficiencies range from ∼0.01%–0.03%. The momentum boosts are ∼0.2. These values are low, indicating that the ionised outflow is not very relevant from the energetic point of view.

Our findings emphasise the importance of utilising integral field unit (IFU) and spatially resolved data to provide valuable insights into the morphology of outflows and better estimate the AGN outflow properties better. The morphology of the outflow is important to understand since it influences how they interact and couple with the ISM, whether it be radiatively driven or kinetic driven outflows. The interplay of interactions, mergers, and ionised gas within our system and the kinematics morphology renders this system interesting for further investigation in other wavelengths such as radio or deeper observation with Chandra to understand the nature of the companions. The discovery of this interactive nature around our quasar points to the importance of utilizing X-ray selection to probe outflows in red and obscured quasars which may be buried in such systems.

Acknowledgments

BM, GS, IEL, BL, IMR and CA are supported by the European Union’s Innovative Training Network (ITN) funded by the Marie Sklodowska-Curie Actions in Horizon 2020 No 860744 (BiD4BEST). MB acknowledges support from PRIN MIUR 2017PH3WAT (‘Black hole winds and the baryon life cycle of galaxies’). CRA acknowledges support from the projects ‘Feeding and feedback in active galaxies’, with reference PID2019-106027GB-C42, funded by MICINN-AEI/10.13039/501100011033 and ‘Quantifying the impact of quasar feedback on galaxy evolution’, with reference EUR2020-112266, funded by MICINN-AEI/10.13039/501100011033 and the European Union NextGenerationEU/PRTR. MP acknowledges support from the research project PID2021- 127718NB-I00 of the Spanish Ministry of Science and Innovation or State Agency of Research (MCIN/AEI/10.13039/501100011033). AL is partly supported by: ‘Data Science methods for MultiMessenger Astrophysics & Multi-Survey Cosmology’, Programmazione triennale 2021/2023 (DM n.2503 dd. 9 December 2019), Programma Congiunto Scuole; Fondazione ICSC, Spoke 3 Astrophysics and Cosmos Observations Project ID CN-00000013; INAF Large Grant 2022 funding scheme, project MeerKAT and LOFAR Team up: a Unique Radio Window on Galaxy or AGN co-Evolution. AL thanks M.V. Zanchettin for useful discussions.

References

Appendix A: Velocity channels considering 100 km s⁻¹ bin size.

The extracted flux maps from various velocity bins of 100 km s⁻¹ from −1200 km s⁻¹ to 1200 km s⁻¹ are shown in Fig. A.1. In this Figure, the different structures towards the NE and SW directions are visible from bins [−600, −500] km s⁻¹ to [200, 300] km s⁻¹. The rest of the bins have the [O III] emission not extending far away from the nucleus of the quasar.

Fig. A.1. Velocity channels considering smaller bins of 100 km s−1 on logarithmic scale. The white contours indicate the 2σ and 3σ.

Appendix B: Flux maps and velocity maps from three-Gaussian fitting and Voronoi binning

B.1. Three-Gaussian fitting without Voronoi binning

We performed spaxel-by-spaxel fitting with a maximum of three Gaussian components (n = 1,2,3) around [O III] doublet ([O III]λ5007 and [O III]λ4949) using the lmfit Python package and the in-house script described in Speranza et al. (2024). The three-Gaussian components are tuned in a such way as to avoid underfitting in the central regions and overfitting in the outer regions of the quasar. In this method, more than one Gaussian is fitted when the χ² improves more than 10% with respect to the previous fit. It is worth noting that the parameters of both emission lines were tied together and the flux ratios between the [OIII]λ4959 and [OIII]λ5007 lines were fixed at a ratio of 1:3 (Osterbrock 1981). We set two narrow components with velocity dispersions (σ) less than 500 km s⁻¹ and one broader component with σ greater than 500 km s⁻¹. The velocity dispersions of each kinematic component can vary within these specified limits. The resulting flux map and kinematic maps are shown in Fig. B.1. These results are similar to the ones presented in Sect. 3.3.2.

B.2. Voronoi bins data fitting with three-Gaussian

Another approach to visualising the kinematics and morphology of outflowing gas in AGNs is by grouping spaxels into bins. A commonly used method for this is adaptive spatial binning using Voronoi tessellations, which was introduced by Cappellari & Copin (2003). The goal of this method is to achieve a constant signal-to-noise ratio (S/N) per bin, which is crucial for accurately analysing IFU data. Compared to analysing spaxel-by-spaxel, this approach allows for more efficient and reliable extraction of kinematic information. Analysing each spaxel separately can be computationally intensive and can produce noisy and biased results if the S/N is not sufficiently high. In contrast, adaptive spatial binning provides a constant S/N per bin, minimising the variance of the S/N ratio within each bin and enabling more reliable and efficient extraction of kinematics.

To implement this method, the desired target S/N for each bin and the minimum S/N required for a spaxel to be included in a bin must be set. Previous studies have employed this approach with various target S/N values, such as 10, 12, 15, 30, 40, and 50, for different analyses, including gas kinematics (e.g. Girdhar et al. 2022; Speranza et al. 2022; Olivares et al. 2022; D’Eugenio et al. 2023; Zheng et al. 2023; Kolcu et al. 2023), stellar populations (Fetherolf et al. 2020, and referenced therein), and resolved optical properties of high-redshift galaxies and others (Duncan et al. 2023; Zhu et al. 2023). We have utilised this approach on our data cube to enhance the S/N. We binned our data cubes to achieve an average S/N of five per wavelength channel in each bin. Spaxels with S/N below three are discarded. We chose a relatively low threshold compared to previous studies to account for the low S/N of our data and to avoid excessively large bins.

After the binning process, we applied three-Gaussian fitting procedures to the [O III]λ4959 and [O III]λ5007 doublet, similar to the one described in Appendix B.1 for each Voronoi bin. The Voronoi bin map and the resulting flux and velocity maps are shown in Fig. B.2.

We compared our results from the three Gaussian fitting with Voronoi binning and without. The results for the velocity profiles are more or less similar giving the same velocity ranges. However, the velocity values in Voronoi binning maps are higher in the SW extended structure with V₀₅ ∼-900 km s⁻¹ and V₀₅ ∼700 km s⁻¹. This also has corresponding high value of velocity dispersion (W₈₀) > 800 km s⁻¹. The Voronoi bin size in this region is big because of the low signal-to-noise (S/N) and as a consequence the high kinematic values. The kinematics of our data are well-defined without the binning. Therefore the results in this study are based on the multi-Guassian fitting described in Sect. 3.3.

Fig. B.1. Results from the spaxel-by-spaxel three-Gaussian fitting of 8″×7.5″≅ 53 kpc × 50 kpc region around the quasar. Top left panel: the flux map. Top right panel: the velocity peak. In the middle panel, we present the V05 and V95 maps indicating the velocities at the 5th and 95th percentile, respectively. Bottom left panel: the line width (W80). Bottom right panel: the number of Gaussians fitted in each pixel.

Fig. B.2. Results from three-Gaussian fitting of Voronoi bins. Top left panel: The flux map. Top right panel: the velocity peak. In the middle panel, we present the V05 and V95 maps indicating the velocities at the 5th and 95th percentile, respectively. Bottom left panel: The line width (W80). Bottom right panel: the Voronoi bins.

All Tables

All Figures

Fig. 1. Full data cube over the whole wavelength range covered in our MUSE observation and the continuum and [O III] emission for size-reduced cubes. Panel a (left) is the image of the full MUSE’s FOV (1′×1′≅402 kpc × 402 kpc) obtained after data reduction. ID 608 is marked by a purple circle. Panel b (centre) shows an inner zoom of the cube (30″ × 30″ ≅ 201 kpc × 201 kpc) with the continuum in blue, extracted from two regions of the spectrum (around [7825, 7900] Å and [8100, 8150] Å observed wavelengths) and the [O III] emission in red extracted from observed wavelengths 7890 Å to 8100 Å. Panel c (right) shows a further zoom (8″ × 7.5″ ≅ 53 kpc × 50 kpc) of the [O III] emission centred on the quasar position. The position of the nucleus of the quasar obtained from the maximum of the [OIII] emission is marked by the red star. The contours in both panels b and c represent 3σ, 5σ, 15σ, 30σ, and 60σ. North is up, East is left.

In the text

Fig. 2. RGB colour image, velocity channels and velocity distribution for the [O III] emission obtained from data without fitting. The different panels are named as a, b, c, d, e, f, and g. Panels a–e: [O III] intensity maps extracted in different velocity channels, as labelled. The arrows indicate the position of companion one (C1). Panel f: RGB colour image constructed from three [O III] velocity bins ([−250, −50] km s−1, [−50, −50] km s−1 and [50, 250] km s−1). Panel g shows the velocity distribution map attained without fitting to show the different regions as will be referred to in this study. The regions are NE, NW, SW, MSW, and FSW. The star shows the position of the nucleus of the quasar. The figures in all the panels are at the same scale and extracted from 8″ × 7.5″ ≅ 53 kpc × 50 kpc FOV.

In the text

	Fig. 3. Spectra of the quasar and the companion galaxies C1, C2, and C3 with their estimated redshifts. The spectrum in red indicates the quasar (ID 608), companion galaxy one in black, companion galaxy two in yellow, and companion galaxy three in blue.
In the text

Fig. 4. Colour map of the continuum emission overlapped with the [O III] emission in 30″ × 30″ ≅ 201 × 201 kpc FOV. The continuum emission is shown in blue with corresponding white contours and the [O III] in red (no contours). The companion galaxies are marked by yellow stars, while ID 608 nucleus is marked by a red star. In the bottom left of this figure, we show the zoom-in map of the continuum (blue contours), [O III] (red without contours), and Hβ (green contours) emission. The positions of C1 and C2 are indicated that do not appear to emit in Hβ. The contours in the zoom in represent 3σ, 5σ, 15σ, 30σ, and 60σ.

In the text

Fig. 5. Spectra relative from the six regions shown in Fig. 2. The spectra are extracted from a small area of a 0.6″ × 0.6″ box at the given position in the FOV similar to Fig. 2f. Within the six regions shown, only the nuclear region and NE are best fitted by more than two Gaussian functions. One or two Gaussian components can adequately describe all the remaining spectra. The legend for all the spectra is represented in the spectrum from the NE region where the black solid spectrum shows the data, the red dotted shows the best fit, the grey dotted spectrum shows the residual and all the other colours (blue, magenta, solid red, green) show the different components fitted. The three vertical lines highlight the centres of the Hβλ4861, [O III]λ4949 and [O III]λ5007 emission lines. The yellow stars indicate the positions of C1 and C2. In the bottom right of the map, the white arrow points towards the position of C3 that is outside this map size. The different regions that is NE, NW, MSW, SW, FSW and ID 608 are indicated in the respective spectrum.

In the text

Fig. 6. Residual spectrum and map after PSF subtraction show that the [O III] emission in our quasar is resolved. The left panel shows the residual spectrum analysis of a randomly selected pixel. The two vertical lines serve as markers highlighting the wavelengths of [O III]λ5007 and [O III]λ4949 emission lines for the centre of ID 608. The right panel showcases a residual map resulting from collapsing the residual spectrum of each pixel. The red square on the right panel corresponds to the specific pixel from which the residue spectrum was extracted in the left panel and the centre is marked by a white star.

In the text

Fig. 7. Results from the spaxel-by-spaxel multi-Gaussian fitting of 8″ × 7.5″ ≅ 53 kpc × 50 kpc region around the quasar. Top left panel: the flux map. Top right panel: the velocity peak. In the middle panel, we present the V05 and V95 maps indicating the velocities at the 5th and 95th percentile, respectively. Bottom left panel: the line width (W80). Bottom right panel: the number of Gaussians fitted in each pixel.

In the text

Fig. 8. Flux maps for the blue-shifted and red-shifted velocities from the multi-Gaussian fitting. We select the region [−2200, −250] km s−1 to extract the flux map for negative velocities (the left panel). In the right panel, the flux map for the positive velocities [250, 2000] km s−1. The light blue, yellow, and green contours indicate the 4, 3, and 2σ levels, respectively. These contours help to determine the significance level of the measured flux values. The blue dashed line represents the maximum distance to the 3σ fitted ellipse in white.

In the text

Fig. 9. Cartoon showing the quasar (ID 608) and the deductions we have made based on analysis of the data, velocity maps and flux maps. The companion galaxies are marked. The outflows are shown in the northeast (NE) and southwest (SW) directions, which we have denoted on the diagram. The extended structure is indicated as an ionised cloud in red and blue colours that indicate the receding and approaching movement of the cloud as seen in Fig. 2. The possible interaction between the quasar and C3 is denoted by the black arrows. The distances are not to scale. There is no meaning attached to the sizes of the galaxies and outflows.

In the text

	Fig. A.1. Velocity channels considering smaller bins of 100 km s−1 on logarithmic scale. The white contours indicate the 2σ and 3σ.
In the text

Fig. B.1. Results from the spaxel-by-spaxel three-Gaussian fitting of 8″×7.5″≅ 53 kpc × 50 kpc region around the quasar. Top left panel: the flux map. Top right panel: the velocity peak. In the middle panel, we present the V05 and V95 maps indicating the velocities at the 5th and 95th percentile, respectively. Bottom left panel: the line width (W80). Bottom right panel: the number of Gaussians fitted in each pixel.

In the text

	Fig. B.2. Results from three-Gaussian fitting of Voronoi bins. Top left panel: The flux map. Top right panel: the velocity peak. In the middle panel, we present the V05 and V95 maps indicating the velocities at the 5th and 95th percentile, respectively. Bottom left panel: The line width (W80). Bottom right panel: the Voronoi bins.
In the text