# GA-NIFS: A galaxy-wide outflow in a Compton-thick mini-broad-absorption-line quasar at $z = 3.5$ probed in emission and absorption

Michele Perna, Santiago Arribas, Xihan Ji, Cosimo Marconcini, Isabella

Lamperti, Elena Bertola, Chiara Circosta, Francesco d'Eugenio, Hannah

Übler, Torsten Böker, et al.

► To cite this version:

Michele Perna, Santiago Arribas, Xihan Ji, Cosimo Marconcini, Isabella Lamperti, et al. GA-NIFS: A galaxy-wide outflow in a Compton-thick mini-broad-absorption-line quasar at $z = 3.5$ probed in emission and absorption. Astronomy & Astrophysics - A&A, 2025, 694, 10.1051/0004-6361/202453090. insu-05101707

HAL Id: insu-05101707

https://insu.hal.science/insu-05101707v1

Submitted on 7 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# GA-NIFS: A galaxy-wide outflow in a Compton-thick mini-broad-absorption-line quasar at $z = 3.5$ probed in emission and absorption

Michele Perna<fnref n="1" />,*, Santiago Arribas<fnref n="1" />, Xihan Ji<fnref n="2" />,<fnref n="3" />, Cosimo Marconcini<fnref n="4" />,<fnref n="5" />, Isabella Lamperti<fnref n="4" />,<fnref n="5" />,<fnref n="1" />, ∗
Elena Bertola<fnref n="5" />, Chiara Circosta<fnref n="6" />, Francesco D'Eugenio<fnref n="2" />,<fnref n="3" />, Hannah Übler<fnref n="7" />,<fnref n="2" />,<fnref n="3" />, Torsten Böker<fnref n="8" />,
Roberto Maiolino<fnref n="2" />,<fnref n="3" />,<fnref n="10" />, Andrew J. Bunker<fnref n="9" />, Stefano Carniani<fnref n="11" />, Stéphane Charlot<fnref n="12" />, Chris J. Willott<fnref n="13" />,
Giovanni Cresci<fnref n="5" />, Alessandro Marconi<fnref n="4" />,<fnref n="5" />, Eleonora Parlanti<fnref n="11" />, Bruno Rodríguez Del Pino<fnref n="1" />, ∗
Jan Scholtz<fnref n="2" />,<fnref n="3" />, and Giacomo Venturi<fnref n="11" />, †

Received 20 November 2024 / Accepted 18 December 2024

## ABSTRACT

**Context.** Studying the distribution and properties of ionised gas in outflows driven by active galactic nuclei (AGN) is crucial for understanding the feedback mechanisms at play in extragalactic environments. These outflows provide key insights into the regulation of star formation and the growth of supermassive black holes.

**Aims.** In this study, we explore the connection between ionised outflows traced by rest-frame ultra-violet (UV) absorption and optical emission lines in GS133, a Compton thick AGN at $z = 3.47$. We combine observations from the *James Webb* Space Telescope (JWST) NIRSpec Integral Field Spectrograph (IFS) with archival Very Large Telescope (VLT) VIMOS long-slit spectroscopic data, as part of the 'Galaxy Assembly with NIRSpec IFS' (GA-NIFS) project.

**Methods.** We performed a multi-component kinematic decomposition of the UV and optical line profiles to derive the physical properties of the absorbing and emitting gas in GS133.

**Results.** Our kinematic decomposition reveals two distinct components in the optical emission lines. The first component likely traces a rotating disc with a dynamical mass of $2 \times 10^{10} M_{\odot}$. The second component corresponds to a galaxy-wide, bi-conical outflow, with a velocity of $\sim \pm 1000 \text{ km s}^{-1}$ and an extension of $\sim 3 \text{ kpc}$. The UV absorption lines show two outflow components, with bulk velocities $v_{\text{out}} \sim -900 \text{ km s}^{-1}$ and $\sim -1900 \text{ km s}^{-1}$, respectively. This characterises GS133 as a mini-broad absorption line (mini-BAL) system. Balmer absorption lines with similar velocities are tentatively detected in the NIRSpec spectrum. Both photoionisation models and outflow energetics suggest that the ejected absorbing gas is located at 1–10 kpc from the AGN. We use 3D gas kinematic modelling to infer the orientation of the [O III] bi-conical outflow, and find that a portion of the emitting gas resides along our line of sight, suggesting that [O III] and absorbing gas clouds are partially mixed in the outflow. The derived mass-loading factor (i.e. the mass outflow rate divided by the star formation rate) of 1–10, and the kinetic coupling efficiency (i.e. the kinetic power divided by $L_{\text{AGN}}$) of 0.1–1% suggest that the outflow in GS133 provides significant feedback on galactic scales.

**Key words.** galaxies: active – galaxies: high-redshift – quasars: absorption lines – quasars: supermassive black holes

## 1. Introduction

Accretion-disc outflows from active galactic nuclei (AGN) are thought to provide significant kinetic-energy feedback to their host galaxies, potentially stopping star formation and suppressing accretion onto the central black hole (e.g. Di Matteo et al. 2005; Hopkins et al. 2008; Perna et al. 2018; Giustini & Proga 2019; Bertola et al. 2024).

The effectiveness of AGN feedback depends on the efficiency with which these outflows couple to the interstellar medium (ISM) in the host galaxies, an issue that remains complex and not fully understood (e.g. Harrison et al. 2017).

The study of AGN outflows often relies on the analysis of blueshifted broad absorption lines (BALs) tracing ionised gas in rest-frame ultraviolet (UV) spectra. These features are

<footnote n="1">¹ Centro de Astrobiología (CAB), CSIC-INTA, Cra. de Ajalvir Km. 4, 28850 Torrejón de Ardoz, Madrid, Spain</footnote>

<footnote n="2">² Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="3">³ Cavendish Laboratory – Astrophysics Group, University of Cambridge, 19 JJ Thomson Avenue, Cambridge CB3 0HE, UK</footnote>

<footnote n="4">⁴ Università di Firenze, Dipartimento di Fisica e Astronomia, Via G. Sansone 1, 50019 Sesto F.no, Firenze, Italy</footnote>

<footnote n="5">⁵ INAF – Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy</footnote>

<footnote n="6">⁶ European Space Agency, ESAC, Villanueva de la Cañada, E-28692 Madrid, Spain</footnote>

<footnote n="7">⁷ Max-Planck-Institut für Extraterrestrische Physik (MPE), Gießenbachstraße 1, 85748 Garching, Germany</footnote>

<footnote n="8">⁸ European Space Agency, c/o STScI, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="9">⁹ Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="10">¹⁰ Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="11">¹¹ Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy</footnote>

<footnote n="12">¹² Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 Bis bd Arago, 75014 Paris, France</footnote>

<footnote n="13">¹³ National Research Council of Canada, Herzberg Astronomy & Astrophysics Research Centre, 5071 West Saanich Road, Victoria, BC V9E 2E7, Canada</footnote>

<footnote>* Corresponding author; mperna@cab.inta-csic.es</footnote>
typically defined by having a full width at half maximum (FWHM) >3000 km s⁻¹ and velocity shifts greater than 5000 km s⁻¹ in the AGN host galaxy rest frame (e.g. Weymann et al. 1981). Mini-BALs, narrower than BALs but still broad enough to be identified as AGN outflows, have velocity shifts between 500 and 2000 km s⁻¹ (e.g. Hamann & Sabra 2004; Moravec et al. 2017; Hamann et al. 2019; Maiolino et al. 2024a). Even slower outflows are observed in the form of narrow absorption lines (NALs), with FWHM = 100–500 km s⁻¹ and velocity offsets of a few 100 km s⁻¹ (e.g. Vestergaard 2003; Hamann & Sabra 2004). There is no consensus about the spatial extent of UV absorption line outflows, as they trace gas along our LOS, making direct measurements of their extent impractical. Indirect methods, such as photoionisation analysis, indicate that UV outflows can reach distances ranging from tens to even thousands of parsecs (Arav et al. 2013; Xu et al. 2020). Combined, the NAL and mini-BAL outflows are more common than the well-studied BALs in AGN spectra. Specifically, ~50% of AGN contain either NAL or mini-BAL outflows, whereas only ~20% contain a C IV BAL outflow (e.g. Giustini & Proga 2019 and references therein).

Another approach to studying AGN outflows involves observations of rest-frame optical emission lines, which may reveal blueshifted broad wings in ionized gas tracers such as [N II] λ6583 and, more frequently, the [O III] λ5007 forbidden lines (e.g. Arribas et al. 2014; Genzel et al. 2014; Perna et al. 2015, 2019; Perrotta et al. 2019; Förster Schreiber et al. 2019; Kakkad et al. 2020; Tozzi et al. 2024). Since forbidden features can originate only from low electron density regions (i.e. $n_e < 10^6$ cm⁻³), their broadening cannot be explained by any contamination from high density AGN broad line regions (BLR). Therefore, forbidden lines are an excellent tracer of ionised outflows, on scales extending up to several kiloparsecs from the AGN (e.g. Harrison et al. 2012; Liu et al. 2013; Cresci et al. 2023; Ulivi et al. 2024). A large fraction of AGN presents signatures of outflows in their [O III] phase, from ~20% to ~70%, depending on the selection criteria and outflow velocity thresholds used (e.g. Véron-Cetty et al. 2001; Mullaney et al. 2013; Perna et al. 2017; Musiimenta et al. 2023).

More generally, ejections of material from the inner regions up to the host galaxy scale can involve different gas phases, and are revealed as parsec-scale BALs (e.g. Bruni et al. 2019; Vietri et al. 2022) and ultrafast outflows (up to ~30% of the speed of light) detected in X-rays (e.g. Pounds et al. 2003; Chartas et al. 2021; Matzeu et al. 2023), to kiloparsec-scale outflows observed in molecular, neutral, and ionised gas (e.g. Cicone et al. 2014; Lamperti et al. 2022; Mehdipour et al. 2023; Parlanti et al. 2024a; Ulivi et al. 2025). Therefore, a comprehensive characterisation of the outflow phenomena requires the use of a range of major facilities that work at different wavelengths and angular resolutions.

The need of multiple observations with different facilities have affected the direct comparison of outflows identified with different tracers, like the kinematically disturbed UV absorption- and optical emission-line components. It is not yet determined whether these phenomena are entirely distinct and are caused during different phases in the AGN fuelling and outflow life cycle, or whether they represent the same outflow episode observed with different methods, possibly depending on the inclination of the systems relative to the line of sight (LOS). In fact, individual absorption lines only reveal the outflow components along the LOS to the AGN; in contrast, optical emission lines can trace outflow emission in any direction. The literature includes only a limited number of studies that compare absorp-

tion and emission lines across AGN samples (e.g. Schulze et al. 2017; Xu et al. 2020; Ahmed et al. 2024; Temple et al. 2024; Kehoe et al. 2024).

Analyses of UV absorption and optical forbidden emission lines in individual AGN sources are also rare: there is only one integral field spectroscopy (IFS) study of two bright AGN at $z \lesssim 0.1$ (Liu et al. 2015), and long-slit works for three individual quasars at $z \approx 2$ (Tian et al. 2019; Dai et al. 2024; Stepney et al. 2024). The use of IFS observations is particularly important as they allow mapping of the ejected [O III] gas in two spatial and one velocity dimension. This capability enables direct measurement of the spatial extent of the outflows, which can only be indirectly inferred from UV absorption line analyses. In turn, this information is essential to quantify the impact these outflows have on the host galaxy.

In this paper, we present the *James Webb* Space Telescope (JWST) Near Infrared Spectrograph (NIRSpec) IFS and Very Large Telescope (VLT) VIMOS observations of GS133, a highly obscured AGN at $z = 3.5$. The NIRSpec IFS data covers the rest-frame optical ionised gas, while the VIMOS long-slit spectrum covers the rest-frame UV features. Therefore, the combination of these datasets allows us to characterise both UV and optical outflows. GS133 is located at 03:32:04.938, −27:44:31.73 in the GOODS-South field (Giavalisco et al. 2004), and is classified as an X-ray AGN in the 7 Ms exposure of the Chandra Deep Field-South (target #133 in Luo et al. 2017). Both X-ray spectroscopy fitting analysis (Luo et al. 2017) and mid-infrared diagnostics (Guo et al. 2021) identify this source as a Compton thick (CT) AGN, having a column density $N_H > 1.5 \times 10^{24}$ cm⁻² (Comastri 2004). It has an absorption-corrected 2–10 keV luminosity $L_X = 9 \times 10^{43}$ erg s⁻¹, an integrated luminosity from UV to infrared $L_{bol} = 4.9 \times 10^{45}$ erg s⁻¹, a black hole mass $M_{BH} = 2 \times 10^7$ M☉, and an Eddington ratio $\lambda_{Edd} = L_{bol}/L_{Edd} = 2.4$ (Guo et al. 2020). Note that because GS133 is a highly obscured AGN with no BLR emission, the BH mass was estimated assuming the $M_{BH}-M_*$ relation with a scaling of 0.003 (Guo et al. 2021).

The paper is outlined as follows. In Section 2 we describe our JWST/NIRSpec IFS observations and data reduction, as well as the archival VLT/VIMOS data. Detailed data analysis of the optical and UV integrated spectra are reported in Sect. 3, while the NIRSpec spatially resolved spectroscopic analysis is reported in Sect. 4. Multi-wavelength diagnostics are discussed in Sect. 5. The calculation of the dynamical mass of the GS133 host galaxy is provided in Sect. 6, followed by an analysis of the outflow energetics in Sect. 7, which incorporates photoionisation and 3D gas kinematic models. The GS133 environment is examined in Sect. 8. A summary of our findings is presented in Sect. 9. Throughout, we adopt a Chabrier (2003) initial mass function (0.1–100 $M_☉$) and a flat $\Lambda$CDM cosmology with $H_0 = 70$ km s⁻¹ Mpc⁻¹, $\Omega_\Lambda = 0.7$, and $\Omega_m = 0.3$. In our analysis of NIRSpec data, we use vacuum wavelengths according to their calibration. However, when discussing rest-frame optical emission lines, we quote their air wavelengths. Similarly, for VIMOS data, we use the air frame based on their wavelength calibration, but refer to rest-frame UV lines using their vacuum wavelengths, consistent with the standard pre-JWST practice.

## 2. Observations and data processing

### 2.1. NIRSpec data

GS133 was observed on September 12, 2022, under programme #1220 (PI: N. Lützgendorf), as part of the NIRSpec IFS GTO programme 'Galaxy Assembly with NIRSpec IFS'
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Integrated spectrum of GS133. The blue curve identifies the VLT/VIMOS spectrum; the black line shows the JWST/NIRSpec spectrum, integrated over a circular aperture of $r = 0.5''$. The most prominent emission lines are marked with grey vertical lines. The gap in the middle of the NIRSpec spectrum ($\lambda_{\text{rest}} \sim 5400 \AA$) is due to the separation between the two NIRSpec detectors.</figcaption></figure>

(GA-NIFS, e.g. Rodríguez Del Pino et al. 2024; Scholtz et al. 2024; Lamperti et al. 2024). The project is based on the use of the NIRSpec’s IFS mode, which provides spatially resolved spectroscopy over a contiguous 3.1” × 3.2” sky area, with a sampling of 0.1”/spaxel and a comparable spatial resolution (Böker et al. 2022; Rigby et al. 2023). The IFS observations were taken with the grating/filter pair G235H/F170LP. This results in a data cube with spectral resolution $R \sim 2700$ over the wavelength range 1.7–3.1 µm (Jakobsen et al. 2022). The observations were taken with the NRSIRS2RAPID readout pattern (Rauscher et al. 2017) with 60 groups per integration and one integration per exposure, using a 4-point medium cycling dither pattern, resulting in a total exposure time of 3560 seconds.

We used v1.8.2 of the JWST pipeline with CRDS context 1105 to create a final cube with drizzle weighting. A patch was included to correct some important bugs that affect this specific version of the pipeline (see details in Perna et al. 2023a). We corrected count-rate images for 1/f noise through a polynomial fit. During stage 2, we removed all data in regions of known failed open MSA shutters. We also masked pixels at the edge of the slices (two pixel wide) to conservatively exclude pixels with unreliable sflat corrections, and implemented the outlier rejection of D'Eugenio et al. (2024). The combination of a dither and drizzle weighting allowed us to sub-sample the detector pixels, resulting in cube spaxels of 0.05”.

## 2.2. VIMOS data

The rest-frame UV spectrum of GS133 was observed as part of the survey 'VANDELS: a VIMOS survey of the CDFS and UDS fields' (McLure et al. 2018). This ESO public spectroscopic survey, conducted with the VIMOS spectrograph on the VLT, aimed to obtain ultra-deep, medium resolution, red-optical spectra of ~2000 galaxies at $1 < z < 7$.

GS133 was observed on May 28, 2017, using the MR grism and the GG475 filter, with a 1” slit width and a slit length of 10” oriented east-west on the sky. The total exposure time was 41 hours. This setup provides a wavelength coverage of 4800–10 000 Å with a nominal resolution $R = 580$, corresponding to a velocity resolution of approximately 500 km s⁻<fnref n="1" />. The spectroscopic analysis presented in this work takes advantage of the VANDELS DR4 fully reduced spectra (Garilli et al. 2021), downloaded from the ESO archive<fnref n="1" />.

A visual inspection of the 2-dimensional spectrum does not reveal any spatially resolved structure. Therefore, the analysis presented in this work focusses on the one-dimensional (1D) VIMOS spectrum.

## 3. Spatially integrated spectra

Figure 1 shows the integrated VIMOS and NIRSpec spectra of GS133. The VIMOS spectrum (blue curve) covers the 1070–2300 Å rest-frame wavelength range, while the NIRSpec spectrum (black curve) covers the range between 3710 and 7085 Å rest-frame; it was extracted from a circular aperture centred at the position of the AGN, with a radius $r = 0.5''$, which matches the width of the VIMOS slit. In this section, we analyse these spectra, and infer the integrated properties of the outflows both in the UV range, which is spatially unresolved in VIMOS data, and in optical, making use of NIRSpec IFS observations. Instead, Sect. 4 investigates the spatially resolved properties of the optical emission.

### 3.1. Optical emission lines in NIRSpec cube

We fitted the most prominent gas emission lines by using the Levenberg-Marquardt least-squares fitting code CAP-MPFIT (Cappellari 2017). In particular, we modelled the $H\alpha$ and $H\beta$ lines, the [O III] $\lambda\lambda4959, 5007$, [N II] $\lambda\lambda6548, 83$, and [S II] $\lambda\lambda6716, 31$ doublets with a combination of Gaussian profiles, applying a simultaneous fitting procedure, so that all line features of a given kinematic component have the same velocity centroid and FWHM (e.g. Perna et al. 2020). We used rest-frame vacuum wavelengths. Moreover, the relative flux of the two [N II] and [O III] doublet components was fixed to 2.99 and the [S II] flux ratio was required to be within the range $0.44 < f(\lambda6716)/f(\lambda6731) < 1.42$ (Osterbrock & Ferland 2006). The final number of kinematic components used to model the spectra was derived on the basis of the Bayesian information criterion (BIC; Schwarz 1978).

Figure 2 shows the best-fit model around the $H\beta$-[O III] and $H\alpha$-[N II] regions. All emission lines exhibit a narrow core, and prominent blue and red wings. Modelling these complex profiles required four Gaussian components. The Gaussian component associated with the line cores likely traces the systemic, unperturbed emission in the GS133 host galaxy, with a broadening (velocity dispersion $\sigma = 137 \pm 5$ km s⁻<fnref n="1" />) attributed solely to the gravitational potential of the system. In contrast, the other

<footnote n="1">¹ https://archive.eso.org/dataset/ADP.2021-02-01T16:09:36.529</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** JWST/NIRSpec integrated spectrum of GS133. The black curve identifies the integrated spectrum; the total, multi-component best-fit curve is in red, while all individual Gaussian components are shown with different colours. The fit residuals are reported in the top panels. The most prominent emission lines are marked with grey vertical lines. Light and dark yellow shaded areas in the right panel mark the [O III] emission within ±2000 km s⁻¹ and ±1000 km s⁻¹, respectively, from the GS133 systemic redshift.</figcaption></figure>

<table><caption>**Table 1.** Properties of GS133 inferred from the optical emission lines of the NIRSpec integrated spectrum, obtained for the systemic ('sys'), outflow ('out'), and total ('tot') line profiles.</caption><thead><tr><th>Measurement</th><th>Value</th></tr></thead><tbody><tr><td>L<sub>[OIII], sys</sub></td><td>(1.21 ± 0.02) × 10<sup>43</sup> erg s<sup>-1</sup></td></tr><tr><td>L<sub>Hα, sys</sub></td><td>(9.97 ± 0.13) × 10<sup>42</sup> erg s<sup>-1</sup></td></tr><tr><td>A<sub>V, sys</sub></td><td>1.3<sup>+0.2</sup><sub>-0.4</sub> mag</td></tr><tr><td>n<sub>e, sys</sub></td><td>3300 ± 500 cm<sup>-3</sup></td></tr><tr><td>z<sub>sys</sub></td><td>3.4732 ± 0.0001</td></tr><tr><td>σ<sub>sys</sub></td><td>137 ± 5 km s<sup>-1</sup></td></tr><tr><td>L<sub>[OIII], out</sub></td><td>1.05<sup>+0.01</sup><sub>-0.04</sub> × 10<sup>43</sup> erg s<sup>-1</sup></td></tr><tr><td>L<sub>Hα, out</sub></td><td>(3.4 ± 0.3) × 10<sup>42</sup> erg s<sup>-1</sup></td></tr><tr><td>A<sub>V, out</sub></td><td>1.3<sup>+0.6</sup><sub>-0.9</sub></td></tr><tr><td>δv<sub>out</sub></td><td>-150 ± 6 km s<sup>-1</sup></td></tr><tr><td>σ<sub>out</sub></td><td>515 ± 3 km s<sup>-1</sup></td></tr><tr><td>W80<sub>out</sub></td><td>1320 ± 60 km s<sup>-1</sup></td></tr><tr><td>v10<sub>out</sub></td><td>-890 ± 40 km s<sup>-1</sup></td></tr><tr><td>v90<sub>out</sub></td><td>425 ± 10 km s<sup>-1</sup></td></tr><tr><td>L<sub>[OIII], tot</sub></td><td>(2.26 ± 0.01) × 10<sup>43</sup> erg s<sup>-1</sup></td></tr><tr><td>L<sub>Hα, tot</sub></td><td>(1.40<sup>+0.1</sup><sub>-0.09</sub>) × 10<sup>43</sup> erg s<sup>-1</sup></td></tr><tr><td>A<sub>V, tot</sub></td><td>1.1<sup>+0.4</sup><sub>-0.6</sub> mag</td></tr><tr><td>Δv<sub>tot</sub></td><td>-85 ± 4 km s<sup>-1</sup></td></tr><tr><td>σ<sub>tot</sub></td><td>370 ± 3 km s<sup>-1</sup></td></tr><tr><td>W80<sub>tot</sub></td><td>780 ± 60 km s<sup>-1</sup></td></tr><tr><td>v10<sub>tot</sub></td><td>-530 ± 30 km s<sup>-1</sup></td></tr><tr><td>v90<sub>tot</sub></td><td>250 ± 20 km s<sup>-1</sup></td></tr></tbody></table>

**Notes.** Emission line luminosities are not corrected for dust extinction. All velocity measurements refer to the [O III] line.

components likely trace a powerful outflow. Therefore, to characterise the outflow, we derived its properties from the integrated spectrum by considering the entire reconstructed profile, combining all but the narrowest and most prominent Gaussian component. This approach is justified by the potential blending of distinct (approaching and receding) outflow components that contribute to the overall profile (see Sect. 4).

As the measurement of kinematic properties of emission line profiles can be obtained using different methods, depending on the specific needs and characteristics of the data,

in Table 1 we reported both non-parametric velocities (e.g. Zakamska & Greene 2014), obtained by measuring the velocity *v* at which a given fraction of the total best-fit line flux is collected using a cumulative function, and the velocity offset Δ*v* and velocity dispersion *σ* defined as the moment-1 and moment-2 measurements, respectively. Table 1 presents all velocity measurements inferred from the [O III] profile, for the reconstructed narrow and outflow components, and for the best-fit total emission line profile. Specifically, we reported the 10th-percentile (*v*10) and the 90th-percentile (*v*90) velocities of the fitted line profile, associated with the highest blueshifted and redshifted gas emission, respectively; W80, defined as *v*90–*v*10 and representing the line width; and the velocity shift and velocity dispersion.

We derived the GS133 redshift from the measured wavelength of the narrow Gaussian reproducing the [O III] core emission in the integrated spectrum shown in Fig. 2: *z* = 3.4732 ± 0.0001. This value is roughly in agreement with previous estimates *z* = 3.462 (Cristiani et al. 2000) and *z* = 3.4739 (Balestra et al. 2010) obtained from (noisier and lower spectral resolution) ground-based rest-frame UV spectra.

The [O III] perturbed gas has a line width $\sigma_{out} = 515 \pm 3 \text{ km s}^{-1}$, which is ~4 times broader than the systemic component, and prominent blue wings extending to $v < -1000 \text{ km s}^{-1}$ (light-yellow region in Fig. 2). The inferred outflow velocity for the [O III] line is consistent with the general $L_{AGN}-v_{out}$ trends reported in the literature (Villar Martín et al. 2020; Musiimenta et al. 2023). Specifically, it appears to be on the higher end of the distribution, suggesting slightly faster outflows compared to other sources with similar luminosities. Therefore, our results support the recent findings by Tozzi et al. (2024) showing that highly obscured and CT AGN tend to exhibit faster outflows than unobscured AGN at a given $L_{AGN}$ (see also Bertola et al., in prep.). A similar conclusion could be reached for another highly obscured QSO in our GA-NIFS survey, ALESS073.1, presented in Parlanti et al. (2024b). The GS133 [O III] outflow is investigated in more detail in Sect. 4.

## 3.2. UV lines in VIMOS spectrum

UV lines were fitted with a multi-step approach, which we summarise here. In the next section we describe our approach in full detail. First, we shifted the wavelength axis to the rest frame, using the redshift inferred from the NIRSpec data. Then, we fitted the continuum emission at >1220 Å with a power-law,
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** VLT/VIMOS spectrum with best-fit results. The flux is normalised to the continuum. The blue curve identifies the integrated spectrum; the total, multi-component best fit curve is in red, while all individual Gaussian (Voigt) components are shown in yellow (green). The most prominent line transitions are marked with grey vertical lines.</figcaption></figure>

by masking all prominent line features (both in absorption and emission). Next, we normalized the continuum to 1 by dividing the original spectrum by the best-fit power-law profile. We subsequently fitted the emission lines with multi-component Gaussian profiles, and the absorption contribution with Voigt profiles. Finally, we obtained a synthetic spectrum without emission line contribution, by subtracting the best-fit emission line Gaussian profiles. This synthetic spectrum, only containing the contribution from absorption features, was finally ingested to VoigtFit (Krogager 2018), to obtain the main physical properties of the various atomic species, namely the column density and kinematics.

Because the VANDELS spectra are calibrated in air and not in vacuum, we used air wavelengths to fit the position of the different lines. The next sections describe our approach in detail.

### 3.2.1. Multi-component fit of emission and absorption lines

To model the low- and high-ionisation lines in the VIMOS spectrum, we used independent fits. This approach accounts for the fact that these species can trace different gas phases, potentially resulting in different systemic velocities and line widths (e.g. Shen et al. 2011; Lanzuisi et al. 2015, 2024). We classified the transitions based on their ionisation potential (IP). High-ionisation lines (HILs) were those with an IP between 30 and 80 eV, and include: N v λλ1238.82, 1242.80, Si iv λλ1393.76, 1402.77, [N IV] λ1483.32, N IV] λ1486.50, C IV λλ1548.19, 1550.77, and He ii λ1640.42. Low-ionisation lines (LILs), with IP < 30 eV, included the following transitions: C II λ1334.53, Si II λ1260.42, Si III λλ1298.9, 1303.32, the fluorescent doublets C II* λλ1335.66, 1335.71 and Si II* λλ1264.74, 1265.02, and the doublet C III] λλ1906.68, 1908.73. We used single Gaussian profiles to fit the emission lines, and Voigt profiles to model the absorption line features of both HILs (N v, C IV, and Si IV) and LILs (C II, C II*, Si II*, Si III). The fitting process used the Levenberg–Markwardt least-squares fitting code CAP-MPFIT (Cappellari 2017).

To achieve better constraints on the properties of transitions with lower S/N and to reduce fit degeneracy, especially in carbon and silicon fine structure lines with potential emission and absorption contributions, we simultaneously fitted transitions within the two individual groups of HILs and LILs (e.g. Perna et al. 2015). Namely, we constrained the wavelength

separation between line transitions in accordance with atomic physics (considering their air wavelengths); moreover, we fixed their widths (in km s⁻¹) to be the same for all emission lines. For the emission line doublets C IV, N v, and Si iv, the flux ratios were allowed to vary between the optically thick (f(B)/f(R) = 1) and thin (f(B)/f(R) = 2) limits, where R and B represent the red and blue transitions of each specific doublet (e.g. Del Zanna et al. 2002). The line ratio C III] f(1907)/f(1909), sensitive to the electron density, was allowed to vary within the physical variation range 0–1.6 (e.g. Keenan et al. 1992). The HIL features show multiple peaks in absorption; therefore, we considered up to two kinematic components to model them. For the LILs in absorption, only one kinematic component was considered. During the fitting procedure, the flux ratio between the absorption lines of each species was treated as free parameter.

Figure 3 shows the best-fit results for both HILs and LILs. The emission lines have a consistent systemic velocity (Δv$_{LIL}$ = −3+27−6 km s⁻¹; Δv$_{HIL}$ = 35+5−40 km s⁻¹). However, the emission line profiles are broader in HILs (σ$_{HIL}$ = 630+40−5 km s⁻¹) than in LILs (σ$_{LIL}$ = 330+60−5 km s⁻¹). Similarly, HILs in absorption are broader and require two kinematic components: the first one is blueshifted by ~−900 km s⁻¹, the second by ~−2000 km s⁻¹; both have a Doppler parameter b ~ 400 km s⁻¹. Conversely, LILs in absorption only require a Voigt kinematic component, blueshifted by ~−900 km s⁻¹ and with b ~ 400 km s⁻¹; therefore, this kinematic component perfectly match the low-velocity one identified in the HILs in absorption. These properties allow us to classify GS133 as a mini-BAL AGN.

As mentioned above, the degeneracy between emission and absorption contributions may affect our results. In particular, N v, Si ii, Si iv, C IV (and C II) show complex P Cygni profiles and are therefore associated with two (three) line transitions with potential contributions both in absorption and emission. Our fitting approach allowed us to reduce this degeneracy in the specific species thanks to the simultaneous modelling of other transitions like Si III, detected only in absorption, and He II and C III], detected only in emission. Nevertheless, a notable limitation is that our model does not identify any absorbing kinematic component at the systemic redshift, although such a component might exist. The combined impact of fit degeneracies and the low spectral resolution of the VIMOS data significantly constrains our ability to accurately model these complex line profiles. These limitations are further explored in the following sections.
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** VoigtFit models in velocity space. The normalised spectra, shown in blue with $3\sigma$ error bars, display the regions in the vicinity of the absorption lines; the best-fit models are overlaid in red. Carbon, nitrogen and silicon lines are fitted simultaneously, with two kinematic components; the vertical lines identify the velocity centroids for each kinematic component. For the doublets, the reported velocities are calculated with respect to the systemic of the blue line of the doublet; the kinematic components associated with the blue transitions are in orange, those of the red transitions are in purple. The narrow top panels show the best-fit residuals, with $3\sigma$ confidence intervals. Emission line contributions were removed prior to the VoigtFit modelling.</figcaption></figure>

### 3.2.2. Modelling of absorption lines with VoigtFit

After removing the emission line contributions inferred from the previous step, we modelled the most prominent absorption lines in the VIMOS spectrum with the Python package VoigtFit (Krogager 2018). Namely, we considered C II $\lambda$1334.53, and the N V $\lambda\lambda$1238.82, 1242.80, Si IV $\lambda\lambda$1393.76, 1402.77, and C IV $\lambda\lambda$1548.19, 1550.77 doublets.

We fitted the transitions from various ions and elements simultaneously, without separating the HILs and LILs, considering two distinct kinematic components. This is justified by the fact that independent fit results presented in Sect. 3.2.1 already showed that the high and low ionisation transitions have similar profiles at the lowest velocities. Figure 4 shows the best-fit results for the four species. The LIL C II is blue-shifted by $-778 \pm 26$ km s$^{-1}$, and is symmetric, with $b = 360 \pm 45$ km s$^{-1}$. The same blue-shifted kinematic component is also detected in the three HILs; these lines, however, also required a kinematic component with more extreme properties, namely $\Delta v = -1910 \pm 50$ km s$^{-1}$ and $b = 490 \pm 80$ km s$^{-1}$. VoigtFit also provides a measurement of the column density for each transition and kinematic components; all best-fit parameters are reported in Table 2.

<table><caption>**Table 2.** VoigtFit best-fit parameters and inferred column densities.</caption><thead><tr><th></th><th>log(N/cm<sup>-2</sup>)</th><th>log(N<sub>H</sub>/cm<sup>-2</sup>)</th><th>log(N<sub>H</sub><sup>c</sup>/cm<sup>-2</sup>)</th></tr><tr><th></th><th>(1)</th><th>(2)</th><th>(3)</th></tr></thead><tbody><tr><td colspan="4">Δv = -778 ± 26 km s<sup>-1</sup>, b = 340 ± 50 km s<sup>-1</sup></td></tr><tr><td>C II 1335</td><td>15.23 ± 0.07</td><td>>18.77</td><td>>18.9-19.0</td></tr><tr><td>C IV 1548, 1551</td><td>15.22 ± 0.05</td><td>>18.76</td><td>>18.8-18.9</td></tr><tr><td>N V 1239, 1243</td><td>15.32 ± 0.07</td><td>>19.42</td><td>>19.53</td></tr><tr><td>Si IV 1394, 1403</td><td>14.88 ± 0.05</td><td>>19.27</td><td>>19.5-20.6</td></tr><tr><td colspan="4">Δv = -1910 ± 50 km s<sup>-1</sup>, b = 410 ± 80 km s<sup>-1</sup></td></tr><tr><td>C II 1335</td><td>14.00 ± 0.70</td><td>-</td><td>-</td></tr><tr><td>C IV 1548, 1551</td><td>14.74 ± 0.07</td><td>>18.28</td><td>>18.4-18.5</td></tr><tr><td>N V 1239, 1243</td><td>14.99 ± 0.09</td><td>>19.09</td><td>>19.2</td></tr><tr><td>Si IV 1394, 1403</td><td>14.60 ± 0.06</td><td>>18.99</td><td>>19.2-20.3</td></tr></tbody></table>

**Notes.** Col. (1): Column density for individual species. Col. (2): Inferred hydrogen column density lower limit from individual species. Col. (3): Inferred hydrogen column density lower limit corrected for dust depletion. C II is well fitted with a single kinematic component (with $\Delta v \sim -800$ km s$^{-1}$); however, the table also reports the tentative column density obtained by VoigtFit for the second kinematic component.

To estimate a conservative lower limit to hydrogen column density $N_H$ in the outflow, we assumed no ionisation correction (e.g. N(C) = N(C IV)), and no correction for saturation. We also assumed solar abundances and the possible effects of depletion onto dust from Jenkins (2009), following Bordoloi et al. (2014). In Table 2 we report for each transition a value without correction for dust depletion (third column) and a range of values for depletion corrected densities (last column) obtained considering the range of values in Jenkins (2009, see their Table 4).

Overall, the column density lower limits obtained from the different species are within a narrow range of $log(N_H/cm^{-2}) = 18.5-19.5$, for both kinematic components. The main discrepancies are observed in the depletion corrected values for the Si IV transitions, which reach values up to $log(N_H/cm^{-2}) = 20.6$: the silicon element, however, has the largest and likely more uncertain depletion corrections, in the range [-0.22, -1.36]. Because of this, we considered in the derivation of the outflow energetics (Sect. 7) an order of magnitude approximation $log(N_H/cm^{-2}) = 19$. This column density is consistent with values generally associated with mini-BALs and NALs (e.g. Laha et al. 2021).

### 3.3. Possible Balmer absorption in the NIRSpec spectrum

The spatially integrated optical spectrum of GS133 in Fig. 2 displays two absorption features on the left side of the $H\alpha$-[N II] complex. Figure 5 shows that they are at $\approx -1000$ and $\approx -2000$ km s$^{-1}$ from the $H\alpha$ peak, hence broadly matching the mini-BAL features observed in the VIMOS spectrum (see Fig. 4). Therefore, these two features could be interpreted as Balmer absorption, similar to those observed in a few broad-line AGN up to $z \sim 2$ in the pre-JWST era (e.g. Hutchings et al. 2002; Schulze et al. 2018; Hamann et al. 2019) and many JWST-selected Type 1 AGN (e.g. Juodžbalis et al. 2024; Kocevski et al. 2024; Matthee et al. 2024; Wang et al. 2024).

These Balmer absorption lines are thought to originate around the edge of the obscuring torus, which is eroded and accelerated by nuclear winds; in fact, their inferred distances from the central engine are larger than the size of the BLR but smaller than the size of the dusty torus based on detailed
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** GS133 spectrum showing tentative Balmer absorption detection in the vicinity of Hα and [N II] emission lines. The black curve shows the integrated spectrum in velocity space, with $v = 0$ km s⁻¹ at the position of the Hα line. The red curve (with $1\sigma$ uncertainties in light-red) displays the residuals obtained by subtracting the best-fit Gaussian emission lines presented in Fig. 2. Two low S/N absorption features are detected close to the orange vertical lines, which mark the velocities of the mini-BAL features observed in the UV (see Table 2).</figcaption></figure>

photoionisation models using CLOUDY (e.g. Zhang et al. 2015; Shi et al. 2016; Juodžbalis et al. 2024). Alternatively, Balmer absorption could be associated with larger distances from the SMBH if the absorption is weak and the hydrogen density and $N_H$ are low (see Sect. 7.1).

However, it is important to note that the absorption features detected in the integrated spectrum of GS133 have relatively low S/N (~2). The subtraction of emission line contributions resulted in negative flux values in some velocity channels, indicating that the absorption feature shapes may be influenced by fit degeneracies. Consequently, we did not pursue further analysis of these tentative Balmer absorption detections.

# 4. Spatially resolved analysis

This section presents the spatially resolved properties of the ionised gas in GS133 as inferred from the NIRSpec data-cube.

## 4.1. Velocity channels and spectro-astrometry

Before performing a spaxel-by-spaxel spectral fit analysis, we generated [O III] velocity channel maps and used spectro-astrometry analysis to understand if the ionised gas is spatially resolved in our NIRSpec data.

Figure 6 shows the [O III] maps obtained by integrating over three velocity ranges: [−800, −200] km s⁻¹ to cover the high-velocity blueshifted (approaching) component, [−200, 200] km s⁻¹ to cover the systemic emission, and [200, 800] km s⁻¹ for the high-velocity redshifted (receding) gas. These maps show an extended morphology covering ~1'' (7.5 kpc) along the NE to SW direction, in particular in the systemic and redshifted components; the blueshifted emission is more compact (contours identify $3\sigma$ and $5\sigma$ levels, in each panel). The NIRSpec PSF (reported in the bottom-right part of the third map) is elongated along the same NE-SW direction, roughly corresponding to the direction of the NIRSpec IFS slices (~60°), but on smaller scales: FWHM ~0.13'' along the slices and ~0.1'' across the slices (D'Eugenio et al. 2024). Therefore, the [O III] emission in GS133 is spatially resolved.

To better understand whether the high-velocity components are spatially separated from the systemic gas in the host galaxy, we also performed a spectro-astrometry analysis. This analysis consists in determining how the centroid position of the [O III] emission changes as a function of velocity. We followed a methodology similar to that of Pereira-Santaella et al. (2018) and Lamperti et al. (2022). We binned together the velocity channels in intervals of 120 km s⁻¹ to increase the S/N, necessary to reliably determine the position of the peak of the emission. Then, we performed a fit with the Photutils package<fnref n="2" /> to determine the peak position in each binned channel. In this way, we could determine the peak positions at sub-pixel scales.

The [O III] spectro-astrometry positions are reported in the right panel of Fig. 6. The low-velocity ($|v| \le 120$ km s⁻¹) channels appear oriented along the NW to SE direction, while higher velocity ($|v| > 120$ km s⁻¹) are roughly perpendicular. This configuration suggests that GS133 hosts a rotating disc with a major axis perpendicular to the extended outflow (see also e.g. Fig. 3 in Lamperti et al. 2022 for a similar configuration in a nearby galaxy). The orientation of the [O III] morphological major axis, which aligns with the outflow rather than the rotational axis, suggests that outflow emission is dominant over the unperturbed gas and that a significant portion of the ejected [O III] may has relatively low projected velocities. This kinematic configuration is further investigated in the next sections.

## 4.2. Spaxel-by-spaxel multi-Gaussian fit

To derive spatially resolved kinematic and physical properties of ionised gas, we fitted the spectra of individual spaxels using the prescriptions already presented in Sect. 3.1. We applied the BIC selection to determine where a multiple-Gaussian fit is required to statistically improve the best-fit model. For each Gaussian component included in the fit, a S/N > 3 was required. This approach ensures that the more complex, and potentially degenerate, multiple-component fits are used only where statistically justified. In most spaxels, two Gaussian components suffice, but a few spaxels closer to the AGN position require up to four Gaussian components.

Figure 7 shows the $\Delta v$–FWHM velocity diagram obtained from our multi-Gaussian models, with the kinematic parameters of all Gaussian components required to fit the datacube. The measurements are coloured by the distance from the nucleus. Although the diagram does not show a clear trend, we note that the highest FWHMs (>400 km s⁻¹) are sometimes associated with significant velocity offsets ($|\Delta v| > 150$ km s⁻¹), as observed in systems hosting AGN outflows (e.g. Woo et al. 2016; Perna et al. 2022, 2023a). In particular, the most external regions (yellow points in the figure) show a positive trend, where both velocity and FWHM increase together; the broadest components are instead blueshifted and closer to the nucleus. In contrast, the narrower Gaussian components have relatively small offsets from the zero velocity.

To check for the presence of rotationally supported motions (see the right panel in Fig. 6), we selected in the $\Delta v$–FWHM plane the components with $|\Delta v| < 120$ km s⁻¹ and FWHM < 400 km s⁻¹ (i.e. within the blue box in Fig. 7), and constructed new datacubes containing only the systemic Hα and [O III] emission (see e.g. Perna et al. 2022 for a similar approach). The flux distribution as well as the velocity offset and dispersion of the narrow Hα and [O III] lines are reported in Fig. 8. The $\Delta v$ maps display a regular pattern along the NW-SE

<footnote n="2">² https://photutils.readthedocs.io/en/stable/</footnote>
<figure><img src="image_7.png" /></figure>

**Fig. 6.** [O III] velocity-channel maps (first-to-third panels) and spectroastrometry [O III] line positions (right panel). First to third panels: The velocity-channel maps were extracted from the ranges labelled in the individual panels. For each flux distribution, we report 3 and $5\sigma$ contour levels. The black cross identifies the GS133 nucleus, as inferred from spectroastrometry position of spectrally integrated [O III] profile (i.e. considering all velocity channels together); the white ellipse in the third panel represents the NIRSpec point spread function (PSF) at ~2.2 µm (the [O III] observed wavelength) estimated by D'Eugenio et al. (2024); a scale bar and a compass are also reported in the first map. Right panel: Spectroastrometry of the [O III] emission, i.e. centroid positions of the line in the different velocity channels, covering a relatively small portion of the spatial extent shown in the velocity-channel maps (see white box in the third panel). The points are colour-coded by the channel velocity. The dot-dashed line roughly identifies the outflow axis, while the dotted line shows the tentative kinematic major axis of the GS133 rotating disc.

<figure><img src="image_8.png" /></figure>

**Fig. 7.** Velocity diagram for the individual Gaussian components used to model the emission line profiles in the GS133 datacube. The measurements are coloured by the distance from the AGN. The blue lines isolate the Gaussian components with FWHM < 400 km s⁻¹ and |Δv| < 120 km s⁻¹ used to separate the narrow emission line kinematics (reported in Fig. 8) from the outflow kinematics (Fig. 9).

direction, consistent with the spectro-astrometry analysis results (Sect. 4.1), with a peak-to-peak velocity shift of ~80 km s⁻¹. The velocity dispersion maps do not display a clear peak at the position of the nucleus, as expected for a rotating disc (e.g. Förster Schreiber et al. 2018). This could be due to the presence of the outflow and fit degeneracies. Indeed, the Δv map shows some redshifted emission towards the NE and blueshifted towards the SW (mostly in [O III]), which do not follow the rotation pattern seen in the more central regions; this component could be similarly due to a portion of the ejected gas with relatively small projected velocities. This scenario is further discussed in Sect. 7.4.

Figure 9 displays the flux, velocity offset, and velocity dispersion obtained for the outflow Gaussian components (outside the box identified in Fig. 7). The [O III] gas is more extended than $H\alpha$, and shows a bi-conical geometry, typical of AGN outflows

(e.g. Venturi et al. 2018; Perna et al. 2022; Falcone et al. 2024). The ionised gas is blueshifted (Δv ~ −280 km s⁻¹) towards the NE, and redshifted (Δv ~ 100 km s⁻¹) towards the SW, up to ~3 kpc. The line widths reach values as high as ~450 km s⁻¹ in the circum-nuclear regions, in line with the integrated fit results (Sect. 3.1).

For completeness, in Fig. A.1 we report the flux and velocity maps obtained from the total, integrated emission lines, hence without separating the systemic from the outflow components. In these maps, the signatures of the rotating disc and the bi-conical outflow configuration are less pronounced, underlining the complex kinematics in GS133.

# 5. Multi-wavelength diagnostics

## 5.1. AGN classification with optical line ratio diagnostics

We investigated the dominant ionisation source for the emitting gas across the GS133 host using the classical 'Baldwin, Phillips & Terlevich' (BPT) diagram (Baldwin et al. 1981). Figure 10 shows the distribution of the flux ratio diagnostics across the GS133 host galaxy extension (green-to-blue squares for increasing x- and y-axis values of the BPT diagram). The figure also shows a few spatially integrated flux ratios: the green and purple large circles refer to the systemic and outflow components obtained from the spectrum in Fig. 2; the light-blue and light-red crescent moon markers display the flux ratios of the outflow component in the integrated spectra extracted from the approaching and receding parts of the biconical outflow, respectively (using circular apertures with $r = 0.15''$, see inset in Fig. 10), and reported in Figs. B.1 and B.2. For reference, the BPT also displays other optical line ratio measurements from the literature, for low-z SDSS galaxies (small grey points), and star-forming galaxies and AGN at $z > 2.6$ recently observed by JWST (from Scholtz et al. 2023; Perna et al. 2023b; Calabrò et al. 2023), as well as the demarcation lines used to separate galaxies and AGN at $z \sim 0$ from Kewley et al. (2001) and Kauffmann et al. (2003).

The spatially resolved line ratio measurements do not show a significant variation across the GS133 host galaxy extension, although a tentative trend can be observed. The lowest [O III]/$H\beta$ and [N II]/$H\alpha$ line ratios are observed in the outskirts of the
Narrow components

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** [O III] (top) and Hα (bottom) flux, moment-1 and moment-2 maps for the systemic kinematic component. A S/N cut of 4 has been applied to generate the maps. The moment-1 maps of both lines show evidence of rotating gas in the QSO host, but neither moment-2 map shows a peak at the nuclear position (identified by a black cross).</figcaption></figure>

Outflow components

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** [O III] (top) and Hα (bottom) flux, moment-1 and moment-2 maps for the outflow kinematic component. A S/N cut of 4 has been applied to generate the maps. The [O III] maps show a typical biconical outflow structure, oriented along the NE-SW direction, consistent with our spectro-astronomy analysis (Sect. 4.1). The Hα flux distribution is less extended, as this Balmer line is fainter than [O III].</figcaption></figure>

galaxy, towards NW and SE (see the map in the inset), and are very close to the BPT locus of star-forming galaxies at high-z (green symbols in the figure). Instead, the highest line ratio measurements are along the outflow axis, likely corresponding to the AGN ionisation cones of GS133. In fact, these high line ratios are in a BPT region free from contamination by star-forming galaxies at any redshift, and only populated by AGN systems (e.g. Scholtz et al. 2023; Perna et al. 2023a,b; Parlanti et al. 2024b; Decarli et al. 2024).

The spatially integrated measurements reported in the figure (large circles and crescent moon symbols) similarly indicate that the systemic narrow emission is consistent with star-formation ionisation, while the outflow components are AGN dominated (see also line ratios reported in Table 3).

In conclusion, both spatially integrated and spatially resolved line ratios support an AGN classification based on the standard BPT diagnostic diagram. In fact, while the BPT diagram may lose sensitivity in distinguishing between AGN and star-forming galaxies in low-mass, low-metallicity systems at high redshift, it remains effective for identifying AGN in more massive and metal-enriched galaxies (e.g. Z > 0.5 Z⊙, Feltre et al. 2016).

## 5.2. AGN classification with UV line ratio diagnostics

The direct observation of high-ionisation emission lines in the UV spectrum of GS133, such as Nvλλ1239, 43, C ivλλ1548, 51, He iiλ1640, also in combination with other UV transitions, can be used to obtain a further confirmation of the presence of AGN emission in this high-z source (e.g. Mascia et al. 2023; Scholtz et al. 2023; Maiolino et al. 2024b).

In Table 4 we report the main UV line ratios commonly used to distinguish between AGN and star-forming galaxies in the literature (e.g. Feltre et al. 2016; Nakajima et al. 2018; Perna et al. 2023b; Topping et al. 2024). All of these line ratios are compatible with AGN ionisation, which is consistent with the X-ray and optical BPT classifications introduced in the previous sections.

## 5.3. CT AGN classification with L([O III])/L(X-ray) diagnostic

X-ray diagnostics are the most reliable methods for identifying CT AGN (e.g. Lanzuisi et al. 2017). Previous analyses by Luo et al. (2017) and Li et al. (2019) of the Chandra X-ray spectrum for GS133 suggested column densities of $N_H \approx 1.5 \times$
<figure><img src="image_11.png" /><figcaption>Fig. 10. Standard BPT diagnostic diagram. The colour-coded squares show GS133 single-spaxel measurements associated with the spatial regions shown in the top-right panel; green-to-blue colours mark increasing line ratios, as indicated with the arrow in the bottom-left part of the diagram; only spaxels with S/N > 3 in all lines are shown. The large green and purple circles display the line ratios measured from r = 0.5" spatially integrated regions, for the systemic and outflow components, as labelled: the large crescent moon symbols refer to the outflow component line ratios measured in circular regions marked in the inset (with the same crescent moon symbols). The solid (Kewley et al. 2001) and dashed (Kauffmann et al. 2003) curves commonly used to separate purely star-forming galaxies (below the curves) from AGN (above the curves) are also reported. Finally, the BPT shows local galaxies from SDSS (Abazajian et al. 2009, indicated in grey), and additional measurements for z = 2.6–6.7 sources from the literature: star-forming galaxies are marked with small green symbols and refer to CEERS (Calabrò et al. 2023) and JADES (Scholtz et al. 2023) surveys, while AGN are marked with small orange symbols and refer to JADES (Scholtz et al. 2023) and GA-NIFS (Perna et al. 2023b) surveys.</figcaption></figure>

10$^{24}$ cm$^{-2}$, which is around the threshold typically used to classify AGN as CT. However, due to the significant uncertainties associated with the low number of X-ray counts typical for high-redshift sources like GS133, where only 70.7 aperture-corrected source counts are available in the 0.5–7 keV band (Luo et al. 2017), the true column density could be lower. To strengthen the case, multi-wavelength diagnostics have been employed by Guo et al. (2021). They compared the X-ray to mid-infrared (6 µm) luminosity, commonly used to identify CT sources as the infrared emission is less affected by obscuration than X-rays; they found obscuration exceeding 1.5×10$^{24}$ cm$^{-2}$. In this section, we tested an additional diagnostic to further support the CT classification of the AGN in GS133.

Maiolino et al. (1998) provided a diagnostic based on the ratio between the observed 2–10 keV and the reddening-corrected [O III] luminosities to infer the column density of obscured AGN. Adopting their methodology, we derived a luminosity ratio log(L$_X$/L$_{[OIII]}$) = −1.4, corresponding to a lower limit log(N$_H$/cm$^{-2}$) > 2 × 10$^{24}$ (see their Fig. 6). Specifically, the X-ray to [O III] luminosity ratio was calculated using the [O III] outflow luminosity reported in Table 1 (hence excluding any possible contribution from star formation). The [O III] luminosity was corrected for dust extinction considering an intrinsic Balmer decrement of 2.86 and assuming the Cardelli et al. (1989) extinction law. Therefore, this result further supports the CT nature of GS133.

## 6. Dynamical mass

Having an estimation of the kinematics and extent of a galaxy allows us to provide a constraint on the dynamical mass of the galaxy, M$_{dyn}$. Under the assumption that the gas is distributed in a flat disc, the dynamical mass enclosed within a radius R is M$_{dyn}$ = 2.33 × 10$^{5}$ v$_{circ}$ R M$_{⊙}$, where v$_{circ}$ is the circular velocity in km s$^{-1}$ at a galactocentric distance R, given in kpc (Walter et al. 2004; Neeleman et al. 2021; Perna et al. 2022). Since we could only marginally resolve the disc structure in GS133, we assumed an extent of R = 1 kpc, typical of z = 3–4 galaxies (Kartaltepe et al. 2023; Costantin et al. 2024), and that the FWHM (i.e. 2.355 × σ) of the systemic Hα line can be used as a proxy for the circular velocity, so that v$_{circ}$ = 0.75 FWHM/sin(i) = 300 km s$^{-1}$, assuming a mid-range inclination of 55° following, for instance, Neeleman et al. (2021) and Decarli et al. (2018).

Under these assumptions, the dynamical mass of GS133 is M$_{dyn}$ = 2 × 10$^{10}$ M$_{⊙}$, broadly consistent with the stellar mass inferred by Guo et al. (2020) and Circosta et al. (in prep.) from spectral-energy-distribution (SED) analysis: M$_{*}$ = 5 × 10$^{9}$ M$_{⊙}$, and 3 × 10$^{10}$ M$_{⊙}$, respectively. The difference in stellar mass estimates can be explained by the common degeneracy between AGN and stellar components in the optical and UV regime (e.g. Ciesla et al. 2015). A more comprehensive comparison of the SED fitting results will be provided in Circosta et al. (in prep.).

## 7. Outflow properties

In this section we measure the mass rate and energetics of the ionised outflow as inferred from the blueshifted outflow components of absorbing (C II, C IV, and N V) and emitting (Hβ and [O III]) line transitions. We begin with standard assumptions commonly used in the literature, then refine these measurements by applying detailed photoionisation modelling for the absorption lines and 3D kinematic modelling for the emission lines. We also compare our results with outflow model predictions and general AGN properties to exclude less plausible assumptions, aiming to derive more accurate outflow properties for both the absorbing and emitting gas.

### 7.1. Energetics of absorbing gas outflow

We measured the mass rate of the absorbing outflow gas as inferred from the blueshifted mini-BAL components of C II, C IV, and N V. Following Bordoloi et al. (2014), we assumed for simplicity a thin-shell approximation (Ω = 4π), and maximised to unity the angular and clumpiness covering factors ($C_{\Omega}C_{f} = 1$):

$$ \dot{M}(\text{Abs}) \approx \frac{N_{\text{H}}}{10^{20} \text{cm}^{-2}} \frac{R_{\text{out}}}{1 \text{kpc}} \frac{v_{\text{out}}}{100 \text{km s}^{-1}} M_{\odot} \text{yr}^{-1}, \quad (1) $$

where N$_{H}$ is the hydrogen column density and R$_{out}$ and v$_{out}$ are the mini-BAL galactocentric radius and velocity, respectively. The kinetic power of the absorbing material can be derived from the relation $\dot{E} = 1/2 \dot{M}v^2$:

$$ \dot{E}(\text{Abs}) \approx 5 \times 10^{39} \frac{N_{\text{H}}}{10^{20} \text{cm}^{-2}} \frac{R_{\text{out}}}{1 \text{kpc}} \left( \frac{v_{\text{out}}}{100 \text{km s}^{-1}} \right)^3 \text{erg s}^{-1} \quad (2) $$

(see also Hamann et al. 2019). We assumed that the absorbing outflow is related to the two distinct components travelling at different velocities (Fig. 4), and that v$_{out}$ is approximated by the velocity shift Δv of a given component, hence ~−800 km s$^{-1}$ and
<table><caption>**Table 3.** Optical emission line ratios from NIRSpec integrated spectra.</caption><thead><tr><th rowspan="2">Line ratio</th><th colspan="3">Integ. spec. (r = 0.5")</th><th colspan="3">NE cone (r = 0.15")</th><th colspan="3">SW cone (r = 0.15")</th></tr><tr><th>tot</th><th>sys</th><th>out</th><th>tot</th><th>sys</th><th>out</th><th>tot</th><th>sys</th><th>out</th></tr></thead><tbody><tr><td>log([O III]/Hβ)</td><td>0.66<sup>+0.01</sup><sub>-0.02</sub></td><td>0.65<sup>+0.04</sup><sub>-0.01</sub></td><td>0.96<sup>+0.07</sup><sub>-0.10</sub></td><td>0.96<sup>+0.03</sup><sub>-0.02</sub></td><td>1.07<sup>+0.10</sup><sub>-0.03</sub></td><td>0.82<sup>+0.09</sup><sub>-0.03</sub></td><td>0.84 ± 0.1</td><td>1.01<sup>+0.09</sup><sub>-0.03</sub></td><td>0.75<sup>+0.01</sup><sub>-0.04</sub></td></tr><tr><td>log([N II]/Hα)</td><td>-0.57<sup>+0.16</sup><sub>-0.03</sub></td><td>-0.57<sup>+0.06</sup><sub>-0.01</sub></td><td>-0.11<sup>+0.35</sup><sub>-0.21</sub></td><td>-0.54<sup>+0.03</sup><sub>-0.04</sub></td><td>-0.45<sup>+0.08</sup><sub>-0.17</sub></td><td>-0.40 ± 0.1</td><td>-0.51 ± 0.01</td><td>-0.84<sup>+0.24</sup><sub>-0.05</sub></td><td>-0.32<sup>+0.04</sup><sub>-0.02</sub></td></tr></tbody></table>

**Notes.** Integrated spectrum emission lines from best-fit shown in Fig. 2; NE and SW cones from best-fit shown in Figs. B.1 and B.2, respectively.

<table><caption>**Table 4.** UV emission line ratios from the VIMOS spectrum.</caption><thead><tr><th>CIII]<br>HeII</th><th>CIV<br>HeII</th><th>CIV<br>CIII]</th><th>CIII]+CIV<br>HeII</th><th>NV<br>HeII</th><th>CIV<br>NV</th><th>OIII]<br>HeII</th><th>log(L(He II)/[erg s<sup>-1</sup>])</th></tr></thead><tbody><tr><td>-0.32<sup>+0.18</sup><sub>-0.02</sub></td><td>0.35 ± 0.01</td><td>0.68<sup>+0.10</sup><sub>-0.06</sub></td><td>0.44<sup>+0.04</sup><sub>-0.01</sub></td><td>0.14<sup>+0.01</sup><sub>-0.01</sub></td><td>0.21<sup>+0.02</sup><sub>-0.01</sub></td><td>&lt;-0.96</td><td>42.2 ± 0.1</td></tr></tbody></table>

**Notes.** O III]λλ1660.81, 1666.15 is not detected in the VIMOS spectrum.

~1900 km s⁻¹ (Table 2). For the extent of the mini-BAL outflow, we initially used a range of values between 1 pc and 10 kpc (following Dunn et al. 2010; Saturni et al. 2016; Moravec et al. 2017; Bruni et al. 2019), as it cannot be inferred without knowing the hydrogen density in the absorbing ejected gas. For the column density estimate, we considered median values inferred from dust-depletion corrected N<sub>H</sub><sup>c</sup> reported in Table 2, excluding the Si IV estimates as they are associated with more uncertain corrections (see Jenkins 2009): 8 × 10<sup>18</sup> cm<sup>-2</sup> and 1 × 10<sup>19</sup> cm<sup>-2</sup> for the low- and high-velocity outflows, respectively. We obtained M̄(Abs) = 0.003–30 M<sub>☉</sub> yr<sup>-1</sup> and Ē(Abs) = (0.0004–4) × 10<sup>43</sup> erg s<sup>-1</sup>, considering the extremely large range of possible values for R<sub>out</sub>, and computing the sum of the contributions from two kinematic components (with ~70% coming from the faster and denser component). In the next subsection, these measurements are refined taking into account detailed photoionisation models.

## 7.2. Photoionisation modelling for ejected absorbing gas

To infer the physical conditions within the mini-BAL, we computed a series of simulations using the photoionisation code CLOUDY (c17.03; Ferland et al. 2017) and compared them with observations. We set up a series of 1D single-cloud models with the ionising SED set to a CLOUDY built-in multi-component AGN spectrum described by the following formula

$$F_v = v^{\alpha_{UV}} \exp(-hv/kT_{BB}) \exp(0.01 \text{Ryd}/hv) + av^{\alpha_{X}}. \quad (3)$$

Specifically, we set the cut-off temperature for the accretion disc to T<sub>BB</sub> = 10<sup>6</sup> K, the optical-to-X-ray ratio to α<sub>ox</sub> = -1.4, the X-ray slope to α<sub>X</sub> = -1, and the UV slope to α<sub>uv</sub> = -0.5 (by adjusting parameter a in the equation above). The SED we adopted roughly reproduces the shape of the observed UV continuum observed by VLT/VIMOS. For the absorbing clouds, we explored a wide range of ionisation parameters<sup>3</sup>, U, from 10<sup>-4.5</sup> to 1. We set the clouds to be isobaric and vary the density in a range of 1 ≤ log(n<sub>H</sub>/cm<sup>-3</sup>) ≤ 7. All models are computed until a cumulative hydrogen column density of N<sub>H</sub> = 10<sup>22</sup> cm<sup>-2</sup>, and cumulative column densities for different atomic/ionic transitions as functions of depth are saved individually for the search of solutions that match the derived values for GS133. For the chemical abundances, we adopted the solar abundance set of

Grevesse et al. (2010) and set the ratio N/C as a free parameter to account for variations due to chemical enrichment history (Maiolino & Mannucci 2019). Finally, since Si is usually highly depleted onto dust grains in the ISM (e.g. Jenkins 2009), we set Si/C as an additional free parameter to probe the level of dust depletion.

We considered separate models for the two kinematic components at ~800 km s⁻¹ and ~1900 km s⁻¹. We used a single-cloud model, where both C II and C IV can be found, but C IV absorption happens in the more inner region of the cloud due to its higher ionisation potential. This means that when C II has reached a large column density at a certain depth in the cloud, the column density of C IV usually has already reached its maximum value. The stopping conditions for the cumulative density at N<sub>H</sub> = 10<sup>22</sup> cm<sup>-2</sup> ensure that we match the values inferred from the VIMOS spectrum for both C II and C IV. Specifically, for the low-velocity component, we searched for models with column densities of the C II that match our derivations from the analysis of the VIMOS spectrum. Then, we allowed the column density of C IV to differ from the measured value by ±0.15 dex (i.e. 3σ, Table 2), and scaled the N/C and Si/C abundances to reproduce the column densities measured for N V and Si IV transitions. We also performed a minimum χ² search using the measured column densities for the above transitions and obtained consistent results.

In Fig. 11 we show the best-fit n<sub>H</sub> and N<sub>H</sub> for the absorbing gas. For the mini-BAL component with v<sub>out</sub> ~ -800 km s⁻¹, we found low-density models with log(N<sub>H</sub>/cm<sup>-2</sup>) = 20.2–20.3, log(n<sub>H</sub>/cm<sup>-3</sup>) = 1–2, and log(U) in the range [-2.7, -2.6] can well reproduce our measurements. High-density models with log(n<sub>H</sub>/cm<sup>-3</sup>) = 4–5 also have a good match with our measurements. All these models require a high log(N/C) = 1 to reproduce N V absorption, consistent with CNO equilibrium ratio in asymptotic giant branch (AGB) stars (Maeder et al. 2015). However, since there is no strong nebular emission of N IV] and N III] in the VIMOS spectrum, we caution the interpretation of the high N/C. As an alternative explanation, due to the high ionisation potential of N V, N V absorption might actually occur in a separate density-bounded cloud closer to the central AGN, which cannot be described with single-cloud models; fit degeneracy between emission and absorption contributions could also explain the high N/C. Based on the best-fit Si/C, we inferred a depletion factor of 0.52 dex for the Si IV. The depletion of Si IV is consistent with a moderate amount of dust depletion found in the compilation of MW sightlines by Jenkins (2009) with a unified depletion strength of F<sub>*</sub> = 0.25. In comparison,

<footnote n="3"><sup>3</sup> Throughout this work U is defined as the dimensionless ratio between the ionising photon flux and the product of hydrogen density and the speed of light.</footnote>
<figure><img src="image_12.png" /><figcaption>Fig. 11. Hydrogen densities and column densities predicted by CLOUDY photoionisation models for the two mini-BAL components. Only models with $\chi^2 < 5$ are shown and the sizes of the symbols are inversely scaled with $\chi^2$. In addition, models are colour coded according to their ionisation parameters. Both outflow components have low-density and high-density solutions, although the high-density solution would imply extremely thin layers of gas.</figcaption></figure>

a stronger depletion strength of $F_*$ = 0.5 is typical for star-forming galaxies at $z \lesssim 0.1$ (Gunasekera et al. 2023). For the high-velocity component, we found low-density models with $\log(N_H/cm^{-2}) = 19.4-19.6$, $\log(n_H/cm^{-3}) = 1.5-3.0$, and $\log(U)$ in the range [-2.8, -2.7] can reproduce our measurements if we leave N/C and Si/C as free parameters. There is a high-density solution at $\log(n_H/cm^{-3}) \approx 6$ as well. If we further constrain N/C and Si/C in the models to the values obtained for the slow component, the preferred models for the fast component would have $\log(N_H/cm^{-2}) \approx 19.6$, $\log(n_H/cm^{-3}) = 2.0-4.0$, and $\log(U)$ in the range [-2.8, -2.7].

For both low- and high-velocity components, our low-density models suggest that the effective thickness of the absorbing gas (i.e. $N_H/n_H$) can span a range from 4 to 0.04 pc if it has a volume filling factor close to unity, consistent with previous results for low-z systems from the literature (e.g. Choi et al. 2022). In comparison, high-density solutions correspond to physical scales of $6 \times 10^{-4}-10^{-3}$ pc for the low-velocity component and $10^{-5}$ pc for the high-velocity component.

Knowing the hydrogen density and ionisation parameter for the cloud from CLOUDY models, we could also infer the distance of the absorbing material from the SMBH, from

$$U = \frac{Q}{4\pi R_{\text{out}}^2 n_H c}, \qquad (4)$$

where Q is the number of ionising photons per unit time, proportional to the AGN bolometric luminosity (see e.g. Eq. 4 in Baron & Netzer 2019), and c is the speed of light. By solving this equation for $R_{\text{out}}$, we obtained distances of 1–10 kpc for the two kinematic components with low-density CLOUDY models reported above. We note that much smaller distances would require significantly higher densities of $n_H$ ($\approx 10^6$ cm$^{-3}$), which, however, would be associated with an unphysical thickness of the cloud, orders of magnitude smaller than typical dimensions of a single cloud in the BLR ($\lesssim 10^{-3}$ pc), and smaller than normally assumed in BAL models (e.g. r ~ 0.3–3 pc in Ishibashi et al. 2024). These considerations favour the low-density models as the most plausible.

Therefore, our CLOUDY models suggest that the mini-BAL material in GS133 is at the same distance of emitting gas and not confined in the nuclear (parsec-scale) regions. While

single-cloud models might be overly simplistic for a system like GS133, which exhibits multiple kinematic components in both low- and high-ionisation lines, as well as an extended, bi-conical [O III] outflow, more complex models with multiple clouds would introduce greater degeneracies (as they require the modelling of >10 line transitions; e.g. Marconi et al. 2024). Therefore, our approach represents a compromise, providing indicative properties of the mini-BALs in GS133.

Using the best-fit CLOUDY parameters, we recalculated the UV outflow energetics. For the low-velocity component, we found mass outflow rates of $14-140 M_{\odot} yr^{-1}$ and kinetic powers of $(5-50) \times 10^{42} erg s^{-1}$, with the range reflecting an outflow extent from 1 to 10 kpc. For the high-velocity component, the corresponding values are $6-60 M_{\odot} yr^{-1}$ and $(1-10) \times 10^{43} erg s^{-1}$, with the same range for the outflow extent.

### 7.3. Energetics of emitting gas outflow

We computed the mass outflow rate of the ionised gas as inferred from the blueshifted outflow component of $H\beta$, assuming a simplified biconical outflow distributed out to a radius $R_{\text{out}}$, following Cresci et al. (2015):

$$\dot{M}_{\text{out}}(Em, H\beta) = 8.6 \times \frac{L_{41}(H\beta) v_{\text{out}}}{n_e R_{\text{out}}} M_{\odot} yr^{-1}, \qquad (5)$$

where $L_{41}(H\beta)$ is the $H\beta$ luminosity associated with the outflow component in units of $10^{41} erg s^{-1}$ (corrected for dust extinction), $n_e$ is the electron density in cm$^{-3}$, and $v_{\text{out}}$ is the outflow velocity in km s$^{-1}$.

For comparison, we also computed the mass outflow rates employing the [O III] luminosity and assuming the same geometrical configuration:

$$\dot{M}_{\text{out}}(Em, [\text{OIII}]) = 0.5 \times \frac{L_{41}([\text{OIII}]) v_{\text{out}}}{n_e R_{\text{out}} 10^{[\text{O/H}]}} M_{\odot} yr^{-1}, \qquad (6)$$

where $L_{41}([\text{OIII}])$ is the [O III] luminosity associated with the outflow component in units of $10^{41} erg s^{-1}$ (corrected for dust extinction), and $10^{[\text{O/H}]}$ is the metallicity of the outflowing gas in solar units (see also Cano-Díaz et al. 2012). The kinetic power
of the ejected emitting gas can be derived as $\dot{E} = 1/2 \dot{M} v^2$:

$$ \dot{E}_{\text{out}}(Em) = 3 \times 10^{35} \left( \frac{\dot{M}(Em)}{M_{\odot} \text{yr}^{-1}} \right) \left( \frac{v_{\text{out}}}{\text{km s}^{-1}} \right)^2 \text{erg s}^{-1}. \quad (7) $$

For simplicity, the outflow energetics of the emitting gas were derived assuming a simple single-radius wind, because of the unknown detailed geometrical configuration of the outflow (i.e. opening angles and inclinations of the conical structures, deprojected sizes and velocities). The fact that the gas is ionised by the AGN did not allow us to derive metallicity measurements for the outflowing gas, but the high [N II]/H$\alpha$ line ratio (>0.3) and stellar mass of the system (~10$^{10}$ M$_{\odot}$) allows us to reasonably assume a solar metallicity to derive [O III]-based outflow energetics. We also assumed an outflow extension of 3 kpc, and an outflow velocity $v_{\text{out}} = 1000 \text{ km s}^{-1}$ for both [O III] and H$\beta$, as order of magnitude estimates (based on the observed distribution of high-$v$ [O III] gas in Fig. 6, and our multi-Gaussian fit decomposition results reported in Table 1). Similarly, we assumed an electron density $n_e = 1000 \text{ cm}^{-3}$, considering that we could not obtain a measurement for the ejected gas because of the faintness and complexity of the [S II] profiles, and the extremely high density derived from the narrow [S II] components (3300 ± 500 cm$^{-3}$) in the integrated spectrum shown in Fig. 2. With these assumptions, we obtained $\dot{M}_{\text{out}}([\text{O III}]) = 60 \text{ M}_{\odot} \text{ yr}^{-1}$ and $\dot{M}_{\text{out}}(H\beta) = 200 \text{ M}_{\odot} \text{ yr}^{-1}$, respectively$^4$. We note that the [O III]-based estimates are a few times lower than those obtained with H$\beta$, in agreement with previous results suggesting that [O III] provides a lower limit on the ionised gas mass (see e.g. Carniani et al. 2015; Perna et al. 2015; Marshall et al. 2023; Venturi et al. 2023). Therefore, in the following, we refer only to the H$\beta$-based outflow energetics for the emitting gas component of the GS133 outflow. The inferred outflow mass rate, kinetic power, and momentum power (defined as $\dot{P} = \dot{M}v$), are reported in Table 5. In Sect. 7.5 we compare these quantities with the AGN bolometric luminosity and some theoretical predictions to infer the driving mechanisms of the mini-BAL and [O III] outflows detected in GS133.

## 7.4. 3D modelling of the outflowing emitting gas

Our CLOUDY models suggest that the mini-BAL absorbing material, traced by C II, C IV, and other atomic species detected in the VIMOS spectrum of GS133, extends to kpc scales. Since we observe the mini-BAL directly along the LOS, we aim to determine whether this absorbing outflow is partially mixed with the emitting [O III] gas.

Despite the fact that the [O III] outflow in GS133 extends over kiloparsec-scales, it only covers a few tens of spaxels in the NIRSpec cube (Fig. 9). A detailed reconstruction of the 3D geometry of the outflow, as done for targets at lower redshifts (e.g. Meena et al. 2021; Cresci et al. 2023; Ulivi et al. 2025), would require finer spatial sampling and higher angular resolution than what NIRSpec IFS provides. Nonetheless, we used our framework ‘Modelling Outflows and Kinematics of AGN in 3D’ (MOKA$^{3D}$; Marconcini et al. 2023) to test whether our NIRSpec data are compatible with the presence of a biconical outflow where the approaching side intersects our LOS.

For simplicity, we adopted constant radial velocity profiles for both the approaching and receding components of

<table><caption>**Table 5.** Mini-BAL and [O III] outflow energetics.</caption><thead><tr><th>Measurement</th><th>Value</th></tr></thead><tbody><tr><td>L<sub>bol</sub>(X-ray)</td><td>1.6 × 10<sup>45</sup> erg s<sup>-1</sup></td></tr><tr><td>L<sub>bol</sub>(Hβ)</td><td>4.5 × 10<sup>46</sup> erg s<sup>-1</sup></td></tr><tr><td>L<sub>bol</sub>([O III])</td><td>3.0 × 10<sup>46</sup> erg s<sup>-1</sup></td></tr><tr><td>L<sub>bol</sub>(SED) (<sup>†</sup>)</td><td>2.0 × 10<sup>45</sup> erg s<sup>-1</sup></td></tr><tr><td>L<sub>bol</sub>(SED) (<sup>•</sup>)</td><td>4.9 × 10<sup>45</sup> erg s<sup>-1</sup></td></tr><tr><td>M&#775;(Hβ)</td><td>200 M<sub>☉</sub> yr<sup>-1</sup></td></tr><tr><td>E&#775;(Hβ)</td><td>3 × 10<sup>43</sup> erg s<sup>-1</sup></td></tr><tr><td>E&#775;(Hβ)/L<sub>bol</sub></td><td>2%</td></tr><tr><td>P&#775;(Hβ)</td><td>10<sup>36</sup> dyne</td></tr><tr><td>P&#775;(Hβ)/(L<sub>bol</sub>/c)</td><td>24</td></tr><tr><td>M&#775;(UV)<sub>1 kpc</sub></td><td>20 M<sub>☉</sub> yr<sup>-1</sup></td></tr><tr><td>E&#775;(UV)<sub>1 kpc</sub></td><td>2 × 10<sup>43</sup> erg s<sup>-1</sup></td></tr><tr><td>E&#775;(UV)<sub>1 kpc</sub>/L<sub>bol</sub></td><td>0.1%</td></tr><tr><td>P&#775;(UV)<sub>1 kpc</sub></td><td>2 × 10<sup>35</sup> dyne</td></tr><tr><td>P&#775;(UV)<sub>1 kpc</sub>/(L<sub>bol</sub>/c)</td><td>4</td></tr><tr><td>M&#775;(UV)<sub>10 kpc</sub></td><td>200 M<sub>☉</sub> yr<sup>-1</sup></td></tr><tr><td>E&#775;(UV)<sub>10 kpc</sub></td><td>2 × 10<sup>44</sup> erg s<sup>-1</sup></td></tr><tr><td>E&#775;(UV)<sub>10 kpc</sub>/L<sub>bol</sub></td><td>1%</td></tr><tr><td>P&#775;(UV)<sub>10 kpc</sub></td><td>2 × 10<sup>36</sup> dyne</td></tr><tr><td>P&#775;(UV)<sub>10 kpc</sub>/(L<sub>bol</sub>/c)</td><td>40</td></tr></tbody></table>

**Notes.** Kinetic power coupling efficiencies and momentum bursts are computed with the bolometric luminosity inferred from X-rays, with a bolometric correction of 260 (see Sect. 7); for the mini-BAL, we considered two potential outflow extents, 1 and 10 kpc, which align with the range determined from the CLOUDY models in Sect. 7.1. (<sup>†</sup>)SED fitting value from Circosta et al. (in prep.). (<sup>•</sup>)SED fitting value from Guo et al. (2021).

the [O III] outflow. We also assumed a bi-conical geometry, with a semi-aperture angle of 45°, consistent with other AGN-driven outflows studied at lower redshifts (e.g. Fischer et al. 2014; Müller-Sánchez et al. 2016; Meena et al. 2021; Perna et al. 2022). For the approaching cone, we required an inclination angle with respect to the LOS in the range [−45°, +45°], to ensure the overlap with our LOS; for the receding cone, we required an inclination angle in the range [180° − 45°, 180° + 45°].

The [O III] outflow is reproduced with MOKA$^{3D}$ with an approaching cone with an inclination angle of 40° with respect to the LOS, a position angle of 50$^\circ$, and an intrinsic deprojected outflow velocity of 900 km s$^{-1}$; the receding cone has an inclination angle of 225°, a position angle of 230°, and an intrinsic deprojected outflow velocity of 800 km s$^{-1}$. The two cones extend out to 5 kpc. Figure 12 shows the 3D representation of the outflow, together with the flux, velocity, and velocity dispersion maps from data, the 3D model, and the residuals; the figure indicates the MOKA$^{3D}$ capabilities in recovering the general properties of [O III] outflowing gas with a bi-conical model. In particular, a hollow conical best-fit model is obtained, consistent with many other AGN-driven outflows reported in the literature (e.g. Marconcini et al. 2023 and references therein). We also note that the MOKA$^{3D}$ model necessitates the presence of emission on kiloparsec scales located close to the plane of the sky. This component exhibits relatively low projected velocities and can explain the slightly redshifted emission observed towards the NE and the blueshifted emission towards the SW in the narrow

<footnote>$^4$ A factor of two higher estimates can be obtained considering $n_e = 400 \pm 90 \text{ cm}^{-3}$, obtained using the ionisation parameter method (Baron & Netzer 2019).</footnote>

<footnote>$^5$ The position angles are measured anti-clockwise from the north.</footnote>
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** MOKA$^{3D}$ model (left) and moment maps from NIRSpec data and MOKA$^{3D}$ model (right panels). In the 3D representations of the GS133 outflow structure, the XY represents the plane of the sky, while Z axis is the LOS. The observer is positioned on the left, and sees as blueshifted (redshifted) all the approaching (receding) gas, according to the colour bar; largest bubbles identify brightest clumps. The three-by-three panels show the comparison between the flux, Mom-1 and Mom-2 maps of [O III] (top panels) and MOKA$^{3D}$ model (middle panels) for the GS133 bi-conical outflow. The bottom panels present the residuals obtained subtracting the model from the data.</figcaption></figure>

component maps shown in Fig. 8 which cannot be attributed to a rotating pattern.

From this simple MOKA$^{3D}$ model, we can derive two main results. On the one hand, the [O III] gas appears to share the same velocity as the low-velocity component in C II, C IV, and other UV lines ($|v_{10\text{out}}| \sim v_{\text{MOKA}} \approx 900 \text{ km s}^{-1}$); this suggests that both the absorbing gas and part of the emitting [O III] could lie along our LOS and may be associated with the same outflow, possibly even physically mixed. On the other hand, the second kinematic component detected in absorption in C IV and other HILs moves at $v_{\text{out}} \sim -1900 \text{ km s}^{-1}$ (see Fig. 4); this component is not detected in the ionised emitting gas, and could therefore have slightly different physical (e.g. composition, ionisation state) or structural (e.g. spatial location) properties. Alternatively, this very high-velocity component may be too faint to be detected in emission.

## 7.5. Physical mechanisms for the wind propagation

We compared our inferred values of the total outflow kinetic power with the AGN bolometric luminosity and the expected kinetic power ascribed to stellar processes. The AGN bolometric luminosity can be constrained from the X-ray emission, applying a Duras et al. (2020) bolometric correction of 260 to the absorption corrected 2–10 keV luminosity from Li et al. (2019), $L_{\text{bol}} = 1.6 \times 10^{45} \text{ erg s}^{-1}$. This value is broadly consistent with the SED-based AGN luminosity inferred by Guo et al. (2021, $L_{\text{bol}} = 4.9 \times 10^{45} \text{ erg s}^{-1}$), and by Circosta et al. (in prep., $2 \times 10^{45} \text{ erg s}^{-1}$). However, these values are significantly lower than those inferred from the outflow component of the optical emission lines: from H$\beta$ flux, applying the Netzer (2019) bolometric correction, $L_{\text{bol}}(\text{H}\beta) = 4.5 \times 10^{46} \text{ erg s}^{-1}$, and from [O III] flux, applying the Heckman et al. (2004) correction, $L_{\text{bol}}([\text{O III}]) = 3 \times 10^{46} \text{ erg s}^{-1}$ (even higher values would be obtained applying a

dust extinction correction, and considering the total line profiles, which however could be contaminated by star-formation processes). Since the X-ray-based value is in between those from the SED analysis, in the following, we refer to the former as our fiducial estimate.

The kinetic power associated with the optical outflow is $3 \times 10^{42} \text{ erg s}^{-1}$, hence ~0.2% of the radiative luminosity of the AGN, consistent with the energetics of other AGN in the literature (e.g. Harrison et al. 2018). Comparable values can be obtained for the mini-BAL considering a kiloparsec-scale extension: $\dot{E}_{\text{out}}/L_{\text{bol}} \sim 0.1\%$ for $R_{\text{out}} = 1 \text{ kpc}$, and ~1% for $R_{\text{out}} = 10 \text{ kpc}$.

Following Veilleux et al. (2005), we assumed an expected kinetic power ascribed to stellar processes proportional to the SFR ($20 \text{ M}_\odot \text{ yr}^{-1}$, from Guo et al. 2020, $<70 \text{ M}_\odot \text{ yr}^{-1}$, from Circosta et al., in prep.), and concluded that the outflow kinetic power of the emitting gas cannot be associated with star formation driven winds, as they would require an extreme coupling between the stellar processes and the observed winds, of the order of ≈100% or more.

The momentum power inferred for the optical emitting gas is $\dot{P}_{\text{H}\beta} = 1 \times 10^{36} \text{ dyne}$, with a momentum rate $\dot{P}_{\text{H}\beta}/(L_{\text{bol}}/c) = 24$. For the UV absorbing gas, we obtained ($2-20 \times 10^{35} \text{ dyne}$, resulting in $\dot{P}_{\text{UV}}/(L_{\text{bol}}/c) = 4-40$, with lowest values occurring for a mini-BAL with $R_{\text{out}} = 1 \text{ kpc}$ and highest for $R_{\text{out}} = 10 \text{ kpc}$ (still compatible with measurements of other AGN in the literature, see e.g. Fig. 3 in Bischetti et al. 2024).

A momentum flux rate obtained under the assumption of a very compact mini-BAL is far smaller than the moment flux rate of the H$\beta$ outflow ($\dot{P}_{\text{UV}}/(L_{\text{bol}}/c) = 0.004$ for $R_{\text{out}} = 1 \text{ pc}$), and hence than what is expected for a momentum conserving wind (e.g. Fiore et al. 2017). If the AGN outflow is energy conserving, it is possible that the mini-BAL is radiatively launched with small momentum flux in the nucleus, and
<figure><img src="image_14.png" /></figure>

**Fig. 13.** JWST/NIRCam F115W/F200W/F444W cutout showing the close environment of GS133. The red dashed box marks the NIRSpec IFS FOV, while the white contours show the Chandra 0.5–7 keV emission from Luo et al. (2017), highlighting the position of two X-ray AGN with a projected separation of ~130 kpc and a velocity offset of ~5000 km s⁻¹ (Δz = 0.08). The apparent offset between the NIR and X-ray positions can be attributed to Chandra’s low angular resolution and the fact that the Chandra image is composed of multiple observations taken with varying aim points and roll angles. The inset shows a zoom-in on GS133. The images are oriented north up, east left.

then boosted to $\dot{P}_{H\beta}/(L_{bol}/c) \sim 24$ at larger distances (e.g. Faucher-Giguère & Quataert 2012). However, such a dramatic increase in momentum flux, by four orders of magnitude, is highly uncommon (see Fig. 8 in Tozzi et al. 2021). This further supports our CLOUDY model results indicating that the mini-BAL is located at kpc scales.

Summarising, the computation of outflow energetics allowed us to obtain an independent (but still indirect) confirmation of the kiloparsec-scale location of the mini-BAL; the modest momentum bursts inferred for optical and UV transitions, of the order of 2 to 40, are more likely compatible with energy conserving winds rather than momentum driven ones (e.g. King & Pounds 2015), also taking into account that with current data we cannot probe the neutral and molecular outflow phases which could significantly contribute to the outflow energetics on kiloparsecscales (e.g. Brusa et al. 2018; Cresci et al. 2023; Davies et al. 2024; D'Eugenio et al. 2024).

## 8. Environment

NIRSpec IFS data do not reveal any close companions or line emitters around GS133, in contrast with several $z > 3$ AGN and QSOs recently observed with JWST, many of which show multiple nearby companions (e.g. Wylezalek et al. 2022; Marshall et al. 2023, 2024; Perna et al. 2023a,b; Übler et al. 2023, 2024a,b; Kashino et al. 2023; Yang et al. 2023). Available JWST/NIRCam images (programme ID 3990, PI: T. Morishita) show some extended emission towards E-SE associated with the GS133 host, possibly indicating merger features or a disturbed morphology, and faint compact sources outside of the NIRSpec FOV. Figure 13 shows the composite F115W/F200W/F444W cutout, highlighting the GS133 extended features. However, these faint features remain undetected in the spectroscopic data,

which are too shallow ($\le 1$ h-integration) in comparison with NIRCam imaging (e.g. 6 h-integration in the red band F444W).

Interestingly, GS133 is located near another AGN, identified as CDFS019505 (in the bottom right of Fig. 13), targeted by VLT/VIMOS as part of the VANDELS project and classified as a narrow-line AGN in Mascia et al. (2023). CDFS019505 was observed with a total integration time of 40 hours, and its spectroscopic redshift, $z = 3.396$, is considered insecure, as it is based on a single emission line at ~5350 Å, assumed to be $Ly\alpha$ (see Fig. 7, bottom panel, in Carnall et al. 2020). On the other hand, this redshift is close to the photometric estimates (e.g. $z_{phot} = 3.417$ in Hsu et al. 2014, $z_{phot} = 3.49_{-0.27}^{+0.09}$ in Carnall et al. 2020), and it is supported by the detection of continuum emission at observed wavelengths >5350 Å. SED analysis by Carnall et al. (2020) classified CDFS019505 as a quiescent galaxy, requiring an escape fraction $f_{esc, Ly\alpha} > 0.1$, with a $Ly\alpha$-based SFR of $0.28 \pm 0.05/f_{esc, Ly\alpha} M_{\odot} yr^{-1}$ (similar SED-based value is reported in Guo et al. 2020, SFR = $1 M_{\odot} yr^{-1}$). This source is also classified as an AGN on the basis of its X-ray luminosity, with 140 full (0.5–7 keV) aperture-corrected source counts, an intrinsic column density $N_H = 10^{24} \text{cm}^{-2}$, and an absorption-corrected intrinsic luminosity $L_{2-10\text{keV}} = 10^{44} \text{ erg s}^{-1}$ (Li et al. 2019).

With a projected separation of 17" (i.e. ~130 kpc) and a velocity offset of ~5000 km s⁻¹ relative to GS133, CDFS019505 displays a significantly larger separation and redshift difference<fnref n="6" /> than is typical for close AGN pairs, or ‘dual AGN’ (e.g. De Rosa et al. 2019; Mannucci et al. 2023; Scialpi et al. 2024; Chen et al. 2024). Perna et al. (2023b) reported the discovery of three dual AGN systems in the GA-NIFS sample, with much closer projected separations of 5–10 kpc and velocity offsets of a few hundred km s⁻¹, suggesting close interactions likely to result in mergers. In contrast, the spatial and redshift separation between GS133 and CDFS019505 makes a direct interaction between these two AGN improbable at the time of observation. Nonetheless, the identification of another AGN pair among the 12 AGN systems at $z \sim 3$ targeted in GA-NIFS further supports the hypothesis that AGN at high redshift tend to reside in dense environments (e.g. Karouzos et al. 2014), potentially hosting multiple accreting SMBHs within close proximity (see also e.g. Coogan et al. 2018; Pensabene et al. 2024; Zamora et al. 2024).

## 9. Conclusions

We have presented JWST NIRSpec IFS data of the Compton thick quasar GS133 at $z = 3.46$. These observations allowed us to map the extension of the quasar host traced by rest-frame optical emission lines. NIRSpec data have been complemented with 1D spectroscopic data from VLT/VIMOS, covering the rest-frame UV lines, and NIRCam images, to explore the GS133 environment. These are our main results.

We identified mini-BAL features in multiple transitions covered by VIMOS data, including N V, Si IV, C IV, C II. The column densities derived from VoigtFit for individual metal species span a narrow range of values, $log(N/cm^{-2}) = 18.3–19.4$. Higher ionisation lines (e.g. C IV and Si IV) are broader than the lower ionisation lines (C II), and show two kinematic components, shifted by −800 km s⁻¹ and −1900 km s⁻¹ respectively. The lower-ionisation lines only reveal the −800 km s⁻¹ component. Both kinematic components are tentatively detected (S/N ~ 2) in absorption in the $H\alpha$ line observed with NIRSpec.

<footnote n="6">⁶ We assumed the spectroscopic redshift $z = 3.396$ for CDFS019505.</footnote>
Using single-cloud CLOUDY models, we derived hydrogen
density and column density estimates for the mini-BALs to con-
strain the distance of the absorbing material from the AGN, and,
therefore, the outflow energetics: log (N<sub>H</sub>/cm<sup>-2</sup>) = 19.4-20.3,
log (n<sub>H</sub>/cm<sup>-3</sup>) = 1-3, with R<sub>out</sub> = 1-10 kpc. Outflow energetics
were derived assuming a thin-shell approximation, and account-
ing for uncertainties in the best-fit CLOUDY parameters: we
obtained a outflow mass rate of 20-200 M<sub>☉</sub> yr<sup>-1</sup>, and a kinetic
power in the range of 2-20 × 10<sup>43</sup> erg s<sup>-1</sup>.

Spaxel-by-spaxel analysis of NIRSpec data allowed us to
perform emission lines kinematic decompositions. We revealed
two distinct components in the host galaxy of GS133: the first
component likely traces a rotating disc with a dynamical mass
~2 × 10<sup>10</sup> M<sub>☉</sub>. The second component corresponds to a galaxy-
wide, bi-conical outflow, with a velocity of ~±1000 km s<sup>-1</sup> and
a projected radius of ~3 kpc.

The inferred outflow velocity for the [O III] line (W80<sub>out</sub> ~ 1300 km s<sup>-1</sup>, v10<sub>out</sub> ~ -900 km s<sup>-1</sup>) aligns with the general L<sub>AGN</sub> - v<sub>out</sub> trends reported in the literature (Musiimenta et al. 2023), placing GS133 at the higher end of the velocity distribution. This is consistent with findings by Tozzi et al. (2024) showing that CT AGN tend to exhibit faster outflows than unobscured AGN at fixed L<sub>bol</sub>.

We derived emitting-gas outflow energetics assuming a
simple bi-conical configuration: the outflow mass rate is
200 M<sub>☉</sub> yr<sup>-1</sup>, and the kinetic power is 3 × 10<sup>43</sup> erg s<sup>-1</sup>, compa-
rable to those inferred for the mini-BALs in the UV regime.
The energetics of both absorbing and emitting gas suggest
that the outflow is driven by energy conserving winds, with
E&#775;/L<sub>bol</sub> = 0.1-2%, P&#775;/ (L<sub>bol</sub>/c) = 4-40, and mass loading factor
M&#775;/SFR = 1-10.

We explored the potential overlap of absorbing and emit-
ting outflows using MOKA<sup>3D</sup>, modelling the [O III] outflow
with a bi-conical structure intersecting our LOS. The MOKA<sup>3D</sup>
model fits the data well with the following best-fit parame-
ters: a semi-aperture angle of 45°, a radius of 5 kpc, a con-
stant radial velocity of 800-900 km s<sup>-1</sup>, and an inclination of
40° with respect to our LOS, hence with an approaching cone
overlapping the LOS. Similarities in velocity, location of the
absorbing gas, and LOS extension of the emitting gas, suggest
at least partial mixing between the mini-BAL and [O III] out-
flows. However, the faster mini-BAL component seen in C IV
and other high-ionisation lines is not present in [O III], sug-
gesting that this UV component might have different prop-
erties, for instance, in terms of gas composition and spatial
location.

NIRSpec IFS data reveal no close companions or line emitters around GS133 within the ~3'' × 3'' FOV. However, JWST/NIRCam images suggest possible merger features that could help explain the Compton thick nature of this AGN and the onset of outflow episodes. We also reported the presence of an X-ray-detected AGN at ~130 kpc and a velocity offset of ~5000 km s⁻¹ from GS133. These results are aligned with recent JWST findings that high-z AGN commonly reside in dense environments.

This comprehensive study demonstrates the power of com-
bining UV and optical data to probe the complex outflows in
high-redshift AGN, providing insights into their origin and asso-
ciated feedback processes. NIRSpec data allowed us to spa-
tially resolve both unperturbed and ejected emitting gas, and
obtain crucial information about the 3D configuration of a bi-
conical outflow at z ~ 3.5. With less than 1 hour of total expo-
sure, JWST NIRSpec IFS provided essential complementary
data to the extremely deep VLT/VIMOS observations (41-hour

integration), demonstrating JWST's ability to efficiently capture
detailed outflow properties.

**Acknowledgements.** We thank the referee for a careful reading of the manuscript and a constructive report. We are grateful to Margherita Giustini for discussing various aspects of this work. MP acknowledges support through the grant RYC2023-044853-I, funded by MICIU/AEI/10.13039/501100011033 and FSE+. MP, SA, BRP, and IL acknowledge support from the research project PID2021-127718NB-I00 of the Spanish Ministry of Science and Innovation/State Agency of Research (MCIN/AEI/10.13039/501100011033). IL acknowledges support from grant PRIN-MUR 2020ACSP5K_002 financed by European Union – Next Generation EU. RM and XJ acknowledge support by the Science and Technology Facilities Council (STFC), by the ERC Advanced Grant 695671 ‘QUENCH’, and by the UKRI Frontier Research grant RISEandFALL; RM is further supported by a research professorship from the Royal Society. AJB acknowledges funding from the ‘First Galaxies’ Advanced Grant from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 789056). MP, GC and EB acknowledge the support of the INAF Large Grant 2022 “The metal circle: a new sharp view of the baryon cycle up to Cosmic Dawn with the latest generation IFU facilities”. HÜ acknowledges support through the ERC Starting Grant 101164796 ‘APEX’.

**References**

Abazajian, K. N., Adelman-McCarthy, J. K., Agüeros, M. A., et al. 2009, ApJS, 182, 543

Ahmed, H., Shemmer, O., Matthews, B., et al. 2024, ApJ, 968, 77

Arav, N., Borguet, B., Chamberlain, C., Edmonds, D., & Danforth, C. 2013, MNRAS, 436, 3286

Arribas, S., Colina, L., Bellocchi, E., Maiolino, R., & Villar-Martín, M. 2014, A&A, 568, A14

Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5

Balestra, I., Mainieri, V., Popesso, P., et al. 2010, A&A, 512, A12

Baron, D., & Netzer, H. 2019, MNRAS, 486, 4290

Bertola, E., Circosta, C., Ginolfi, M., et al. 2024, A&A, 691, A178

Bischetti, M., Choi, H., Fiore, F., et al. 2024, ApJ, 970, 9

Bökér, T., Arribas, S., Lützgendorf, N., et al. 2022, A&A, 661, A82

Bordoloi, R., Lilly, S. J., Hardmeier, E., et al. 2014, ApJ, 794, 130

Bruni, G., Piconcelli, E., Misawa, T., et al. 2019, A&A, 630, A111

Brusa, M., Cresci, G., Daddi, E., et al. 2018, A&A, 612, A29

Calabrò, A., Pentericci, L., Feltre, A., et al. 2023, A&A, 679, A80

Cano-Díaz, M., Maiolino, R., Marconi, A., et al. 2012, A&A, 537, L8

Cappellari, M. 2017, MNRAS, 466, 798

Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245

Carnall, A. C., Walker, S., McLure, R. J., et al. 2020, MNRAS, 496, A695

Carniani, S., Marconi, A., Maiolino, R., et al. 2015, A&A, 580, A102

Chabrier, G. 2003, PASP, 115, 763

Chartas, G., Cappi, M., Vignali, C., et al. 2021, ApJ, 920, 24

Chen, Y. C., Gross, A. C., Liu, X., et al. 2024, ApJ, submitted
[arXiv:2409.16351]

Choi, H., Leighly, K. M., Dabbieri, C., et al. 2022, ApJ, 936, 110

Cicone, C., Maiolino, R., Sturm, E., et al. 2014, A&A, 562, A21

Ciesla, L., Charmandaris, V., Georgakakis, A., et al. 2015, A&A, 576, A10

Comastri, A. 2004, Astrophys. Space Sci. Lib., 308, 245

Coogan, R. T., Daddi, E., Sargent, M. T., et al. 2018, MNRAS, 479, 703

Constantin, L., Gillman, S., Boogaard, L. A., et al. 2024, A&A, submitted
[arXiv:2407.00153]

Cresci, G., Mainieri, V., Brusa, M., et al. 2015, ApJ, 799, 82

Cresci, G., Tozzi, G., Perna, M., et al. 2023, A&A, 672, A128

Cristiani, S., Appenzeller, I., Arnouts, S., et al. 2000, A&A, 359, 489

Dai, X., Hao, L., Jiang, P., et al. 2024, ApJ, 969, 99

Davies, R. L., Belli, S., Park, M., et al. 2024, MNRAS, 528, 4976

Decarli, R., Walter, F., Venemans, B. P., et al. 2018, ApJ, 854, 97

Decarli, R., Loiacono, F., Farina, E. P., et al. 2024, A&A, 689, A219

Del Zanna, G., Landini, M., & Mason, H. E. 2002, A&A, 385, 968

De Rosa, A., Vignali, C., Bogdanović, T., et al. 2019, New Astron. Rev., 86, 101525

D'Eugenio, F., Pérez-González, P. G., Maiolino, R., et al. 2024, Nat. Astron., 8, 1443

Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604

Dunn, J. P., Bautista, M., Arav, N., et al. 2010, ApJ, 709, 611

Duras, F., Bongiorno, A., Ricci, F., et al. 2020, A&A, 636, A73

Falcone, J., Crenshaw, D. M., Fischer, T. C., et al. 2024, ApJ, 971, 17

Faucher-Giguère, C.-A., & Quataert, E. 2012, MNRAS, 425, 605

Feltre, A., Charlot, S., & Gutkin, J. 2016, MNRAS, 456, 3354
Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, *Rev. Mex. Astron. Astrofis.*, 53, 385

Fiore, F., Feruglio, C., Shankar, F., et al. 2017, A&A, 601, A143

Fischer, T. C., Crenshaw, D. M., Kraemer, S. B., Schmitt, H. R., & Turner, T. J. 2014, ApJ, 785, 25

Förster Schreiber, N. M., Renzini, A., Mancini, C., et al. 2018, ApJS, 238, 21

Förster Schreiber, N. M., Übler, H., Davies, R. L., et al. 2019, ApJ, 875, 21

Garilli, B., McLure, R., Pentericci, L., et al. 2021, A&A, 647, A150

Genzel, R., Förster Schreiber, N. M., Rosario, D., et al. 2014, ApJ, 796, 7

Giavalisco, M., Ferguson, H. C., Koekemoer, A. M., et al. 2004, ApJ, 600, L93

Giustini, M., & Proga, D. 2019, A&A, 630, A94

Grevesse, N., Asplund, M., Sauval, A. J., & Scott, P. 2010, Ap&SS, 328, 179

Gunasekera, C. M., Ji, X., Chatzikos, M., Yan, R., & Ferland, G. 2023, MNRAS, 520, 4345

Guo, X., Gu, Q., Ding, N., Contini, E., & Chen, Y. 2020, MNRAS, 492, 1887

Guo, X., Gu, Q., Ding, N., Yu, X., & Chen, Y. 2021, ApJ, 908, 169

Hamann, F., & Sabra, B. 2004, ASP Conf. Ser., 311, 203

Hamann, F., Herbst, H., Paris, I., & Capellupo, D. 2019, MNRAS, 483, 1808

Harrison, C. M., Alexander, D. M., Swinbank, A. M., et al. 2012, MNRAS, 426, 1073

Harrison, C. M., Johnson, H. L., Swinbank, A. M., et al. 2017, MNRAS, 467, 1965

Harrison, C. M., Costa, T., Tadhunter, C. N., et al. 2018, Nat. Astron., 2, 198

Heckman, T. M., Kauffmann, G., Brinchmann, J., et al. 2004, ApJ, 613, 109

Hopkins, P. F., Hernquist, L., Cox, T. J., & Kereš, D. 2008, ApJS, 175, 356

Hsu, L.-T., Salvato, M., Nandra, K., et al. 2014, ApJ, 796, 60

Hutchings, J. B., Crenshaw, D. M., & Kraellup, S. B., et al. 2002, AJ, 124, 2543

Ishibashi, W., Fabian, A. C., & Hewett, P. C. 2024, MNRAS, 533, 4384

Jakobsen, P., Ferruit, P., de Oliveira, C. A., et al. 2022, A&A, 661, A80

Jenkins, E. B. 2009, ApJ, 700, 1299

Juoždžalis, I., Ji, X., Maiolino, R., et al. 2024, MNRAS, 535, 853

Kakkad, D., Mainieri, V., Vietri, G., et al. 2020, A&A, 642, A147

Karouzos, M., Jarvis, M. J., & Bonfield, D. 2014, MNRAS, 439, 861

Kartaltepe, J. S., Rose, C., Vanderhoof, B. N., et al. 2023, ApJ, 946, L15

Kashino, D., Lilly, S. J., Matthee, J., et al. 2023, ApJ, 950, 66

Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055

Keenan, F. P., Feibelman, W. A., & Berrington, K. A. 1992, ApJ, 389, 443

Kehoe, E., Shapley, A. E., Schreiber, N. M. F., et al. 2024, ApJ, 976, 28

Kewley, L. J., Dopita, M. A., Sutherland, R. S., Heisler, C. A., & Trevena, J. 2001, ApJ, 556, 121

King, A., & Pounds, K. 2015, ARA&A, 53, 115

Kocevski, D. D., Finkelstein, S. L., Barro, G., et al. 2024, ApJ, submitted [arXiv:2404.03576]

Krogager, J. K. 2018, ArXiv e-prints [arXiv:1803.01187]

Laha, S., Reynolds, C. S., Reeves, J., et al. 2021, Nat. Astron., 5, 13

Lamperti, I., Pereira-Santaella, M., Perna, M., et al. 2022, A&A, 668, A45

Lamperti, I., Arribas, S., Perna, M., et al. 2024, A&A, 691, A153

Lanzuisi, G., Perna, M., Delvecchio, I., et al. 2015, A&A, 578, A120

Lanzuisi, G., Delvecchio, I., Berta, S., et al. 2017, A&A, 602, A123

Lanzuisi, G., Matzeu, G., Baldini, P., et al. 2024, A&A, 689, A247

Li, J., Xue, Y., Sun, M., et al. 2019, ApJ, 877, 5

Liu, G., Zakamska, N. L., Greene, J. E., Nesvadba, N. P. H., & Liu, X. 2013, MNRAS, 436, 2576

Liu, G., Arav, N., & Rupke, D. S. N. 2015, ApJS, 221, 9

Luo, B., Brandt, W. N., Xue, Y. Q., et al. 2017, ApJS, 228, 2

Maeder, A., Meynet, G., & Chiappini, C. 2015, A&A, 576, A56

Maiolino, R., & Mannucci, F. 2019, A&ARv, 27, 3

Maiolino, R., Salvati, M., Bassani, L., et al. 1998, A&A, 338, 781

Maiolino, R., Scholtz, J., Witstok, J., et al. 2024a, Nature, 627, 59

Maiolino, R., Übler, H., Perna, M., et al. 2024b, A&A, 687, A67

Mannucci, F., Scialpi, M., Ciurlo, A., et al. 2023, A&A, 680, A53

Marconcini, C., Marconi, A., Cresci, G., et al. 2023, A&A, 677, A58

Marconi, A., Amiri, A., Feltre, A., et al. 2024, A&A, 689, A78

Marshall, M. A., Perna, M., Willott, C. J., et al. 2023, A&A, 678, A191

Marshall, M. A., Yue, M., Eilers, A. C., et al. 2024, A&A, submitted [arXiv:2410.11035]

Mascia, S., Pentericci, L., Saxena, A., et al. 2023, A&A, 674, A221

Matthee, J., Naidu, R. P., Brammer, G., et al. 2024, ApJ, 963, 129

Matzeu, G. A., Brusa, M., Lanzuisi, G., et al. 2023, A&A, 670, A182

McLure, R. J., Pentericci, L., Cimatti, A., et al. 2018, MNRAS, 479, 25

Meena, B., Crenshaw, D. M., Schmitt, H. R., et al. 2021, ApJ, 916, 31

Mehdipour, M., Kriss, G. A., Brusa, M., et al. 2023, A&A, 670, A183

Moravec, E. A., Hamann, F., Capellupo, D. M., et al. 2017, MNRAS, 468, 4539

Mullaney, J. R., Alexander, D. M., Fine, S., et al. 2013, MNRAS, 433, 622

Müller-Sánchez, F., Comerford, J., Stern, D., & Harrison, F. A. 2016, ApJ, 830, 50

Musiimenta, B., Brusa, M., Liu, T., et al. 2023, A&A, 679, A84

Nakajima, K., Schaerer, D., Le Fèvre, O., et al. 2018, A&A, 612, A94

Neeleman, M., Novak, M., Venemans, B. P., et al. 2021, ApJ, 911, 141

Netzer, H. 2019, MNRAS, 488, 5185

Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei (Sausalito: University Science Books)

Parlanti, E., Carniani, S., Venturi, G., et al. 2024a, A&A, submitted [arXiv:2407.19008]

Parlanti, E., Carniani, S., Übler, H., et al. 2024b, A&A, 684, A24

Pensabene, A., Cantalupo, S., Cicone, C., et al. 2024, A&A, 684, A119

Pereira-Santaella, M., Colina, L., García-Burillo, S., et al. 2018, A&A, 616, A171

Perna, M., Brusa, M., Cresci, G., et al. 2015, A&A, 574, A82

Perna, M., Lanzuisi, G., Brusa, M., Mignoli, M., & Cresci, G. 2017, A&A, 603, A99

Perna, M., Sargent, M. T., Brusa, M., et al. 2018, A&A, 619, A90

Perna, M., Cresci, G., Brusa, M., et al. 2019, A&A, 623, A171

Perna, M., Arribas, S., Catalán-Torrecilla, C., et al. 2020, A&A, 643, A139

Perna, M., Arribas, S., Colina, L., et al. 2022, A&A, 662, A94

Perna, M., Arribas, S., Marshall, M., et al. 2023a, A&A, 679, A89

Perna, M., Arribas, S., Lamperti, I., et al. 2023b, A&A, submitted [arXiv:2310.03067]

Perrotta, S., Hamann, F., Zakamska, N. L., et al. 2019, MNRAS, 488, 4126

Pounds, K. A., Reeves, J. N., King, A. R., et al. 2003, MNRAS, 345, 705

Rauscher, B. J., Arendt, R. G., Fixsen, D. J., et al. 2017, PASP, 129, 105003

Rigby, J., Perrin, M., McElwain, M., et al. 2023, PASP, 135, 048001

Rodríguez Del Pino, B., Perna, M., Arribas, S., et al. 2024, A&A, 684, A187

Saturni, F. G., Revese, D., Vagnetti, F., Perna, M., & Dadina, M. 2016, A&A, 587, A43

Scholtz, J., Maiolino, R., D'Eugenio, F., et al. 2023, A&A, submitted [arXiv:2311.18731]

Scholtz, J., Curti, M., D'Eugenio, F., et al. 2024, MNRAS, submitted [arXiv:2411.07695]

Schulze, A., Schramm, M., Zuo, W., et al. 2017, ApJ, 848, 104

Schulze, A., Misawa, T., Zuo, W., & Wu, X.-B. 2018, ApJ, 853, 167

Schwarz, U. J. 1978, A&A, 65, 345

Scialpi, M., Mannucci, F., Marconcini, C., et al. 2024, A&A, 690, A57

Shen, Y., Richards, G. T., Strauss, M. A., et al. 2011, ApJS, 194, 45

Shi, X., Zhou, H., Shu, X., et al. 2016, ApJ, 819, 99

Stepney, M., Banerji, M., Tang, S., et al. 2024, MNRAS, 533, 2948

Temple, M. J., Rankine, A. L., Banerji, M., et al. 2024, MNRAS, 532, 424

Tian, Q., Shi, X., Lu, H., et al. 2019, ApJ, 877, 72

Topping, M. W., Stark, D. P., Senchyna, P., et al. 2024, ApJ, submitted [arXiv:2407.19009]

Tozzi, G., Cresci, G., Marasco, A., et al. 2021, A&A, 648, A99

Tozzi, G., Cresci, G., Perna, M., et al. 2024, A&A, 690, A141

Übler, H., Maiolino, R., Curtis-Lake, E., et al. 2023, A&A, 677, A145

Übler, H., Maiolino, R., Pérez-González, P. G., et al. 2024a, MNRAS, 531, 355

Übler, H., D'Eugenio, F., Perna, M., et al. 2024b, MNRAS, 533, 4287

Ulivi, L., Venturi, G., Cresci, G., et al. 2024, A&A, 685, A122

Ulivi, L., Perna, M., Lamperti, I., et al. 2025, A&A, 693, A36

Veilleux, S., Cecil, G., & Bland-Hawthorn, J. 2005, ARA&A, 43, 769

Venturi, G., Nardini, E., Marconi, A., et al. 2018, A&A, 619, A74

Venturi, G., Treister, E., Finlez, C., et al. 2023, A&A, 678, A127

Véron-Cetty, M. P., Véron, P., & Gonçalves, A. C. 2001, A&A, 372, 730

Vestergaard, M. 2003, ApJ, 599, 116

Vietri, G., Misawa, T., Piconcelli, E., et al. 2022, A&A, 668, A87

Villar Martín, M., Perna, M., Humphrey, A., et al. 2020, A&A, 634, A116

Walter, F., Carilli, C., Bertoldi, F., et al. 2004, ApJ, 615, L17

Wang, B., de Graaff, A., Davies, R. L., et al. 2024, ArXiv e-prints [arXiv:2403.02304]

Weymann, R. J., Carswell, R. F., & Smith, M. G. 1981, ARA&A, 19, 41

Woo, J.-H., Bae, H.-J., Son, D., & Karouzos, M. 2016, ApJ, 817, 108

Wylezalek, D., Vayner, A., Rupke, D. S. N., et al. 2022, ApJ, 940, L7

Xu, X., Zakamska, N. L., Arav, N., Miller, T., & Benn, C. 2020, MNRAS, 495, 305

Yang, J., Wang, F., Fan, X., et al. 2023, ApJ, 951, L5

Zakamska, N. L., & Greene, J. E. 2014, MNRAS, 442, 784

Zamora, S., Venturi, G., Carniani, S., et al. 2024, A&A, submitted [arXiv:2412.02751]

Zhang, S., Zhou, H., Shi, X., et al. 2015, ApJ, 815, 113
## Appendix A: Spaxel-by-spaxel multi-Gaussian fit of the NIRSpec datacube

## Appendix B: Spatially integrated spectra of off-nuclear approaching and receding gas

<table><caption>**Table B.1.** Optical emission line properties from the NIRSpec integrated spectrum of NW cone (blueshifted outflow).</caption><thead><tr><th>Measurement</th><th>Value</th></tr></thead><tbody><tr><td>Δv<sub>sys</sub></td><td>49 ± 3</td></tr><tr><td>σ<sub>sys</sub></td><td>97 ± 7</td></tr><tr><td>Δv<sub>out</sub></td><td>-75<sup>+32</sup><sub>-26</sub> km s<sup>-1</sup></td></tr><tr><td>σ<sub>out</sub></td><td>265<sup>+28</sup><sub>-18</sub> km s<sup>-1</sup></td></tr><tr><td>W80<sub>out</sub></td><td>688<sup>+58</sup><sub>-30</sub> km s<sup>-1</sup></td></tr><tr><td>v10<sub>out</sub></td><td>-440 ± 60 km s<sup>-1</sup></td></tr><tr><td>v90<sub>out</sub></td><td>250 ± 10 km s<sup>-1</sup></td></tr><tr><td>Δv<sub>tot</sub></td><td>-2 ± 4 km s<sup>-1</sup></td></tr><tr><td>σ<sub>tot</sub></td><td>196 ± 10 km s<sup>-1</sup></td></tr><tr><td>W80<sub>tot</sub></td><td>480 ± 15 km s<sup>-1</sup></td></tr><tr><td>v10<sub>tot</sub></td><td>-290 ± 10 km s<sup>-1</sup></td></tr><tr><td>v90<sub>tot</sub></td><td>190 ± 10 km s<sup>-1</sup></td></tr></tbody></table>

**Notes.** All velocity measurements refer to the [O III] line.

<table><caption>**Table B.2.** Optical emission line properties from the NIRSpec integrated spectrum of SW cone (redshifted outflow).</caption><thead><tr><th>Measurement</th><th>Value</th></tr></thead><tbody><tr><td>Δv<sub>sys</sub></td><td>-14 ± 2</td></tr><tr><td>σ<sub>sys</sub></td><td>137 ± 3</td></tr><tr><td>Δv<sub>out</sub></td><td>96<sup>+13</sup><sub>-3</sub> km s<sup>-1</sup></td></tr><tr><td>σ<sub>out</sub></td><td>311<sup>+8</sup><sub>-3</sub> km s<sup>-1</sup></td></tr><tr><td>W80<sub>out</sub></td><td>840<sup>+18</sup><sub>-10</sub> km s<sup>-1</sup></td></tr><tr><td>v10<sub>out</sub></td><td>-310<sup>+30</sup><sub>-3</sub> km s<sup>-1</sup></td></tr><tr><td>v90<sub>out</sub></td><td>485 ± 10 km s<sup>-1</sup></td></tr><tr><td>Δv<sub>tot</sub></td><td>53<sup>+6</sup><sub>-2</sub> km s<sup>-1</sup></td></tr><tr><td>σ<sub>tot</sub></td><td>260 ± 5 km s<sup>-1</sup></td></tr><tr><td>W80<sub>tot</sub></td><td>660 ± 20 km s<sup>-1</sup></td></tr><tr><td>v10<sub>tot</sub></td><td>-290 ± 10 km s<sup>-1</sup></td></tr><tr><td>v90<sub>tot</sub></td><td>360 ± 10 km s<sup>-1</sup></td></tr></tbody></table>

**Notes.** All velocity measurements refer to the [O III] line.
<figure><img src="image_15.png" /></figure>

**Fig. A.1.** Total profiles [O III] (top) and Hα (bottom) flux, moment-1 and moment-2 maps. A S/N cut of 4 has been applied to generate the maps. The disc and biconical outflow components identified in Figs. 8 and 9, respectively, are here less evident. Hα is less affected by the outflow, and its moment-1 still reveal a NW-SE gradient likely associated with a disc rotation.

<figure><img src="image_16.png" /></figure>

**Fig. B.1.** NIRSpec integrated spectrum of NE cone, showing the blueshifted outflow. The black curve identifies the integrated spectrum, integrated over a circular region identified by the light-blue crescent moon in the inset in Fig. 10; the total, multi-component best fit curve is in red, while the green and grey Gaussian components show the systemic and outflow, respectively. The fit residuals are reported in the top panels. The most prominent emission lines are marked with grey vertical lines.

<figure><img src="image_17.png" /></figure>

**Fig. B.2.** NIRSpec integrated spectrum of SW cone, showing the redshifted outflow. The black curve identifies the integrated spectrum, integrated over a circular region identified by the light-red crescent moon in the inset in Fig. 10; the total, multi-component best fit curve is in red, while the green and grey Gaussian components show the systemic and outflow, respectively. The fit residuals are reported in the top panels. The most prominent emission lines are marked with grey vertical lines.