11institutetext: Centro de Astrobiología (CAB), CSIC–INTA, Departamento de
Astrofísica, Cra. de Ajalvir Km. 4, 28850 – Torrejón de Ardoz, Madrid, Spain
22institutetext: National Research Council of Canada, Herzberg Astronomy &
Astrophysics Research Centre, 5071 West Saanich Road, Victoria, BC V9E 2E7,
Canada 33institutetext: Kavli Institute for Cosmology, University of
Cambridge, Madingley Road, Cambridge, CB3 0HA, UK 44institutetext: Cavendish
Laboratory - Astrophysics Group, University of Cambridge, 19 JJ Thomson
Avenue, Cambridge, CB3 0HE, UK 55institutetext: Department of Physics,
University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK
66institutetext: Sorbonne Université, CNRS, UMR 7095, Institut d’Astrophysique
de Paris, 98 bis bd Arago, 75014 Paris, France 77institutetext: Scuola Normale
Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy 88institutetext: Cosmic
Dawn Center (DAWN), Copenhagen, Denmark Niels Bohr Institute, University of
Copenhagen, Jagtvej 128, DK-2200, Copenhagen, Denmark 99institutetext: Niels
Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200, Copenhagen,
Denmarkk 1010institutetext: Department of Physics and Astronomy, University
College London, Gower Street, London WC1E 6BT, UK 1111institutetext: European
Space Agency, c/o STScI, 3700 San Martin Drive, Baltimore, MD 21218, USA
1212institutetext: European Space Agency (ESA), European Space Astronomy
Centre (ESAC), Camino Bajo del Castillo s/n, 28692 Villanueva de la Cañada,
Madrid, Spain 1313institutetext: INAF - Osservatorio Astrofisco di Arcetri,
largo E. Fermi 5, 50127 Firenze, Italy 1414institutetext: European Southern
Observatory, Karl-Schwarzschild-Straße 2, 85748, Garching, Germany
1515institutetext: EUMETSAT - Fertigung für Verteidigung und Raumfahrt,
Darmstadt, Germany 1616institutetext: AURA for European Space Agency, Space
Telescope Science Institute, 3700 San Martin Drive. Baltimore, MD, 21210

# GA-NIFS: The ultra-dense, interacting environment of a dual AGN at $z\sim
3.3$ revealed by JWST/NIRSpec IFS

M. Perna e-mail<EMAIL_ADDRESS>The ultra-dense, interacting
environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS:
The ultra-dense, interacting environment of a dual AGN at $z\sim 3.3$ revealed
by JWST/NIRSpec IFS S. Arribas GA-NIFS: The ultra-dense, interacting
environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS:
The ultra-dense, interacting environment of a dual AGN at $z\sim 3.3$ revealed
by JWST/NIRSpec IFS M. Marshall GA-NIFS: The ultra-dense, interacting
environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS:
The ultra-dense, interacting environment of a dual AGN at $z\sim 3.3$ revealed
by JWST/NIRSpec IFS F. D’Eugenio GA-NIFS: The ultra-dense, interacting
environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS:
The ultra-dense, interacting environment of a dual AGN at $z\sim 3.3$ revealed
by JWST/NIRSpec IFSGA-NIFS: The ultra-dense, interacting environment of a dual
AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS: The ultra-dense,
interacting environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec
IFS H. Übler GA-NIFS: The ultra-dense, interacting environment of a dual AGN
at $z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS: The ultra-dense,
interacting environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec
IFSGA-NIFS: The ultra-dense, interacting environment of a dual AGN at $z\sim
3.3$ revealed by JWST/NIRSpec IFSGA-NIFS: The ultra-dense, interacting
environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFS A.
Bunker GA-NIFS: The ultra-dense, interacting environment of a dual AGN at
$z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS: The ultra-dense, interacting
environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFS S.
Charlot GA-NIFS: The ultra-dense, interacting environment of a dual AGN at
$z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS: The ultra-dense, interacting
environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFS S.
Carniani GA-NIFS: The ultra-dense, interacting environment of a dual AGN at
$z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS: The ultra-dense, interacting
environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFS P.
Jakobsen GA-NIFS: The ultra-dense, interacting environment of a dual AGN at
$z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS: The ultra-dense, interacting
environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS:
The ultra-dense, interacting environment of a dual AGN at $z\sim 3.3$ revealed
by JWST/NIRSpec IFSGA-NIFS: The ultra-dense, interacting environment of a dual
AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFS R. Maiolino GA-NIFS: The
ultra-dense, interacting environment of a dual AGN at $z\sim 3.3$ revealed by
JWST/NIRSpec IFSGA-NIFS: The ultra-dense, interacting environment of a dual
AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS: The ultra-dense,
interacting environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec
IFSGA-NIFS: The ultra-dense, interacting environment of a dual AGN at $z\sim
3.3$ revealed by JWST/NIRSpec IFSGA-NIFS: The ultra-dense, interacting
environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS:
The ultra-dense, interacting environment of a dual AGN at $z\sim 3.3$ revealed
by JWST/NIRSpec IFS B. Rodríguez Del Pino GA-NIFS: The ultra-dense,
interacting environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec
IFSGA-NIFS: The ultra-dense, interacting environment of a dual AGN at $z\sim
3.3$ revealed by JWST/NIRSpec IFS C. J. Willott GA-NIFS: The ultra-dense,
interacting environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec
IFSGA-NIFS: The ultra-dense, interacting environment of a dual AGN at $z\sim
3.3$ revealed by JWST/NIRSpec IFS T. Böker GA-NIFS: The ultra-dense,
interacting environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec
IFSGA-NIFS: The ultra-dense, interacting environment of a dual AGN at $z\sim
3.3$ revealed by JWST/NIRSpec IFS C. Circosta GA-NIFS: The ultra-dense,
interacting environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec
IFSGA-NIFS: The ultra-dense, interacting environment of a dual AGN at $z\sim
3.3$ revealed by JWST/NIRSpec IFSGA-NIFS: The ultra-dense, interacting
environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS:
The ultra-dense, interacting environment of a dual AGN at $z\sim 3.3$ revealed
by JWST/NIRSpec IFS G. Cresci GA-NIFS: The ultra-dense, interacting
environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS:
The ultra-dense, interacting environment of a dual AGN at $z\sim 3.3$ revealed
by JWST/NIRSpec IFS M. Curti GA-NIFS: The ultra-dense, interacting environment
of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS: The ultra-
dense, interacting environment of a dual AGN at $z\sim 3.3$ revealed by
JWST/NIRSpec IFS B. Husemann GA-NIFS: The ultra-dense, interacting environment
of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS: The ultra-
dense, interacting environment of a dual AGN at $z\sim 3.3$ revealed by
JWST/NIRSpec IFS N. Kumari GA-NIFS: The ultra-dense, interacting environment
of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS: The ultra-
dense, interacting environment of a dual AGN at $z\sim 3.3$ revealed by
JWST/NIRSpec IFS I. Lamperti GA-NIFS: The ultra-dense, interacting environment
of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS: The ultra-
dense, interacting environment of a dual AGN at $z\sim 3.3$ revealed by
JWST/NIRSpec IFS P. G. Pérez-González GA-NIFS: The ultra-dense, interacting
environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS:
The ultra-dense, interacting environment of a dual AGN at $z\sim 3.3$ revealed
by JWST/NIRSpec IFS J. Scholtz GA-NIFS: The ultra-dense, interacting
environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS:
The ultra-dense, interacting environment of a dual AGN at $z\sim 3.3$ revealed
by JWST/NIRSpec IFSGA-NIFS: The ultra-dense, interacting environment of a dual
AGN at $z\sim 3.3$ revealed by JWST/NIRSpec IFSGA-NIFS: The ultra-dense,
interacting environment of a dual AGN at $z\sim 3.3$ revealed by JWST/NIRSpec
IFS

(Received September 15, 1996; accepted March 16, 1997)

###### Abstract

Context. LBQS 0302$-$0019 is a blue quasar (QSO) at $z\sim 3.3$ that hosts
powerful outflows and resides in a complex environment consisting of an
obscured active galactic nucleus (AGN) candidate and multiple companions, all
within 30 kpc in projection.

Aims. We aim to characterise this complex system using JWST NIRSpec Integral
Field Spectrograph (IFS) observations obtained as part of the NIRSpec IFS GTO
programme ‘Galaxy Assembly with NIRSpec IFS’ (GA-NIFS); these data cover the
QSO rest-frame optical emission lines with a spatial resolution of $\sim
0.1\arcsec$ and a sampling of 0.05″($\sim 380$ pc) over a contiguous sky area
of $\sim 3\arcsec\times 3\arcsec$ ($23\times 23$ kpc2).

Methods. We developed a procedure to correct for the spurious oscillations (or
‘wiggles’) in NIRSpec single-spaxel spectra caused by the spatial under-
sampling of the point spread function. We performed a QSO–host decomposition
with the QDeblend3D tools. We used multi-component kinematic decomposition of
the optical emission line profiles to infer the physical properties of the
emitting gas in the QSO environment.

Results. The QSO–host decomposition allows us to identify both a low- and a
high-velocity component. The former possibly traces a warm rotating disk with
a dynamical mass $M_{dyn}\sim 10^{11}$ M⊙ and a rotation-to-random motion
ratio $v_{\rm{rot}}/\sigma_{0}\sim 2$. The other kinematic component traces a
spatially unresolved ionised outflow with a velocity of $\sim 1000$ km s-1 and
an outflow mass rate of $\sim 10^{4}$ M⊙ yr-1. We clearly detect eight
companion objects close to LBQS 0302$-$0019\. For two of them, we detect a
regular velocity field that likely traces rotating gas, and we infer
individual dynamical masses of $\approx 10^{10}$ M⊙. Another companion shows
evidence of gravitational interaction with the QSO host. Optical line ratios
confirm the presence of a second, obscured AGN $\sim 20$ kpc from the primary
QSO; the dual AGN dominates the ionisation state of the gas in the entire
NIRSpec field of view.

Conclusions. This work has unveiled in unprecedented detail the complex
environment of LBQS 0302$-$0019, which includes its host galaxy, a close
obscured AGN, and nine interacting companions (five of which were previously
unknown), all within 30 kpc of the QSO. Our results support a scenario where
mergers can trigger dual AGN and can be important drivers of rapid early
supermassive black hole growth.

###### Key Words.:

quasars: supermassive black holes – quasars: emission lines – Galaxies: high-
redshift – Galaxies: interactions – Galaxies: active – ISM: jets and outflows

Figure 1: Original (top) and PSF-subtracted (bottom) images of the QSO LBQS
0302$-$0019 and its close neighbouring galaxies as observed from space- and
ground-based telescopes. In the left panels, we show the HST WFC3 near-
infrared images from Husemann et al. (2021), with contours in the bottom panel
showing the Jil1-to-4 galaxies discovered by Husemann et al. The middle panels
present the MUSE Ly$\alpha$ emission before (top) and after (bottom) the QSO
PSF subtraction, from Husemann et al. (2018b), with contours from the PSF-
subtracted HST image. The right panels show the [O iii]$\lambda 5007$ emission
from JWST/NIRSpec observations; see Sect. 5.3 for details on the QSO PSF
subtraction. North is up, and east is to the left.

##  1 Introduction

The James Webb Space Telescope (JWST) promises to reveal a new view of galaxy
formation in the early Universe. Thanks to its unprecedented sensitivity and
spectroscopic capability in the near- and mid-infrared wavelengths, the rest-
frame optical nebular emission lines (e.g. H$\beta$, [O
iii]$\lambda\lambda$4959,5007, H$\alpha$, and [N ii]$\lambda\lambda$6548,6583)
of star-forming galaxies and active galactic nuclei (AGN) can, for the very
first time, be directly detected and resolved across early cosmic epochs, from
cosmic noon ($z\sim 2-3$; e.g. Förster Schreiber & Wuyts 2020) to the epoch of
re-ionisation ($z\gtrsim 7$; e.g. Robertson et al. 2022; Curtis-Lake et al.
2022). Early Release Observations and Cycle 1 General Observer and Guaranteed
Time Observations (GTO) programme results have clearly demonstrated the power
of JWST’s spectroscopic observations (e.g. Brinchmann 2022; Bunker et al.
2023; Cameron et al. 2023; Cresci et al. 2023; Curti et al. 2023; Kocevski et
al. 2023; Tacchella et al. 2022; Vayner et al. 2023), promising many exciting
discoveries over the coming years.

All cosmological models of hierarchical structure formation predict the
existence of multiple supermassive black holes (SMBHs) inside many galaxies,
consequences of previous merging events (Hopkins et al. 2007; Colpi 2014;
Volonteri et al. 2021). These events can be revealed by the detection of dual
AGN separated by up to a few kiloparsecs. The observational search for close
dual quasars (QSOs) at $1<z<3$ (i.e. at the peak of QSO activity) is
particularly important for constraining the merger process in cosmological
models because the effects of mergers are believed to be the most significant
in the high-luminosity, close-separation regime (e.g. Hopkins et al. 2008; Van
Wassenhove et al. 2012). Unfortunately, only very few dual AGN have been
confirmed observationally at such high $z$ (e.g. Chen et al. 2022a, b; Lemon
et al. 2022; Mannucci et al. 2022); whether these systems are intrinsically
rare or are simply undiscovered is not yet known. The study of the few dual
AGN known so far at high $z$ is therefore of paramount importance for testing
the predictions of the cosmological models in these early epochs of the
Universe. In this paper we use data from the JWST/NIRSpec Integral Field
Spectrograph (IFS; Jakobsen et al., 2022; Böker et al., 2022) of the optically
luminous QSO LBQS 0302$-$0019, one of the rare QSOs at high $z$ with a close
AGN (Husemann et al. 2018b).

The QSO LBQS 0302$-$0019 (RA $3^{\rm h}4^{\rm m}49.93^{\rm s}$, Dec.
$-0^{\circ}8^{\prime}13.10\arcsec$, J2000) at $z\sim 3.3$ has been intensively
targeted for studies of the intergalactic medium along our line of sight
(LOS). It is one of the rare ultraviolet-transparent luminous QSOs that allows
the He ii Ly$\alpha$ absorption of the intergalactic medium to be investigated
in detail: Worseck et al. (2021) inferred for LBQS 0302$-$0019 a large
proximity zone, 13.2 Mpc, caused by the enhanced ionising photon flux around
the QSO (e.g. Jakobsen et al. 1994), which implies a long active phase of more
than $11$ Myr for this QSO.

Analysing archival observations from the Multi Unit Spectroscopic Explorer
(MUSE; Bacon et al. 2010) on the Very Large Telescope (VLT), Husemann et al.
(2018b) report the detection of a Ly$\alpha$ nebula surrounding LBQS
0302$-$0019 out to tens of kiloparsecs that is associated with various high
ionisation lines. In particular, these authors report the serendipitous
discovery of an obscured AGN – dubbed Jil (Klingon for neighbour) – about 20
kpc from the QSO, inferred from Ly$\alpha$, C iv$\lambda$1549, He
ii$\lambda$1640, and C iii]$\lambda$1909 ultraviolet emission-line
diagnostics. The He ii line luminosity, L(He ii) $\sim 1.7\times 10^{42}$ erg
s-1, was inconsistent with being induced by LBQS 0302$-$0019 given the
compact, point-like spatial distribution of this line emission and its
corresponding small cross-section. The He ii luminosity can more easily be
explained by the presence of an AGN of about 1/500–1/1000 the luminosity of
LBQS 0302$-$0019 (corresponding to a bolometric luminosity of
L${}_{\rm{AGN}}\sim 10^{45}$ erg s-1), if located within the compact region
emitting He ii.

Follow-up ground-based Ks-band imaging and near-infrared spectroscopy are
presented in Husemann et al. (2018a), who successfully detected Jil’s host
galaxy emission, with an estimated stellar mass of $\sim 10^{11}$ M⊙, and the
optical [O iii]$\lambda 5007$ line ([O iii] hereinafter), with $L($[O
iii]$)\sim 2.5\times 10^{42}$ erg/s. However, no other rest-frame optical
lines were detected. Finally, Husemann et al. (2021) present Hubble Space
Telescope (HST) Wide-Field Camera 3 (WFC3) near-infrared imaging of the QSO,
revealing the presence of close multiple companion objects: emission from Jil
was resolved into two sources separated by $\sim 1\arcsec$ ($\sim 8$ kpc),
Jil1 and Jil2, while two additional sources were dubbed Jil3 and Jil4. They
also constrained stellar ages and masses for the two most prominent
companions, Jil1 with $t_{*}=252_{-109}^{+222}$ Myr and log(M∗/M⊙) =
$11.2_{-0.1}^{+0.3}$, and Jil2, associated with the compact He ii emission,
with $t_{*}=19_{-14}^{+74}$ Myr and log(M∗/M⊙) = $9.4_{-0.4}^{+0.9}$. These
early near-infrared (HST) and optical (MUSE) observations are presented in
Fig. 1 to display the complex environment of LBQS 0302$-$0019.

LBQS 0302$-$0019 also hosts a powerful outflow: Shen (2016), after analysing
near-infrared slit spectroscopy, reported the presence of an ionised outflow
traced by [O iii], with a velocity of $\gtrsim 1000$ km s-1. A velocity offset
of C iv relative to the centroid of the H$\beta$ broad line region (BLR) and
[O iii] narrow line region (NLR) of 400-600 km s-1 is also reported by Coatman
et al. (2017) and Zuo et al. (2020); such a significant displacement of the C
iv to the blue suggests the presence of strong nuclear outflows in the BLR of
LBQS 0302$-$0019 (see also e.g. Vietri et al. 2020).

In this manuscript we present the JWST/NIRSpec IFS observations of LBQS
0302$-$0019 to study the rest-frame optical lines and characterise its
intergalactic and interstellar medium. NIRSpec data enable us to shed light on
the gravitational interaction between the Jil sources and the QSO host galaxy,
as well as the possible accretion onto the QSO host through the circumgalactic
medium and the ejection of material through powerful outflows. The paper is
outlined as follows. In Sect. 2 we describe the JWST NIRSpec observations, and
our data reduction is outlined in Sect. 3. Detailed data analysis of the
integrated QSO spectrum and the spatially resolved spectroscopic analysis are
reported in Sects. 4 and 5, respectively. Section 5 also presents the new
procedure developed to model and subtract the wiggle artefacts in NIRSpec IFS
cubes. Finally, we present a discussion of our results in Sect. 6, before
concluding with a summary of our findings in Section 7.

Throughout, we adopt a Chabrier (2003) initial mass function
($0.1-100\leavevmode\nobreak\ M_{\odot}$) and a flat $\Lambda$ cold dark
matter cosmology with $H_{0}=70$ km s-1 Mpc-1, $\Omega_{\Lambda}=0.7$, and
$\Omega_{m}=0.3$. In the analysis we use vacuum wavelengths, but when
referring to emission lines we quote their rest-frame air wavelengths if not
specified otherwise.

##  2 Observations

LBQS 0302$-$0019 was observed on August 8, 2022, as part of the NIRSpec IFS
GTO programme ‘Galaxy Assembly with NIRSpec IFS’´´ (GA-NIFS) under programme
#1220 (PI: N. Luetzgendorf). The project is based on the use of the NIRSpec’s
IFS mode, which provides spatially resolved spectroscopy over a contiguous
3.1${}^{\prime\prime}\times$ 3.2′′ sky area, with a sampling of 0.1′′/spaxel
and a spatial resolution from $\sim 0.04$″(at $\sim 1\mu$m) to $\sim 0.15$″
(at $\sim 5\mu$m; see Böker et al. 2022; Rigby et al. 2022). 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 $\mu$m. The observations were taken with the
IRS2RAPID readout pattern with 60 groups, using a 4-point medium cycling
dither pattern, resulting in a total exposure time of 3560 seconds.

##  3 Data reduction

The raw data were reduced with the JWST calibration pipeline version 1.8.2,
using the context file jwst_1041.pmap. All of the individual raw images were
first processed for detector-level corrections using the Detector1Pipeline
module of the pipeline (Stage1 hereinafter). Then, the individual products
(count-rate images) were calibrated through Calwebb_spec2 (Stage2
hereinafter), where wcs-correction, flat-fielding, and the flux-calibrations
are applied to convert the data from units of count-rate to flux density. The
individual Stage2 images were then resampled and co-added onto a final data
cube through the Calwebb_spec3 processing (Stage3 hereinafter). A number of
additional steps (and corrections in the pipeline code) were applied to
improve the data reduction quality; different configurations were also used to
obtain additional data products and test the pipeline robustness (e.g. of flux
and spatial resolution recovery). In particular:

Figure 2: NIRSpec spectra obtained with internal (orange curve) and external
(green) flux calibration and integrated over a region of $r=1.5$″. The two
NIRSpec spectra are extracted from the drizzle cubes, with 0.05″spaxels.
Vertical lines indicate the main emission line features detected in the
NIRSpec spectrum. The inset shows the same NIRSpec spectra compared with the
Magellan/FIRE (magenta) and the SDSS (purple) spectra, rescaled by a factor of
$1.6$ to match the NIRSpec spectra in the vicinity of the H$\beta$ and [O iii]
lines.

  * •

In order to correct for the artefacts known as a ‘snowballs’, caused by large
cosmic ray impacts, we applied the snowball flagging for the jump during Stage
1. Sometimes this step incorrectly flags elongated streaks (due to cosmic ray
impacts) as snowballs. Even though these streaks affect only a narrow region
of the detector, the algorithm flags an entire circle containing the streak.
This results in extended, circular regions with signal over-subtraction in the
final count-rate images. To address this issue, we patched the pipeline to fit
ellipses to all flagged regions consisting of five or more adjacent pixels;
regions with best-fit ellipses having axis ratio smaller than 0.1 are removed
from the list of snowballs.

  * •

The individual count-rate frames were further processed at the end of Stage 1,
to correct for different zero levels in the dithered frames: for each image,
we subtracted the median value (computed considering the entire image) to get
a base level consistent with zero counts per second. This step is particularly
important for the very first frame obtained for LBQS 0302$-$0019, showing
(unrealistic) negative ramps in the raw (level 1b) data, and resulting in
negative counts at the end of Stage 1.

  * •

We further processed these count-rate images to subtract the $1/f$ noise (e.g.
Kashino et al. 2022). This correlated vertical noise is modelled in each
column (i.e. along the spatial axis) with a low-order polynomial function,
after removing all bright pixels (e.g. associated with the observed target)
with a $\sigma$-clipping algorithm. The modelled $1/f$ noise is then
subtracted before proceeding with Stage 2 of the pipeline.

  * •

The flux calibration was performed using two different approaches: the first
uses the photom step of Stage 2, and the second takes advantage of the
commissioning observations of the standard star TYC 4433-1800-1 (PID 1128,
o009). In the latter case, the flux calibration is performed as a post-
processing correction: we reduced the star with the same pipeline version and
context file, and obtained the response curve of the instrument required to
convert count rates into flux densities. Hereinafter, we refer to the first
approach as internal flux calibration, and to the second as external flux
calibration.

  * •

The outlier_detection step of Stage 3 is required to identify and flag all
remaining cosmic rays and other artefacts left over from previous calibration
steps, resulting in a significant number of spikes in the reduced data.
Unfortunately, with the current version of the pipeline, this step cannot be
used, because it tends to identify too many false positives and seriously
compromises the data quality111At the time of this writing, the newest version
of the pipeline, v1.9.4, and the latest context file, jwst_1063.pmap  are
still affected by these issues.. We therefore decided to follow two different
approaches to remove the spikes: the first one uses an algorithm similar to
lacosmic (van Dokkum, 2001) to remove outliers in individual exposures (at the
end of Stage 2): because our sources are under-sampled in the spatial
direction, we calculated the derivative of the count-rate maps only along the
(approximate) dispersion direction. The derivative was then normalised by the
local flux (or by 3$\times$ the noise, whichever was highest) and we rejected
the 95th percentile of the resulting distribution (see D’Eugenio et al. 2023
for details). The second approach consists of a post-processing correction,
and is done applying a $\sigma$ clipping to exclude all spikes in the reduced
data cubes (at spaxel level).

  * •

Finally, we applied the cube_build step to produce two combined data cubes:
one with a spaxel size of 0.1″, obtained with the emsm weighting (with higher
signal-to-noise at spaxel level), and a second with a spaxel size of 0.05″,
obtained with the drizzle weighting; the latter has a higher spatial
resolution but is more affected by point spread function (PSF) effects (see
Sect. 5.1). We manually rescaled the drizzle cubes by a factor of
$(0.05\arcsec/0.1\arcsec)^{2}$ to ensure the flux conservation.

We patched the cube_build script, fixing a bug affecting the drizzle algorithm
as implemented in the version 1.9.0222cube_build code changes in
https://github.com/spacetelescope/jwst/pull/7306; we also patched the photom
script, applying the corrections implemented in the same version333photom code
changes in https://github.com/spacetelescope/jwst/pull/7319, which allows more
reasonable flux densities to be inferred (i.e. a factor of $\sim 100$ smaller
with respect to those obtained with standard pipeline 1.8.2).

###  3.1 Astrometric registration

We obtained a bona fide astrometric registration matching the QSO nucleus
position with that in the HST image shown in Fig. 1, that is, applying a
correction of $\Delta$RA $=-0.492$″and $\Delta$DEC $=-0.062$″. This offset is
due to an error in the reference files responsible for the coordinate
transformation, then partially solved with the release of the context file
jwst_1063.pmap444 jwst_1063.pmap corrects for a $\sim 4$ pixels systematic
offset associated with the coordinate transformation between the ‘OTEIP’ and
the world systems, but not for a smaller offset ($\sim$0.2–0.4 pixels) between
the ‘GWA’ and the ‘virtual slit’ frame (see Dorner et al. 2016)..

Figure 3: Multi-component simultaneous best-fit results for the continuum-
subtracted spectrum of LBQS 0302$-$0019, around the H$\beta$-[O iii] (left)
and H$\alpha$-[N ii] regions (right; integrated over a circular region with
$r=0.5\arcsec$). The blue curve represents the rest-frame NIRSpec spectrum,
and the red curve indicates the best fit, with individual kinematic components
shown with different colours (as labelled in the right panel). For the outflow
(1) and (2) components, we also show the contribution from the only [N ii]
lines with dashed lines, as the grey and black curves do not allow a clear
distinction between the H$\alpha$ and [N ii] line transitions. Vertical red
lines indicate the most prominent emission lines, as in Fig. 2. The lower
panels show the residual to the model fit, that is, the difference between the
observed spectrum and the model.

###  3.2 Recovery of the QSO flux

Figure 2 shows the integrated NIRSpec spectra of LBQS 0302$-$0019, obtained
from the drizzle cubes, reduced with the internal (orange curve) and external
(green) flux calibration. The NIRSpec spectra were extracted from a circular
aperture centred at the position of the QSO nucleus, with $r=1.5$″, hence
matching the Sloan Digital Sky Survey (SDSS) fibre radius (see below). These
spectra are compared in the inset with the near-infrared Magellan/FIRE
spectrum (magenta, from Shen 2016) and the SDSS spectrum (in purple); the
latter is rescaled by a factor of 1.6 to match the fluxes in the vicinity of
the H$\beta$ and [O iii] lines.

The agreement between NIRSpec and the spectra from other facilities is
remarkable, and the small differences can be explained by taking flux
calibration uncertainties into account. The small mismatch between the two
integrated NIRSpec spectra (obtained with external and internal flux
calibrations) is of the order of $\sim 2-3\%$, well within the nominal
uncertainties of the JWST calibration pipeline (Böker et al. 2023). Being in
the very early stages of pipeline development, we avoided investigating the
discovered discrepancies further; nevertheless, we note that a larger mismatch
would be present without applying all corrections reported in Sect. 3 for the
internal calibration.555Systematic, wavelength-dependent discrepancies (up to
$\sim 25\%$) would also appear using context files older than jwst_1023.pmap,
because of placeholder flat field corrections used during the photom step of
the internal calibration.

All results described in this paper refer to the drizzle data cubes, which we
preferred over the emsm cubes as the former better preserve the NIRSpec
spatial resolution. Moreover, we preferred the internal to the external flux
calibration, as the former also allow corrections for the flat field. Finally,
we used the cubes obtained with the modified outlier detection method (see
above), although there are no major differences between these and those
corrected with a $\sigma$-clipping method.

##  4 Spectral analysis of the integrated LBQS 0302$-$0019spectrum

###  4.1 Spectral fit

We fit 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\lambda$4959,5007, [N ii]$\lambda\lambda$6548,83, and [S
ii]$\lambda\lambda$6716,31 doublets with a combination of Gaussian profiles,
applying a simultaneous fitting procedure (e.g. Perna et al. 2020), so that
all line features of a given kinematic component have the same velocity
centroid and full width at half maximum (FWHM). The modelling of the H$\alpha$
and H$\beta$ BLR emission requires the use of broken power-law components
(e.g. Nagao et al. 2006; Cresci et al. 2015; Trefoloni et al. 2023): they are
preferred over a combination of extremely broad Gaussian profiles because the
former tend to minimise the degeneracy between NLR and BLR emission. Finally,
we used the theoretical model templates of Kovacevic et al. (2010) to
reproduce the iron (Fe ii) emission in the wavelength region $4000-5500\AA$.
The final number of kinematic components used to model the spectra is derived
on the basis of the Bayesian information criterion (BIC; Schwarz 1978).

Figure 3 shows the best-fit model around the H$\beta$-[O iii] and H$\alpha$-[N
ii] regions. The BLR emission is fitted with a broken power law; iron emission
is fitted with the S and G group lines (Kovacevic et al. 2010). The [O iii]
doublet shows a narrow core, and prominent blue and red wings, and requires
three Gaussian components. To reduce the degeneracy between BLR and NLR, we
simultaneously fit four additional spectra extracted from circular regions
with radius of 0.2″(4 spaxels) and centred at different positions within a few
spaxels from the peak emission of the QSO: BLR profiles are tied, assuming
that these emission components originate from the same unresolved region,
while all other components are free to vary as originating from more extended
(and likely resolved) regions. The outcomes of this simultaneous fit (reported
in Fig. 15) are therefore used to fix the BLR parameters during the fit of the
integrated spectrum shown in Fig. 3.

We note that the integrated spectra reported in Fig. 15 show additional peaks
and/or inflection points in the H$\alpha$-[N ii] complex, due to the presence
of strong [N ii] emission line components; these nitrogen features are not
resolved in the integrated spectrum in Fig. 3, although they are still
definable from our fit decomposition. The absence of inflection points in Fig.
3 is likely due to the more prominent BLR emission, and the stronger
degeneracy between BLR and NLR kinematic components.

###  4.2 Systemic redshift

We derived the LBQS 0302$-$0019 redshift from the measured wavelength of the
narrow [O iii] emission in the integrated spectrum shown in Fig. 3:
$z=3.2870\pm 0.0003$, which is in agreement with Zuo et al. (2015, 2020) but
at odds with other redshift measurements from the literature. Husemann et al.
(2018b) reported values in the range $3.2882-3.2887$ (for different
ultraviolet lines); Coatman et al. (2019) reported $z=3.2856\pm 0.0002$ for
the [O iii], and $z=3.2868\pm 0.0012$ for the H$\beta$. All these previous
measurements are within $\approx\pm 100$ km s-1 of our zero velocity (assuming
$z=3.2870$). These small discrepancies are likely due to the presence of
powerful outflows, affecting all of the most prominent ultraviolet-to-optical
emission line profiles (Sect. 4.5).

###  4.3 Velocity offset between BLR and NLR

As shown in Fig. 3, the BLR emission line components are blueshifted with
respect to the [O iii] core component, by $480\pm 60$ km s-1. Relative
redshiftings (and blueshiftings) of the peaks of the broad Balmer line
emission are quite common in AGN (e.g. Gaskell 1983). Different explanations
for these offsets have been proposed: they could due to the orbital motion of
a SMBH binary (e.g. Ju et al. 2013), to recoiling SMBHs (e.g. Komossa et al.
2008), or to a perturbed accretion disk around a SMBH (e.g. Gaskell 2010). We
did not investigate these scenarios further as they go beyond the goals of
this paper; however, we note that each explanation is plausible given the
complex environment of LBQS 0302$-$0019.

###  4.4 Black hole mass

Assuming that the gas in the BLR is virialised, we calculated the central
black hole mass from the spectral properties of the H$\alpha$ and H$\beta$ BLR
region following the single-epoch calibrations from Dalla Bontà et al. (2020),

$\begin{split}M_{\rm{BH}}(H\beta)&=1.87\times 10^{6}\\\
&\times\leavevmode\nobreak\ \left(\frac{L_{\rm{H\beta}}}{10^{42}\
\rm{erg\leavevmode\nobreak\
s}^{-1}}\right)^{0.703}\left(\frac{\rm{\sigma}_{\rm{H\beta}}}{10^{3}\
\rm{km\leavevmode\nobreak\ s}^{-1}}\right)^{2.183}\rm{M}_{\odot}\end{split},$
(1)

with an intrinsic scatter of $\sim 0.3$ dex, and from Greene & Ho (2006):

$\begin{split}M_{\rm{BH}}(H\beta)&=3.6\times 10^{6}\\\
&\times\leavevmode\nobreak\ \left(\frac{L_{\rm{H\beta}}}{10^{42}\
\rm{erg\leavevmode\nobreak\
s}^{-1}}\right)^{0.56}\left(\frac{\rm{FWHM}_{\rm{H\beta}}}{10^{3}\
\rm{km\leavevmode\nobreak\ s}^{-1}}\right)^{2}\rm{M}_{\odot}\end{split}$ (2)

$\begin{split}M_{\rm{BH}}(H\alpha)&=2.0\times 10^{6}\\\ &\leavevmode\nobreak\
\times\leavevmode\nobreak\ \left(\frac{L_{\rm{H\alpha}}}{10^{42}\
\rm{erg\leavevmode\nobreak\
s}^{-1}}\right)^{0.55}\left(\frac{\rm{FWHM}_{\rm{H\alpha}}}{10^{3}\
\rm{km\leavevmode\nobreak\ s}^{-1}}\right)^{2.06}\rm{M}_{\odot},\end{split}$
(3)

with larger intrinsic scatters of $\sim 0.4$ dex. Aside from the small
differences in the intrinsic scatters in the chosen relations, we stress that
all single epoch relations reported in the literature have been inferred for
low-$z$ and low-luminosity AGN; as a result, significant extrapolations are
required for the measurement of the LBQS 0302$-$0019 black hole mass.

We find a H$\alpha/$H$\beta$ flux ratio of $3.88_{-0.12}^{+0.19}$ for the BLR
components. Taking as a reference the distribution of values of BLR Balmer
ratios obtained by Dong et al. (2008) for a large, homogeneous sample of
$\sim$ 500 low-$z$ Seyfert 1 and QSOs with minimal dust extinction effects,
H$\alpha/$H$\beta=3.06\pm 1.11$ (see also Baron et al. 2016), our Balmer
decrement measurement does not suggest significant extinction in the BLR of
LBQS 0302$-$0019\. Therefore, we did not perform any extinction correction for
the Balmer line luminosities required to compute the MBH.

The Balmer line luminosities and widths are measured from our best-fit BLR
profiles shown in Fig. 3 (i.e. the broken power-law components); we obtain
estimates of the black hole mass of the order of $\sim 2\times 10^{9}$ M⊙.
These values are broadly consistent with those previously reported in the
literature and are based on H$\beta$ and C iv BLR measurements (with the
latter being slightly larger, as commonly reported in the literature; e.g.
Coatman et al. 2017).

We calculated the bolometric luminosity of LBQS 0302$-$0019 following Dalla
Bontà et al. (2020), hence using the H$\beta$ BLR luminosity (their Eq. 25):
log(Lbol/ [erg s-1]) $\sim 47.2$ (consistent with the value inferred from the
continuum luminosity at $5100\AA$, and using the Netzer 2019 bolometric
correction, $\sim 47.3$). This bolometric luminosity is also consistent with
the one obtained starting from the intrinsic X-ray luminosity reported by
Nardini et al. (2019), and computed applying a bolometric correction
$k_{bol}=21$ (from Eq. 3 by Duras et al. 2020), log(Lbol/ [erg s-1]) $\sim
46.7$.

Using our black hole mass estimate from the H$\beta$ BLR (Eq. 1, which has
smaller scatter than Eqs. 2 or 3), we find an Eddington ratio of $\lambda_{\rm
Edd}=0.9\pm 0.1$. This value indicates that the accretion onto the central
black hole is close to the Eddington limit. All measurements so far inferred,
and the quantities required for their computation are reported in Table 1.

Table 1: Measurements of central black hole mass (with errors that include the intrinsic scatter of the single-epoch relations mentioned in the text), [O iii] luminosity (corrected for extinction), and outflow velocity from the integrated nuclear spectrum (see Sect. 4). Measurement | Value
---|---
$\log(L_{\rm H\alpha}/(\rm erg/s))$ | $45.45\pm 0.01$
FWHMHα [km/s] | $3650_{-60}^{+90}$
$\sigma_{\rm H\alpha}$ [km/s] | $3260_{-120}^{+40}$
$M_{\rm BH}({\rm H\alpha})_{GH06}[M_{\odot}]$ | $2.28^{+0.06}_{-0.09}\times 10^{9}$
$\log(L_{\rm H\beta}/(\rm erg/s))$ | $44.86\pm 0.02$
FWHMHβ [km/s] | $3560_{-120}^{+170}$
$\sigma_{\rm H\beta}$ [km/s] | $2800_{-70}^{+120}$
$M_{\rm BH}({\rm H\beta})_{GH06}[M_{\odot}]$ | $1.99^{+0.20}_{-0.16}\times 10^{9}$
$M_{\rm BH}({\rm H\beta})_{DB20,\ \sigma}[M_{\odot}]$ | $1.82_{-0.11}^{+0.19}\times 10^{9}$
$\log(L_{\rm bol}/(\rm erg/s))_{DB20}$ | $47.22\pm 0.01$
$\lambda_{\rm Edd}$ | $0.9\pm 0.1$
$\log(L_{[OIII]}/(\rm erg/s))$ | $45.71\pm 0.03$
[O iii] $W80$ [km/s] | $1080\pm 30$
[O iii] $W90$ [km/s] | $1600\pm 55$
[O iii] $V10$ [km/s] | $-760\pm 50$
Figure 4: Comparison of the BLR profiles of the C iv (from Shen 2016) and
H$\alpha$and H$\beta$ (from the integrated NIRSpec spectrum), in velocity
space. For the Balmer lines, we also report the best-fit BLR profiles. All
line profiles have been normalised to the flux of the BLR component in the
reddest parts, which are likely less affected by BLR and NLR outflows. The C
iv shows a significant excess in the blue part, at velocities of a few
thousand km s-1, which is not observed in the Balmer lines or in the [O iii]
line. This excess possibly indicates strong BLR winds.  Figure 5: Sinusoidal-
type patterns in single-spaxel spectra extracted from the drizzle data cube
with a spaxel size of 0.05″. Top panel: LBQS 0302$-$0019 spectrum integrated
over an aperture of $r=0.5$″(orange curve), in comparison with the spectrum of
the brightest spaxel (blue curve). Both spectra are normalised to their
maximum values, for visualisation purposes. The wiggles affecting the single-
spaxel spectrum are reported in grey and are obtained as the difference
between the blue and orange curves (after subtracting a low-order polynomial
function that takes the differences in the continuum level into account).
Bottom panel: Wiggles obtained from the eight pixels closest to the brightest
one. These wiggles strongly affect the shape of the continuum and, in
particular, the H$\beta$ profile and the wings of the [O iii] lines. See Fig.
16 for analogous effects in the emsm cube.

###  4.5 Outflow properties

Before analysing the [O iii] profile, we investigated the possible presence of
winds in the BLR of LBQS 0302$-$0019\. In Fig. 4 we compare the Balmer line
profiles with the C iv$\lambda 1549$ line (from ground-based observations;
Shen 2016). All profiles are normalised to the emission in the red wing,
unambiguously associated with the BLR emission for all transitions (see also
Fig. 1 in Zuo et al. 2020). Figure 4 highlights the presence of a blue excess
in the high-ionisation C iv line, up to a few thousand km s-1, likely due to
the presence of BLR winds. Such C iv outflows are commonly observed in
luminous high-z QSOs, and often associated with large-scale [O iii] outflows
(e.g. Coatman et al. 2019; Vietri et al. 2020).

From the best-fit model shown in Fig. 3, we inferred the [O iii] outflow
velocity, considering different tracers commonly used in the literature:
$V10\sim-760$ km s-1, the velocity at the 10th percentile of the overall
emission-line profile, $W80\sim 1080$ km s-1, defined as the line width
containing 80% of the emission line flux (obtained as the difference between
the velocities at 90th and 10th percentiles), and $W90\sim 1600$ km s-1,
containing 90% of the line flux (and obtained as the difference between the
velocities at 95th and 5th percentiles). All of these measurements are
consistent with previous values obtained for LBQS 0302$-$0019 from ground-
based observations (e.g. Villar Martín et al. 2020).

We anticipate here that the ionised outflow is not spatially resolved in our
NIRSpec observations; hence, outflow energetics, reported in Sect. 6.2, have
been derived on the basis of spatially integrated quantities.

Figure 6: Modelling of the wiggles in single-spaxel spectra. Top panel:
integrated LBQS 0302$-$0019spectrum (orange curve), single-spaxel spectrum
(blue), and wiggles (grey) as already reported in Fig. 5. The red curve
represents the best-fit model of the wiggles. Central panel: Single-spaxel
spectrum after the correction for the wiggles (dark blue), in comparison with
the integrated spectrum (orange); the grey curve represents the new residuals
with respect to the integrated spectrum. Bottom panel: Best parameter for the
frequency of the sinusoidal functions used to model the wiggles (blue points);
a low-order polynomial function fitting these points is also reported. All
panels display red shaded regions (associated with the QSO emission lines)
that are excluded during the fit.

##  5 Spatially resolved spectroscopy

###  5.1 Sinusoidal-type patterns in NIRSpec IFS

The spatial under-sampling in the NIRSpec IFS may result in apparent wiggles
in the single-spaxel spectra close to the position of bright point sources,
such as stars and QSOs. This effect is inherent to the cube building process,
and is more pronounced in data cubes with better spatial sampling (i.e. in
data cubes with spaxels of 0.05″, and constructed with the drizzle weighting
method). Further details about this effect, also known as ‘resampling noise’
can be found for instance in Smith et al. (2007) and Law et al. (2023). There
is currently no correction in the pipeline for this; large spatial extraction
regions are hence required to reduce the amplitude of the effect in extracted
1-D spectra. For isolated point sources, for which the extraction of spatially
resolved information is not possible, this effect is irrelevant, as when the
flux is integrated over a large aperture the wiggles disappear. However, there
are situations where a point source overlaps with extended emission, thus
requiring to disentangle the flux from both sources. This is the case, for
instance, in studies of QSO hosts and their close environment.

Figure 5 (top panel) displays the LBQS 0302$-$0019 spectrum integrated over an
aperture of $0.5\arcsec$ (in radius), in comparison with the spectrum of the
brightest spaxel extracted from the data cube constructed with the drizzle
weighting method (with spaxels of $0.05\arcsec$). The wiggles affecting the
single-spaxel spectrum are reported in the same panel with a grey curve, and
are obtained as the difference between the integrated and the single-spaxel
spectra (after subtracting a low-order polynomial function taking the
differences in the continuum levels into account). Similar sinusoidal-type
patterns are observed in all spaxels close to the brightest one, as shown in
the bottom panel of Fig. 5: they can affect a region as large as $r\sim
0.2-0.5$″.

The wiggles strongly limit the reconstruction and modelling of the target
spectrum at single-spaxel level. In particular, they affect the determination
of the continuum shape, and the modelling of permitted (e.g. Balmer) and
forbidden (e.g. [O iii]) emission lines. All of these components are required
to remove the signal from the nuclear point source (especially its PSF wings)
from the underlying extended emission (see e.g. Husemann et al. 2013; Marasco
et al. 2020).

These limitations also affect the single-spaxel spectra extracted from emsm
cubes with spaxels of 0.1″(see Fig. 16), although the amplitude of their
wiggles is $\approx 2-3$ times smaller than in the drizzle cubes. Moreover,
the use of emsm implies a decrease in spatial resolution, down to $\sim
0.2\arcsec$ (Vayner et al. 2023). In the next section we describe our approach
for modelling and subtracting these wiggles from NIRSpec data cubes; this
algorithm, written in python, is available for
download666https://github.com/micheleperna/JWST-NIRSpec_wiggles.

Figure 7: Reconstructed PSF from the spatially unresolved BLR emission. Left:
PSFs measured from the drizzle cube for the H$\alpha$ and H$\beta$ BLR
emission, respectively, as described in Sect. 5.3. The reconstructed H$\beta$
PSF is less extended than H$\alpha$, being at shorter wavelengths and
therefore associated with a smaller FWHM. Right: Visualisation of the BLR
subtraction in an individual spaxel at 0.14″north-east of the nucleus, using
the H$\beta$ (bottom) and H$\alpha$ (top) BLR template. The blue spectrum is
the original continuum-subtracted spectrum in the spaxel. The orange line is
the BLR model. Using the broad spectral windows marked in grey, the BLR model
is scaled to fit the original spectrum. The black curve shows the residual to
that fit, which is the BLR-subtracted spectrum.

###  5.2 Modelling of the wiggles

Figures 5 and 16 show sinusoidal-type patterns with relatively constant
amplitudes across the entire wavelength range, and significant variations for
the phase shift and the frequency within the $3\times 3$ innermost nuclear
spaxels. We note that the frequency changes smoothly along the whole
wavelength range, being almost constant in relatively narrow ranges; we took
advantage of this behaviour to model the wiggles.

As a first step, we fit the wiggles of the spectrum extracted from the
brightest spaxel, the one with highest signal-to-noise ratio (S/N). We used a
sinusoidal function to model the wiggles, $y(w)=A\sin(2\pi f_{w}w+\phi)+B$,
where A is the amplitude, $f_{w}$ is the frequency in $1/\mu$m, $w$ is the
wavelength, $\phi$ is the phase shift, and B is the continuum level; we
repeated the process in small portions of the wavelength range ($\sim
0.1\mu$m) as many times as necessary to cover the entire spectrum. The
combination of all best-fit sinusoidal functions is shown in the top panel of
Fig. 6 (red curve). The high spectral resolution and the small number of
parameters to fit the wiggles allow us to get a good representation of the
wiggles across the entire wavelength range, after masking the channels
associated with the most prominent emission lines and the gap between
detectors.

In the central panel of Fig. 6, we compare the integrated spectrum (orange)
with the corrected one (dark blue), obtained after subtracting the best-fit
model for the wiggles. The new residuals with respect to the integrated
spectrum are significantly smaller than the original ones (reported in grey in
the top panel).

By modelling the wiggles, we discover that the wiggle frequency, $f_{w}$,
changes smoothly as a function of the wavelength, as shown in the bottom panel
of Fig. 6: $f_{w}\sim 40\ \mu$m-1 at shortest and longest wavelengths, and
$f_{w}\sim 5\ \mu$m-1 in the central part of the spectrum. This $f_{w}$ trend
is common to all single-spaxel spectra around the QSO peak, and can be used to
better constrain the shape of the wiggles even for lower S/N spectra, or in
masked regions (associated with strong emission lines, and the gap between the
two detectors). As a final step, therefore, we fit all neighbouring spaxels
using the inferred $f_{w}$ as a prior for the modellisation of the wiggles.
Figures 17 and 18 show the same residuals presented in Figs. 5 and 16, but
after the correction described above. In Appendix B we also present some
caveats of our procedure.

We stress here that the wiggles behave similarly in all the data cubes of
bright point-like sources analysed so far within the GTO programme; the
procedure we described above is perfectly capable of modelling and correcting
for them. As an example, we report in Fig. 19 the wiggles modelling for
another target from our GTO programme, VDES J0020-3656, a QSO at $z=6.86$
observed with NIRSpec IFS with the grating–filter pair G395H–F290LP and
presented in Marshall et al. (2023). We note that data cube of this target has
been obtained by combining two datasets, observed with different telescope
position angles (on October 1, 2022, with $62^{\circ}$ and on October 16,
2022, with $160^{\circ}$). Nevertheless, the wiggles are very similar to those
in LBQS 0302$-$0019, consistent with the fact that these artefacts are
inherent to the cube building process.

###  5.3 QSO subtraction

Having corrected for the wiggles at single-spaxel level, we proceeded with the
separation between the host and QSO emission, making use of the QDeblend3D
routines (Husemann et al., 2013, 2014), which is optimised to subtract the PSF
emission from NIRSpec IFS data.

QDeblend3D  considers the relative strength of the BLR lines in each spaxel to
map out the spatial PSF, as the BLR is spatially unresolved. Due to the
NIRSpec PSF dependence with wavelength, we performed the QSO subtraction
twice: one for the wavelength channels around the H$\alpha$ line, taking as a
reference the H$\alpha$ BLR emission, and one for those in the vicinity of
H$\beta$, taking as a reference the H$\beta$ broad wings.

A PSF subtraction was performed following the procedure described in detail in
Marshall et al. (2023), also illustrated in Fig. 7. Briefly, we used the
previously built model for the BLR (and iron) emission (Sec. 4.1 and Fig. 3)
as a template, rescaled in each spaxel to fit the BLR emission in broad
spectral windows covering the wings of the Balmer lines (see Fig. 7). These
broad spectral windows are free from any narrow and outflow component
contributions, to avoid any bias in the measurement of the BLR strength.
Finally, we subtracted this rescaled template from each spaxel spectrum and
generated a new BLR-subtracted data cube.

Figure 8: PSFs measured from the drizzle cube for the nuclear emission around
H$\alpha$ (left) and H$\beta$ (right) and including both the BLR and outflow
components (as described in Sect. 5.3). With respect to Fig. 7, these panels
better reproduce the 2D distribution of the unresolved emission.

A fractional map of the relative brightness of the spatially unresolved BLR,
that is, the 2D PSF, is shown in the left part of Fig. 7, for both H$\alpha$
and H$\beta$. We note that the described subtraction does not take the NLR
emission into account, which is similarly spread according to the PSF shape.
To take this further contribution into account, we performed a different QSO
subtraction, this time using i) the integrated nuclear spectrum as a template,
and ii) broader spectral windows at both sides of the Balmer lines, including
the emission from high-velocity gas associated with the outflow (which is
unresolved in LBQS 0302$-$0019; see Sect. 6). This new reconstructed PSF is
shown in Fig. 8, and better reproduces the 2D distribution of unresolved
emission (as the NLR outflow wings have higher S/N than the BLR wings). The
cubes obtained from the subtraction of this high-velocity components (from
both NLR and BLR) are not used in the analysis described in the next sections,
but have been used to generate the [O iii] map shown in Fig. 1.

###  5.4 Line fitting

To derive spatially resolved kinematic and physical properties of ionised gas,
we fit the spectra of individual spaxels using the prescriptions already
presented in Sect. 4. We applied the BIC selection to determine where a
multiple-Gaussian fit is required to statistically improve the best-fit model.
This choice allows us to use the more degenerate multiple-component fits only
where they are really needed. For the spatially resolved analysis, we used two
Gaussian components at maximum, as they are perfectly capable of reproducing
the line profile variations in the field of view (FOV); this limited number of
components is also required to reduce the degeneracy in the fit.

Figure 9 shows the LBQS 0302$-$0019 velocity diagram, with all kinematic
parameters of the Gaussian components required to fit the BLR-subtracted data
cube. There is a clear trend in the figure, with the highest FWHMs ($>500$ km
s-1) associated with significant blueshifts ($\Delta v<-100$ km s-1), as
usually observed in systems hosting AGN outflows (e.g. Woo et al. 2016; Perna
et al. 2022). The Gaussian components with smaller FWHMs have relatively small
offsets from the zero velocity (up to a few hundred km s-1). In this figure,
we use different colours to distinguish between different regions (targets) in
the FOV: while the LBQS 0302$-$0019 host (black points) is often associated
with extreme kinematic parameters, all other companions (see the labels) show
narrower profiles possibly associated with rotation. A detailed
characterisation of the individual kinematic systems is reported in the next
sections.

Figure 9: Velocity diagram for the individual Gaussian components used to
model the emission line profiles in the data cube. Different colours are used
to identify different targets in the NIRSpec FOV, as labelled. Figure 10:
From left to right: H$\alpha$ and [O iii] flux distributions, and Moment 1 and
Moment 2 maps of the QSO host galaxy, obtained from the narrow components of
our best-fit models. Both lines show evidence of rotating gas in the QSO host.

##  6 Results

###  6.1 QSO host disk

Figure 10 shows an overview of the flux distribution and kinematics of the
narrow component in the LBQS 0302$-$0019 host galaxy, as derived from our
modelling of the [O iii] line (top panels) and H$\alpha$ (bottom) in the BLR-
subtracted data cube. The flux distribution of the two lines is dominated by
the nuclear emission, which spreads according to the PSF (see Fig. 8),
although a few clumps towards the east and south-east as well as an extended
plume towards the north-east (in the [O iii] map) can be easily recognised.
All of these features, reasonably associated with different sources in the QSO
host environment, are discussed in the next section.

The velocity distribution, traced by the Moment 1, shows evidence for a
velocity gradient along the north-east–south-west direction, with a velocity
amplitude of $\sim\pm 120$ km s-1, possibly associated with a rotating disk.
The most significant deviations from this gradient are found in the external
regions, in correspondence with the clumps and the plume identified in the
flux distribution panel. We also note that the H$\alpha$ velocity field is
noisier than the [O iii] one, because of the BLR subtraction step, and the
degeneracy between H$\alpha$ and [N ii] lines.

The [O iii] and H$\alpha$ line widths, traced by the Moment 2 map, do not show
significant variations across the host. However, elevated dispersions in the
central region of the galaxy in both the H$\alpha$ and [O iii] maps might be
present.

As the H$\alpha$ maps are probably more affected by PSF artefacts and BLR-
subtraction, we decided to use the [O iii] line to model the gas kinematics
with 3D-Barolo (Di Teodoro & Fraternali 2015), following the procedure
described in Perna et al. (2022), to test whether the QSO host kinematics are
compatible with a rotation-supported system and to infer the host dynamical
mass. The main assumption of the 3D-Barolo model is that all the emitting
material of the galaxy is confined to a geometrically thin disk, and its
kinematics are dominated by pure rotational motion. The possible presence of
residual components associated with the outflow, as well as the presence of
additional kinematic components associated with close companions might affect
the modelling. Nevertheless, this model enables us to assess the presence of
such disks and to infer a simple kinematic classification through the standard
$v_{\rm{rot}}/\sigma_{0}$ ratio, where $v_{\rm{rot}}$ is the intrinsic maximum
rotation velocity (corrected for inclination, $v_{\rm{rot}}=v_{LOS}/sin(i)$)
and $\sigma_{0}$ is the intrinsic velocity dispersion of the rotating disk,
related to its thickness. In this work, we define $\sigma_{0}$ as the measured
line width in the outer parts of the galaxy, corrected for the instrumental
spectral resolution (e.g. Förster Schreiber et al. 2018). The 3D-Barolo best-
fit plots are shown in Fig. 22. From them, we infer an inclination $i=10\pm
7^{\circ}$, a $v_{\rm{rot}}=360\pm 80$ km s-1, and a $\sigma_{0}=190\pm 45$ km
s-1 (at $\sim 0.3\arcsec$, i.e. $\sim 2.4$ kpc from the nucleus, as the more
external regions are more affected by noise). The rotation-to-random motion
ratio $v_{\rm{rot}}/\sigma_{0}\approx 2$ indicates that this galaxy is
associated with a dynamically warm disk, consistent with $z\sim 2$ galaxies
presented in Förster Schreiber et al. (2018), with $v_{\rm{rot}}/\sigma_{0}$
spanning the range from 0.97 to 13 (with a median of 3.2), as inferred from
H$\alpha$ gas kinematics (see also e.g. Wisnioski et al. 2019).

The 3D-Barolo best-fit velocity maps also show significant residuals in the
receding part, at $\sim 0.15\arcsec$ south-west of the nucleus, with
velocities $\approx 100$ km s-1; they might be associated with a plume, or a
further companion on the LOS. This kinematic component might also be present
in the integrated spectrum in Fig. 3: the significant residuals in the red
part of the H$\alpha$ line, if due tp H$\alpha$ line, would correspond to
L(H$\alpha$) $\approx 10^{43}$ erg/s, consistent with the luminosity of other
Jil companions (see Table 2).

From the 3D-Barolo best fit, we also inferred a tentative estimate for the
dynamical mass, assuming that the source of the gravitational potential is
spherically distributed (following e.g. Perna et al. 2022): Mdyn = $(14\pm
6)\times 10^{10}$ M⊙, within a radius of $2.4\pm 0.6$ kpc (corrected for the
PSF, and containing 85% of the [O iii] total flux, as inferred from the QSO-
subtracted cube). Combining this measurement with the MBH derived in Sect. 4,
we obtained a M${}_{\rm{BH}}/$M${}_{\rm{dyn}}\approx 0.014$. This places the
LBQS 0302$-$0019 host galaxy slightly above the local black hole-–host mass
relation (Kormendy & Ho 2013), consistent with other high-$z$ QSOs reported in
the literature (see e.g. Marshall et al. 2023 and references therein).

We finally investigated the dominant ionisation source for the emitting gas
across the LBQS 0302$-$0019 host, using the classical ‘Baldwin, Phillips &
Terlevich’ (BPT) diagram (Baldwin et al. 1981). The distributions of the flux
ratio diagnostics are almost constant across the host galaxy extension, with
log([NII]/H$\alpha$) $=-0.3\pm 0.2$ and log([O iii]/H$\beta$) $=0.7\pm 0.2$.
These values place the LBQS 0302$-$0019 host in the AGN-dominated region of
the BPT diagram (Baldwin et al. 1981; Kewley et al. 2013; Fig. 14).

Figure 11: From left to right: H$\alpha$ and [O iii] flux distributions, and
Moment 1 and Moment 2 maps of the Jil companions, obtained from the total
profiles of our best-fit models. Both lines show evidence of rotating gas in
the north-east companions.

###  6.2 QSO outflow energetics

The outflow component used to model the QSO host is not spatially resolved,
and is therefore not reported in the figures. In this section we measure the
mass of the ionised outflow as inferred from the blueshifted outflow component
of H$\beta$. We used the equation

$\dot{M}_{\rm{out}}(H\beta)=8.6\times\frac{L_{41}(H\beta)\
v_{\rm{out}}}{n_{e}\ R_{\rm{out}}}M_{\odot}yr^{-1}$ (4)

from Cresci et al. (2015), where L41(H$\beta$) is the H$\beta$ luminosity
associated with the outflow component in units of $10^{41}$ erg s-1, $n_{e}$
is the electron density, $v_{\rm{out}}$ is the outflow velocity, and
$R_{\rm{out}}$ is the radius of the outflowing region in units of kiloparsecs.

In general, $n_{e}$ can be estimated from the [S ii] doublet ratio (e.g.
Osterbrock & Ferland 2006), using the high-velocity components of the [S ii]
lines. Unfortunately, these components are only barely detected in our
integrated spectra, and cannot be used to infer the outflow electron density.
We therefore conservatively considered an electron density of 1000 cm-3,
inferred from the study of large samples of AGN both at low redshift ($z<0.8$,
Perna et al. 2017) and at $0.6<z<2.7$ (Förster Schreiber et al. 2019). A
factor of $\sim 3$ higher mass rate would be obtained for instance using the
electron density measured in the outflowing gas of the QSO XID2028 at $z\sim
1.5$ (i.e. $360\pm 180$ cm-3), as measured from recent JWST/NIRSpec IFS
observations (Cresci et al. 2023).

Here we consider the L41(H$\beta$) to be the luminosity of the H$\beta$
outflow component as measured from our full integrated spectral fit described
in Sect. 4 and shown in Fig. 3, as the outflow is not resolved in our data
cube: log(L(H$\beta$)/[erg/s]) $=45.17_{-0.13}^{+0.09}$. The luminosity has
been corrected for the extinction considering the colour excess for the same
outflow component, inferred from the Balmer decrement and assuming a Milky Way
extinction law (Cardelli et al. 1989): $E(B-V)=0.58_{-0.15}^{+0.07}$.

The identification of the BLR component in LBQS 0302$-$0019 suggests that the
outflow could be primarily orientated towards us; this could also explain why
the ejected gas is not spatially resolved, regardless the exquisite NIRSpec
resolution ($\sim 800$ pc). Under this assumption, the observed velocity
offset of the outflow components with respect to the BLR systemic is close to
the true outflow velocity (e.g. Harrison et al. 2012); as the outflow
component in the integrated spectrum requires the use of two Gaussian
components, we decided to use as velocity offset the v50 inferred from the
total outflow profile. We therefore derive a $v_{\rm{out}}=930_{-110}^{+60}$
km s-1.

The last ingredient required for the computation of the mass rate is the
outflow extension; as this component is not spatially resolved in our NIRSpec
cube, we assumed that the outflow is propagating at constant velocity (e.g.
Brusa et al. 2015; Fiore et al. 2017), and that its dynamical time ($t_{d}$)
is equal to the AGN phase inferred by Worseck et al. (2021), $t_{d}>11$ Myr.
This is very close to the $t_{d}$ usually inferred from observations of
ionised outflows (e.g. Greene et al. 2012; Perna et al. 2015a). We therefore
estimate $R_{\rm{out}}=t_{d}\times v_{\rm{out}}\gtrsim 9$ kpc. This estimate
is compatible with the extension of ionised outflows observed in other QSOs at
high $z$, in the range $\approx 2-15$ kpc (Carniani et al. 2015; Kakkad et al.
2020; Cresci et al. 2023). Because of that, we considered the inferred lower
limit as an order of magnitude estimate for the outflow extension.

We therefore obtain an outflow mass rate $\dot{M}_{\rm{out}}(H\beta)\sim
10^{4}$ M⊙ yr-1. This value, although significantly larger than other mass
rates reported in the literature, is still consistent with the general
expectations inferred from the scaling relations presented, for instance, in
Fiore et al. (2017) and Fluetsch et al. (2019). The inferred value is $\sim 3$
times larger than the one obtained from the [O iii] gas, following Carniani et
al. (2015); similar discrepancies are often reported in the literature (e.g.
Carniani et al. 2015; Perna et al. 2015a, 2019; Marshall et al. 2023) and are
probably due to the ionisation structure of the [O iii] and H$\beta$ clouds in
the NLR of an AGN. The kinetic and momentum powers are
$\dot{E}_{\rm{out}}=1/2\dot{M}_{\rm{out}}v_{\rm{out}}^{2}\sim 4\times 10^{45}$
erg s-1 and $\dot{P}_{\rm{out}}=\dot{M}_{\rm{out}}v_{\rm{out}}\sim 8\times
10^{37}$ dyne, respectively. Hence, the kinetic power is $\sim 2$% of the
radiative luminosity of the AGN, while the momentum rate is in excess of $\sim
15$ times the radiative momentum flux (Lbol/c), consistent with the energetics
of other QSOs in the literature (see e.g. Perna et al. 2015b; Bischetti et al.
2017; Tozzi et al. 2021).

###  6.3 Further considerations of the outflow extension

We report here two further arguments to better justify the assumed outflow
extension ($\approx 9$ kpc). On the one hand, greater extensions would be at
odds with the fact that the outflow is unresolved in our data cube: high
collimation (with a half opening angle $\alpha_{out}$ of a few degrees) would
be required to explain the presence of a spatially unresolved ($\lesssim 0.8$
kpc, i.e. below the spatial resolution of our data) and highly extended
outflow ($>9$ kpc) along our LOS, at odds with the reconstructed geometry of
other outflows at lower redshifts (with $\alpha_{out}\approx 10-60^{\circ}$;
e.g. Müller-Sánchez et al. 2011; Meena et al. 2021; Cresci et al. 2023). On
the other hand, by assuming that the outflow has an extension $<0.8$ kpc, we
would obtain outflow energetics that are ten times higher (e.g.
$\dot{M}_{\rm{out}}\sim 10^{5}$ M⊙ yr-1). Both the scenarios are quite
unlikely. We therefore conclude that the measurements reported in the previous
section can represent rough estimates of the outflow energetics for LBQS
0302$-$0019.

###  6.4 QSO environment: The Jil objects

Figure 11 shows an overview of the flux distribution and kinematics of the
ionised gas in LBQS 0302$-$0019 companion sources, as derived from our
modelling of the [O iii] line (top panels) and H$\alpha$ (bottom). The flux
distribution shows multiple clumps in the north-east regions, as well as
plumes and irregular structures within $\approx 1\arcsec$ of the LBQS
0302$-$0019 nucleus. All of these sources have relative velocity shifts up to
a few hundred km s-1 with respect to the QSO systemic; this implies that they
are not artefacts induced by the nuclear PSF. The velocity distribution,
traced by the Moment 1 of the total fitted profiles, shows evidence for
gradients with velocity amplitudes of $\sim\pm 200$ km s-1. The velocity width
(Moment 2) in these companions is significantly smaller than the ones in the
QSO host.

In order to better identify all possible companions around LBQS 0302$-$0019,
in Fig. 12 we show a few narrow-band images for the best-fit [O iii] emission
line, with overlaid contours from the HST image (already reported in Fig. 1,
left): these narrow-band images clearly show several clumps at different
velocities. Some of them are associated with the Jil companions already
identified by Husemann et al. (2021): Jil1, Jil2, and Jil3. However, Jil1 is
barely detected in our data as it resides on the very edge of the NIRSpec FOV,
where the noise is higher and the data reduction generates unreliable spectral
features. We also note that NIRSpec [O iii] emission slightly differs from the
flux distribution in the near-infrared HST, the former being more extended and
clumpier; this also makes it difficult to separate the Jil sources. Additional
[O iii] clumps not detected in the HST image are here dubbed Jil5, Jil6, Jil7,
Jil8, and Jil9, following Husemann et al. Their integrated spectra are shown
in the left part of Fig. 12.

The emission line properties of each companion, inferred from spectroscopic
analysis, are reported in Table 2. Here we give a brief description of the
specific properties inferred from each companion.

Figure 12: Jil spectra and spatial distribution. Left panels: Jil companions’
spectra. Jil2, Jil3, and Jil8 spectra have been obtained by integrating over
the entire extension of the targets (with [O iii] detected above $8\sigma$).
For the remaining sources, we considered $3\times 3$ spaxel integration (Jil5,
Jil6, Jil7, and Jil9), or $5\times 5$ spaxels (Jil1). All profiles are
relatively narrow and redshifted with respect to the QSO systemic; this
ensures that we are not affected by PSF contamination. The continuum emission
is never detected; the continuum observed in the Jil1 spectrum is likely due
to data reduction artefacts. Right panels: [O iii] flux distributions,
obtained by integrating over different velocity channels. The top-left panel
has been obtained by integrating over a large velocity range, using the QSO
NLR$+$BLR-subtracted cube (Sect. 5.3); all other velocity channel maps have
been extracted from the best-fit [O iii] profiles. The top-left panel shows
the HST contours presented in Husemann et al. (2021); all remaining panels
show the Jil companions, as detected in NIRSpec. The dashed red lines identify
the regions from which the integrated Jil2, Jil3, and Jil8 spectra, reported
on the left part of the figure, have been extracted.

Table 2: Properties of the companion sources in the LBQS 0302$-$0019 environment.  target |  $z$ ($\Delta v$ [km s-1]) | log($\frac{\rm{L([OIII])}}{\rm{erg\ s}^{-1}}$) | log($\frac{\rm{L(H\alpha)}}{\rm{erg\ s}^{-1}}$) | E(B$-$V) | log([O iii] / H$\beta$) | log([N ii] / H$\alpha$) |  $W80$ (km s-1)
---|---|---|---|---|---|---|---
Jil1 |  $3.2898\pm 0.0003$ (195) | $41.5\pm 0.1$ | $<40.9$ | $-$ | $>0.72$ | $-$ | $390\pm 30$
Jil2 |  $3.2877\pm 0.0002$ (50) | $42.9\pm 0.1$ | $42.6\pm 0.1$ | $0.2\pm 0.1$ | $0.81\pm 0.05$ | $<-1.04$ | $360\pm 30$
Jil3 |  $3.2904\pm 0.0002$ (240) | $43.6\pm 0.2$ | $43.0\pm 0.2$ | $0.5\pm 0.2$ | $1.00\pm 0.04$ | $<-0.82$ | $470\pm 30$
Jil5 |  $3.2886\pm 0.0002$ (110) | $40.9\pm 0.1$ | $<40.5$ | $-$ | $>0.7$ | $-$ | $220_{-30}^{+140}$
Jil6 |  $3.2866\pm 0.0004$ (-30) | $41.2\pm 0.1$ | $40.5\pm 0.2$ | $-$ | $0.67_{-0.11}^{+0.37}$ | $0.11_{-0.37}^{0.25}$ | $165_{-55}^{+30}$
Jil7 |  $3.2836\pm 0.0004$ (-240) | $41.4\pm 0.3$ | $<40.6$ | $-$ | $0.52_{-0.04}^{+0.36}$ | $-$ | $580\pm 80$
Jil8 |  $3.2848\pm 0.0004$ (-150) | $42.1\pm 0.2$ | $41.3\pm 0.2$ | $-$ | $0.69_{-0.07}^{+0.14}$ | $<-0.31$ | $440_{-55}^{+30}$
Jil9 |  $3.2858\pm 0.0004$ (-80) | $41.0\pm 0.1$ | $40.4\pm 0.2$ | $-$ | $>0.7$ | $<-0.10$ | $170_{-5}^{+30}$

Notes: For each target, in the second column we report the redshift and the
velocity offset with respect to the LBQS 0302$-$0019 host galaxy. Integrated
[O iii] and H$\alpha$ luminosities have been corrected for extinction, when
E(B$-$V) could be estimated, assuming a Milky Way extinction law (Cardelli et
al. 1989). For targets with no H$\beta$ detection, we measured the log([O
iii]/H$\beta$) lower limit assuming that the H$\beta$ upper limit is three
times smaller than H$\alpha$. The non-parametric velocity $W80$ refers to the
[O iii] line profile.

Jil1 is detected in [O iii] at $\sim 3\sigma$. This is the only companion for
which we could detect continuum emission; however, it has to be considered a
spurious measure, because of its position at the edges of the FOV and the
known issues with the data reduction.

Jil2 is detected in [O iii], H$\alpha$, and H$\beta$, but not in [N ii] or [S
ii] lines. It shows a clumpy morphology, with an extension over $\sim 4$ kpc
in projection. A velocity gradient with amplitude of $\sim\pm 100$ km s-1 is
observed. To test whether Jil2 is compatible with a rotationally supported
system, we modelled the [O iii] line with 3D-Barolo, as done for the QSO host.
The Jil2 best-fit models are shown in Fig. 23. The significant residuals in
the maps are likely due to the clumpy morphology of this system, as well as
the superposition with Jil3. These arguments likely explain the measured
rotation-to-random motion ratios, $v_{\rm{rot}}/\sigma_{0}\approx 0.7$, and
question the presence of a rotating disk. Nevertheless, we infer a tentative
dynamical mass for this system, log(Mdyn) = $(8_{-6}^{+16})\times 10^{9}$ M⊙,
considering a circular velocity of 75 km s-1(corrected for an inclination
$i=80^{\circ}\pm 15^{\circ}$), as measured with 3D-Barolo, and a radius of 1
kpc (as order of magnitude size, given the clumpy morphology of this source).
This dynamical mass is $\approx 3\times$ the stellar mass estimate inferred by
Husemann et al. (2021) (from spectral energy distribution analysis), but still
consistent within the errors.

Figure 13: Jil5 companion spectrum (light blue), together with two additional
spectra extracted from the region between Jil5 and the QSO host galaxy
(labelled as Jil5a and Jil5b, at 7.4 and 4.8 kpc from the QSO nucleus,
respectively, also indicated in Fig. 12 with red crosses). The dark blue
spectrum has been extracted from a region at a distance of 4.8 kpc from the
nucleus (as for Jil5b) but covering the PSF wing extending towards the north.
In order to ease the visualisation, we added vertical offsets to the spectra.
This figure highlights a velocity gradient of a few hundred km s-1 across a
few kiloparsecs (see also Fig. 11), possibly indicating feeding processes or a
tidal tail due to the interaction between Jil5 and the QSO host galaxy.

The presence of an obscured AGN in Jil2, initially proposed by Husemann et al.
(2018b) on the basis of the presence of bright ultraviolet lines in the MUSE
cube, will be discussed in the next section. Here we briefly mention that the
measured log([O iii]/H$\beta$) $\sim 0.8$, slightly higher than the value
obtained for the NLR gas associated with LBQS 0302$-$0019, is compatible with
the presence of an AGN in this companion.

Jil3 is detected in [O iii], H$\alpha$, and H$\beta$, but not in [N ii] or [S
ii] lines. It shows an elongated morphology, with a clear and regular gradient
with amplitude $\pm 200$ km s-1 over $\sim 8$ kpc. This extension, and the
clumpy morphology, could suggest the presence of multiple systems;
nevertheless, we also provide a 3D-Barolo model for the [O iii] emission for
this source. The best-fit results are reported in Fig. 24. We infer a
tentative dynamical mass log(Mdyn) = $(1.3\pm 0.5)\times 10^{10}$ M⊙, assuming
a circular velocity of 70 km s-1(corrected for an inclination $i=80\pm 10$, as
measured with 3D-Barolo) and a radius of 1 kpc (as for Jil2). From the
integrated spectrum of Jil3, we measure a high log([O iii]/H$\beta$) $\sim 1$,
consistent with possible presence of an AGN in its vicinity (i.e. in Jil2).

No further characterisation can be obtained for Jil4, which falls outside the
NIRSpec FOV.

Jil5 is located at $\sim 10$ kpc north-east of the LBQS 0302$-$0019 nucleus,
and is connected with the QSO host by a filamentary structure showing a clear
velocity gradient (see Fig. 11). To highlight the presence of such a gradient,
in Fig. 13 we show the Jil5 spectrum in comparison with those extracted from
intermediate positions along this elongated structure (identified by red
crosses in the velocity channels at $\sim 200$ km s-1in Fig. 12). We can
therefore speculate that this companion is contributing to the feeding of the
QSO host. For this companion we measure a log([O iii]/H$\beta$) $>0.7$,
consistent with flux ratios measured in the QSO host, and hence likely ionised
by the QSO radiation.

Jil6 is located at $\sim 10$ kpc south-east of the QSO, and is detected in [O
iii] and H$\alpha$ (and in H$\beta$ and [N ii] at S/N $\sim 2-3$). Both log([O
iii]/H$\beta$) $\sim 0.7$ and log([N ii]/H$\alpha$) $\sim 0.1$ suggest a QSO
ionisation.

Jil7 is located at $\sim 10$ kpc south-east of the QSO, and is detected in [O
iii] and H$\beta$. It shows a prominent blue wing in the [O iii]
($V10\sim-550\pm 50$ km s-1), likely due to the superposition of different
kinematic components along the LOS, and relatively high line ratios (log([O
iii]/H$\beta$) $\sim 0.52$).

Jil8 is located at $\sim 8$ kpc east of the QSO nucleus, with an extension of
$\sim 2$ kpc. It is detected in [O iii], H$\alpha$, and H$\beta$. The broad
components in the emission lines are due to PSF artefacts. For this companion,
log([O iii]/H$\beta$) $\sim 0.7$ suggests a QSO ionisation.

Jil9 is located at $\sim 23$ kpc north-east of the blue QSO, and is detected
in [O iii] and H$\alpha$. It shows a velocity offset of $\sim 300$ km s-1 from
Jil3, and narrow line profiles ($W80\sim 170$ km s-1). In this case, we detect
a lower limit for log([O iii]/H$\beta$) $>0.7$, consistent with the presence
of high ionisation.

####  6.4.1 Dual QSO with 20 kpc separation

All flux ratios so far inferred for the Jil targets (and for the QSO host
galaxy) are reported in Fig. 14. These constraints locate almost all Jil
sources in the AGN regions of the BPT diagram; for the remaining sources not
included in the diagram, Jil1, Jil5, and Jil7, for which we cannot detect
H$\alpha$ or [N ii], we can likely assume physical conditions similar to those
in the other Jil companions, because of the similarly high [O iii]/H$\beta$
ratios.

The line ratio diagram also shows that Jil2 and Jil3 galaxies are associated
with very stringent upper limits for the log([N ii]/H$\alpha$), of the order
of $\lesssim-1$. This may indicate that they are metal-poor AGN or galaxies,
consistent with model predictions ($Z\lesssim 0.5Z_{\odot}$; e.g. Groves et
al. 2004; Baron & Netzer 2019; see e.g. the predicted ratios from Nakajima &
Maiolino 2022 reported in the figure), and the ultraviolet diagnostics
(Husemann et al. 2018b). On the other hand, the LBQS 0302$-$0019 host might be
associated with a higher metallicity ($Z\approx Z_{\odot}$; according to the
same grid models), because of the higher [N ii]/H$\alpha$.

The relative proximity of Jil5, Jil6, Jil7, and Jil8 to the QSO is a likely
explanation for the high [O iii]/H$\beta$ in such targets. On the other hand,
the high [O iii]/H$\beta$ in Jil1, Jil2, Jil3, and Jil9 can be explained by
the presence of an AGN in Jil2, inferred by Husemann et al. (2018b) on the
basis of ultraviolet diagnostics. In support of this scenario, we used the He
ii$\lambda 4686$ diagnostics (Shirazi & Brinchmann 2012; Nakajima & Maiolino
2022; Übler et al. 2023; Tozzi et al. 2023). Since the He ii$\lambda 4686$ is
undetected in NIRSpec, we used the ratio He ii$\lambda 1640$/He ii$\lambda
4686=7.2$, expected for recombination (Seaton 1978), to infer the He
ii$\lambda 4686$ flux in Jil2 (correcting for extinction). This gives for Jil2
a log(He ii$\lambda 4686/$H$\beta$) $=-0.22$, consistent with AGN ionisation
(see Fig. 7 in Übler et al. 2023). We stress that the detection of He
ii$\lambda 1640$ emission line in the surroundings of QSOs (i.e. at scales
$>10$ kpc) is not common: for instance, this line has been tentatively
detected (at $\sim 2\sigma$) by stacking MUSE data cubes of 27 bright QSOs at
$z=3-4.5$ by Fossati et al. (2021, to be compared with the $>10\sigma$
detection in Jil2).

For Jil2, we also report a $\sim 4\sigma$ detection of the [S ii] doublet, and
hence a log([S ii]/H$\alpha$) $=-0.52\pm 0.07$. This value places Jil2 in the
Seyfert-like region of the line ratio diagnostic diagram [O iii]/H$\beta$
versus [S ii]/H$\alpha$ (Veilleux & Osterbrock 1987).

We infer for Jil2 an AGN bolometric luminosity log(Lbol/ [erg/s]) $\sim 45.8$,
from the narrow H$\beta$ luminosity (corrected for extinction; see Table 2),
following Netzer 2019. This result is consistent with the predictions reported
in Husemann et al. (2018b, 2021), to explain the presence of He ii$\lambda
1640$ in the Jil2 spectrum. All the arguments raised so far therefore further
support the scenario of a dual QSO in this complex system at $z\sim 3.3$.

Figure 14:  BPT diagnostic diagram. Red points represent flux ratios inferred
from the integrated Jil spectra, while the blue point indicates the QSO host
ratios. For Jil1, Jil5, and Jil7, [O iii]/H$\beta$ ratios are reported outside
of the BPT, as H$\alpha$ and [N ii] are undetected for these companions. Local
galaxies from SDSS DR7 (Abazajian et al. 2009) are indicated in grey, while
small stars represent model predictions for low-metallicity AGN from Nakajima
& Maiolino (2022, see this paper for a plethora of physical parameters related
to gas and AGN properties, such as ionisation and accretion disk temperature),
as labelled. The dashed line indicates the demarcation by Kauffmann et al.
(2003) between star-forming galaxies (left) and AGN (right) at low $z$; the
solid line from Kewley et al. (2001) includes more extreme starbursts and
composite objects among the star-forming galaxies at low $z$; the dot-dashed
grey line from Strom et al. (2017) shows the locus of star-forming galaxies at
$z\sim 2$.

####  6.4.2 Mergers as drivers for rapid SMBH growth?

Although the detailed physical connections among the eight companions - and
with the QSO host - is difficult to establish with the present data, it is
remarkable that LBQS 0302$-$0019 has this set of Jil galaxies within a
(projected) distance of $\sim 20$ kpc, all within a velocity range of $\sim\pm
250$ km s-1 from the QSO host systemic velocity. A blank field at $z\sim 3$ is
expected to have a space density of $\sim 0.01$ [O iii] emitters (with L([O
iii]) $>10^{41}$ erg/s) per Mpc-3 (Khostovan et al. 2015; Hirschmann et al.
2022); this corresponds to $5\times 10^{-4}$ expected galaxies within a $\sim
3\arcsec\times 3\arcsec$ region (the NIRSpec FOV), and within the narrow
redshift range associated with the Jil companions ($z=3.286-3.290$). We
conclude therefore that LBQS 0302$-$0019 is sitting in a ultra-dense
environment, being its space density many orders of magnitude higher than the
general field.

Interestingly, both ground- and space-based observations of $z>3$ QSOs have
shown that the presence of companions is common: for instance, sub-millimetre
galaxies and Ly$\alpha$ emitters in the vicinity of high-$z$ QSOs have been
identified with ALMA (e.g. Trakhtenbrot et al. 2017; Venemans et al. 2020;
Bischetti et al. 2021; García-Vergara et al. 2022) and MUSE (e.g. Fossati et
al. 2021), respectively. Indeed, almost all luminous high-$z$ QSOs so far
observed with JWST/NIRSpec IFS (LBQS 0302$-$0019; SDSS J1652+1728 in Wylezalek
et al. 2022; DELS J0411$-$0907 and VDES J0020$-$3653 in Marshall et al. 2023;
GS_3073 in Übler et al. 2023) and JWST/NIRCam WFSS (SDSS J0100+2802 in Kashino
et al. 2022) are surrounded by newly discovered companions.

These results clearly support the idea that mergers can be important drivers
for rapid early SMBH growth (e.g. Hopkins et al. 2008; Zana et al. 2022).
Indeed, NIRSpec IFS, thanks to its high sensitivity and angular resolution
($\sim 0.8$ kpc in a FOV of $25\times 25$ kpc2 at $z\sim 3$), is revealing
tidal bridges and tails at kiloparsec scales connecting such companions, hence
allowing the study of galaxy interactions at such high redshifts.

##  7 Conclusions

We have presented JWST/NIRSpec integral field spectroscopy of the blue QSO
LBQS 0302$-$0019 at $z=3.2870$. These observations cover a contiguous sky area
of $\sim 3\arcsec\times 3\arcsec$ (23 $\times$ 23 kpc2), which allowed us to
map the extension of the QSO host as well as characterise its environment with
a spatial sampling of $\sim 0.4$ kpc. The main results of our analysis
focussed on the QSO host are summarised below.

  * •

By analysing the integrated QSO spectrum, we measured the black hole mass from
the H$\beta$ and H$\alpha$ broad lines: $M_{\rm{BH}}\approx 2\times 10^{9}$
M⊙. With a bolometric luminosity of log(Lbol/ [erg s-1]) $\sim 47.2$, this QSO
is accreting material close to the Eddington limit ($\lambda_{\rm Edd}=0.9\pm
0.1$).

  * •

We have presented and make available for download a new procedure to model and
subtract the apparent wiggles in single-spaxel spectra due to the spatial
under-sampling of the PSF in NIRSpec IFS observations (see Figs. 5 and 6).
This correction is essential for performing spatial analyses of extended
emission sitting below a point source, such as for studies of QSO hosts and
close environments.

  * •

We performed a QSO–host decomposition using models of the QSO broad lines, and
used multi-component kinematic decomposition of the optical emission lines to
infer the physical properties of the emitting gas in the LBQS 0302$-$0019
host, as well as in its environment.

  * •

We revealed a broadly regular velocity field in the QSO host, which is
possibly tracing a warm rotating disk with $v_{\rm{rot}}/\sigma_{0}\approx 2$,
as inferred from 3D-Barolo modelling. We also derived a tentative dynamical
mass for the host, M${}_{\rm{dyn}}=(14\pm 6)\times 10^{10}$ M⊙; this places
our galaxy slightly above the local black hole–host mass relation (Kormendy &
Ho, 2013), consistent with other high-$z$ QSOs.

  * •

We identified a powerful outflow, with a velocity $v_{\rm{out}}\sim 1000$ km
s-1and a mass rate $\dot{M}_{\rm{out}}\sim 10^{4}$ M⊙ yr-1. Its kinetic and
momentum powers are compatible with the general predictions of AGN feedback
models (e.g. Harrison et al. 2018).

  * •

Standard BPT line ratios indicate that the central QSO dominates the
ionisation state of the gas, with no obvious sign of a contribution from young
stars in the host galaxy.

We also studied the complex, ultra-dense environment of LBQS 0302$-$0019
thanks to the large FOV of our IFS observations, covering three out of the
four companions already discovered by Husemann et al. (2021). Our main results
are as follows.

  * •

We detected eight Jil companion objects close to LBQS 0302$-$0019, three of
which were already discovered with MUSE and HST observations (Husemann et al.
2018b, 2021), for a total of nine companions within 30 kpc of the QSO. All of
these companions are within $\pm 250$ km s-1 of the QSO systemic velocity.

  * •

Regular velocity gradients, possibly tracing rotating gas, were detected in
Jil2 and Jil3. For these targets, we derived tentative dynamical masses of the
order of $10^{10}$ M⊙. However, we caution that the observed velocity
gradients may also be due to merger processes between different companions.

  * •

Though difficult to determine, some morpho-kinematic structures suggest that
the Jil companions may be connected with the QSO LBQS 0302$-$0019, so we can
speculate that they contribute to its feeding. In particular, Jil5 shows
evidence of gravitational interaction with the QSO host.

  * •

All BPT line ratios measured for Jil companions are compatible with AGN
ionisation.

  * •

We provide further evidence for the presence of an obscured QSO at $\sim 20$
kpc from LBQS 0302$-$0019 on the basis of [O iii]/H$\beta$, [S ii]/H$\alpha$,
and He ii/H$\beta$ line ratios. This QSO is likely responsible for the gas
ionisation in the surroundings of Jil2.

This work has explicitly demonstrated the exceptional capabilities of the
JWST/NIRSpec IFS to study the QSO environments in the early Universe. With a
total exposure time of $\sim 1$ hour, we unveiled in unprecedented detail the
interstellar properties of the LBQS 0302$-$0019 host galaxy and those of its
multiple companions in its immediate vicinity.

The study of the LBQS 0302$-$0019 host galaxy was limited by PSF artefacts;
before we could subtract them, we had to address the wiggles. We have shown
that wiggles can be modelled and subtracted, taking advantage of the fact that
their frequency, $f_{w}$, changes smoothly as a function of the wavelength
and, most importantly, that $f_{w}$ does not show spaxel-to-spaxel variations.
However, this step adds further difficulties in the analysis of the NIRSpec
data cubes. We note that the amplitude of these artefacts decreases as the
number of exposures increases. This information should be taken into
consideration by observers when planning NIRSpec IFS observations.

###### Acknowledgements.

This work is based on observations made with the NASA/ESA/CSA James Webb Space
Telescope. The data were obtained from the Mikulski Archive for Space
Telescopes at the Space Telescope Science Institute, which is operated by the
Association of Universities for Research in Astronomy, Inc., under NASA
contract NAS 5-03127 for JWST. These observations are associated with program
#1220, as part of the NIRSpec Galaxy Assembly IFS GTO program. We are grateful
to the anonymous referee for a constructive report that helped to improve the
quality of this manuscript. We thank Kimihiko Nakajima for providing the
theoretical model grids published by Nakajima & Maiolino (2022), and David Law
and Bartolomeo Trefoloni for helpful comments on an earlier version of this
manuscript. MP, SA, and BRP acknowledge support from the research project
PID2021-127718NB-I00 of the Spanish Ministry of Science and Innovation/State
Agency of Research (MICIN/AEI). MP also acknowledges support from the Programa
Atracción de Talento de la Comunidad de Madrid via grant 2018-T2/TIC-11715.
MAM acknowledges the support of a National Research Council of Canada Plaskett
Fellowship, and the Australian Research Council Centre of Excellence for All
Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number
CE170100013. RM, JS and FDE acknowledge support by the Science and Technology
Facilities Council (STFC), from the ERC Advanced Grant 695671 “QUENCH”. RM and
JS also acknowledge funding from a research professorship from the Royal
Society. GC acknowledges 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Ü gratefully acknowledges support by the Isaac
Newton Trust and by the Kavli Foundation through a Newton-Kavli Junior
Fellowship. AJB, GCJ and AJC acknowledge funding from the “FirstGalaxies”
Advanced Grant from the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation program (Grant agreement No.
789056). SC acknowledges support from the European Union (ERC,
WINGS,101040227). PGP-G acknowledges support from Spanish Ministerio de
Ciencia e Innovación MCIN/AEI/10.13039/501100011033 through grant
PGC2018-093499-B-I00. IL acknowledges support from PID2022-140483NB-C22 funded
by AEI 10.13039/501100011033 and BDC 20221289 funded by MCIN by the Recovery,
Transformation and Resilience Plan from the Spanish State, and by
NextGenerationEU from the European Union through the Recovery and Resilience
Facility.
This research has made use of NASA’s Astrophysics Data System, QFitsView, and
SAOImageDS9, developed by Smithsonian Astrophysical Observatory. It also had
made use of Python packages and software AstroPy (Astropy Collaboration et
al., 2013), Matplotlib (Hunter, 2007), NumPy (van der Walt et al., 2011),
OpenCV (Bradski, 2000), opencv-python, Photutils (Bradley et al., 2018),
Regions (Bradley et al., 2022), QDeblend3D (Husemann et al., 2013, 2014).

## References

  * Abazajian et al. (2009) Abazajian, K. N., Adelman-McCarthy, J. K., Agüeros, M. A., et al. 2009, ApJS, 182, 543
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Bacon et al. (2010) Bacon, R., Accardo, M., Adjali, L., et al. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III, ed. I. S. McLean, S. K. Ramsay, & H. Takami, 773508
  * Baldwin et al. (1981) Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5
  * Baron & Netzer (2019) Baron, D. & Netzer, H. 2019, MNRAS, 486, 4290
  * Baron et al. (2016) Baron, D., Stern, J., Poznanski, D., & Netzer, H. 2016, ApJ, 832, 8
  * Bischetti et al. (2021) Bischetti, M., Feruglio, C., Piconcelli, E., et al. 2021, A&A, 645, A33
  * Bischetti et al. (2017) Bischetti, M., Piconcelli, E., Vietri, G., et al. 2017, A&A, 598, A122
  * Böker et al. (2022) Böker, T., Arribas, S., Lützgendorf, N., et al. 2022, A&A, 661, A82
  * Böker et al. (2023) Böker, T., Beck, T. L., Birkmann, S. M., et al. 2023, PASP, 135, 038001
  * Bradley et al. (2022) Bradley, L., Deil, C., Ginsburg, A., et al. 2022, astropy/regions: v0.7
  * Bradley et al. (2018) Bradley, L., Sipocz, B., Robitaille, T., et al. 2018, astropy/photutils: v0.5
  * Bradski (2000) Bradski, G. 2000, Dr. Dobb’s Journal of Software Tools
  * Brinchmann (2022) Brinchmann, J. 2022, arXiv e-prints, arXiv:2208.07467
  * Brusa et al. (2015) Brusa, M., Bongiorno, A., Cresci, G., et al. 2015, MNRAS, 446, 2394
  * Bunker et al. (2023) Bunker, A. J., Saxena, A., Cameron, A. J., et al. 2023, arXiv e-prints, arXiv:2302.07256
  * Cameron et al. (2023) Cameron, A. J., Saxena, A., Bunker, A. J., et al. 2023, arXiv e-prints, arXiv:2302.04298
  * Cappellari (2017) Cappellari, M. 2017, MNRAS, 466, 798
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245
  * Carniani et al. (2015) Carniani, S., Marconi, A., Maiolino, R., et al. 2015, A&A, 580, A102
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Chen et al. (2022a) Chen, Y.-C., Hwang, H.-C., Shen, Y., et al. 2022a, ApJ, 925, 162
  * Chen et al. (2022b) Chen, Y.-C., Liu, X., Foord, A., et al. 2022b, arXiv e-prints, arXiv:2209.11249
  * Coatman et al. (2017) Coatman, L., Hewett, P. C., Banerji, M., et al. 2017, MNRAS, 465, 2120
  * Coatman et al. (2019) Coatman, L., Hewett, P. C., Banerji, M., et al. 2019, MNRAS, 486, 5335
  * Colpi (2014) Colpi, M. 2014, Space Sci. Rev., 183, 189
  * Cresci et al. (2015) Cresci, G., Mainieri, V., Brusa, M., et al. 2015, ApJ, 799, 82
  * Cresci et al. (2023) Cresci, G., Tozzi, G., Perna, M., et al. 2023, arXiv e-prints, arXiv:2301.11060
  * Curti et al. (2023) Curti, M., D’Eugenio, F., Carniani, S., et al. 2023, MNRAS, 518, 425
  * Curtis-Lake et al. (2022) Curtis-Lake, E., Carniani, S., Cameron, A., et al. 2022, arXiv e-prints, arXiv:2212.04568
  * Dalla Bontà et al. (2020) Dalla Bontà, E., Peterson, B. M., Bentz, M. C., et al. 2020, ApJ, 903, 112
  * D’Eugenio et al. (2023) D’Eugenio, F., Perez-Gonzalez, P., Maiolino, R., et al. 2023, arXiv e-prints, arXiv:2308.06317
  * Di Teodoro & Fraternali (2015) Di Teodoro, E. M. & Fraternali, F. 2015, MNRAS, 451, 3021
  * Dong et al. (2008) Dong, X., Wang, T., Wang, J., et al. 2008, MNRAS, 383, 581
  * Dorner et al. (2016) Dorner, B., Giardino, G., Ferruit, P., et al. 2016, A&A, 592, A113
  * Duras et al. (2020) Duras, F., Bongiorno, A., Ricci, F., et al. 2020, A&A, 636, A73
  * Fiore et al. (2017) Fiore, F., Feruglio, C., Shankar, F., et al. 2017, A&A, 601, A143
  * Fluetsch et al. (2019) Fluetsch, A., Maiolino, R., Carniani, S., et al. 2019, MNRAS, 483, 4586
  * Förster Schreiber et al. (2018) Förster Schreiber, N. M., Renzini, A., Mancini, C., et al. 2018, ApJS, 238, 21
  * Förster Schreiber et al. (2019) Förster Schreiber, N. M., Übler, H., Davies, R. L., et al. 2019, ApJ, 875, 21
  * Förster Schreiber & Wuyts (2020) Förster Schreiber, N. M. & Wuyts, S. 2020, ARA&A, 58, 661
  * Fossati et al. (2021) Fossati, M., Fumagalli, M., Lofthouse, E. K., et al. 2021, MNRAS, 503, 3044
  * García-Vergara et al. (2022) García-Vergara, C., Rybak, M., Hodge, J., et al. 2022, ApJ, 927, 65
  * Gaskell (1983) Gaskell, C. M. 1983, in Liege International Astrophysical Colloquia, Vol. 24, Liege International Astrophysical Colloquia, ed. J.-P. Swings, 473–477
  * Gaskell (2010) Gaskell, C. M. 2010, Nature, 463, E1
  * Greene & Ho (2006) Greene, J. E. & Ho, L. C. 2006, ApJ, 641, 117
  * Greene et al. (2012) Greene, J. E., Zakamska, N. L., & Smith, P. S. 2012, ApJ, 746, 86
  * Groves et al. (2004) Groves, B. A., Dopita, M. A., & Sutherland, R. S. 2004, ApJS, 153, 75
  * Harrison et al. (2012) Harrison, C. M., Alexander, D. M., Swinbank, A. M., et al. 2012, MNRAS, 426, 1073
  * Harrison et al. (2018) Harrison, C. M., Costa, T., Tadhunter, C. N., et al. 2018, Nature Astronomy, 2, 198
  * Hirschmann et al. (2022) Hirschmann, M., Charlot, S., Feltre, A., et al. 2022, arXiv e-prints, arXiv:2212.02522
  * Hopkins et al. (2008) Hopkins, P. F., Hernquist, L., Cox, T. J., & Kereš, D. 2008, ApJS, 175, 356
  * Hopkins et al. (2007) Hopkins, P. F., Richards, G. T., & Hernquist, L. 2007, ApJ, 654, 731
  * Hunter (2007) Hunter, J. D. 2007, CiSE, 9, 90
  * Husemann et al. (2018a) Husemann, B., Bielby, R., Jahnke, K., et al. 2018a, A&A, 614, L2
  * Husemann et al. (2014) Husemann, B., Jahnke, K., Sánchez, S. F., et al. 2014, MNRAS, 443, 755
  * Husemann et al. (2013) Husemann, B., Wisotzki, L., Sánchez, S. F., & Jahnke, K. 2013, A&A, 549, A43
  * Husemann et al. (2021) Husemann, B., Worseck, G., Arrigoni Battaia, F., Sander, A. A. C., & Shanks, T. 2021, A&A, 653, A122
  * Husemann et al. (2018b) Husemann, B., Worseck, G., Arrigoni-Battaia, F., & Shanks, T. 2018b, A&A, 610, L7
  * Jakobsen et al. (1994) Jakobsen, P., Boksenberg, A., Deharveng, J. M., et al. 1994, Nature, 370, 35
  * Jakobsen et al. (2022) Jakobsen, P., Ferruit, P., de Oliveira, C. A., et al. 2022, A&A, 661, A80
  * Ju et al. (2013) Ju, W., Greene, J. E., Rafikov, R. R., Bickerton, S. J., & Badenes, C. 2013, ApJ, 777, 44
  * Kakkad et al. (2020) Kakkad, D., Mainieri, V., Vietri, G., et al. 2020, A&A, 642, A147
  * Kashino et al. (2022) Kashino, D., Lilly, S. J., Matthee, J., et al. 2022 [arXiv:2211.08254]
  * Kauffmann et al. (2003) Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055
  * Kewley et al. (2001) Kewley, L. J., Dopita, M. A., Sutherland, R. S., Heisler, C. A., & Trevena, J. 2001, ApJ, 556, 121
  * Kewley et al. (2013) Kewley, L. J., Maier, C., Yabe, K., et al. 2013, ApJ, 774, L10
  * Khostovan et al. (2015) Khostovan, A. A., Sobral, D., Mobasher, B., et al. 2015, MNRAS, 452, 3948
  * Kocevski et al. (2023) Kocevski, D. D., Onoue, M., Inayoshi, K., et al. 2023, arXiv e-prints, arXiv:2302.00012
  * Komossa et al. (2008) Komossa, S., Zhou, H., & Lu, H. 2008, ApJ, 678, L81
  * Kormendy & Ho (2013) Kormendy, J. & Ho, L. C. 2013, ARA&A, 51, 511
  * Kovacevic et al. (2010) Kovacevic, J., Popovic, L. C., & Dimitrijevic, M. S. 2010, ApJS, 189, 15
  * Law et al. (2023) Law, D. D., Morrison, J. E., Argyriou, I., et al. 2023, AJ, 166, 45
  * Lemon et al. (2022) Lemon, C., Millon, M., Sluse, D., et al. 2022, A&A, 657, A113
  * Mannucci et al. (2022) Mannucci, F., Pancino, E., Belfiore, F., et al. 2022, Nature Astronomy, 6, 1185
  * Marasco et al. (2020) Marasco, A., Cresci, G., Nardini, E., et al. 2020, A&A, 644, A15
  * Marshall et al. (2023) Marshall, M. A., Perna, M., Willott, C. J., et al. 2023, arXiv e-prints, arXiv:2302.04795
  * Meena et al. (2021) Meena, B., Crenshaw, D. M., Schmitt, H. R., et al. 2021, ApJ, 916, 31
  * Müller-Sánchez et al. (2011) Müller-Sánchez, F., Prieto, M. A., Hicks, E. K. S., et al. 2011, ApJ, 739, 69
  * Nagao et al. (2006) Nagao, T., Marconi, A., & Maiolino, R. 2006, A&A, 447, 157
  * Nakajima & Maiolino (2022) Nakajima, K. & Maiolino, R. 2022, MNRAS, 513, 5134
  * Nardini et al. (2019) Nardini, E., Lusso, E., Risaliti, G., et al. 2019, A&A, 632, A109
  * Netzer (2019) Netzer, H. 2019, MNRAS, 488, 5185
  * Osterbrock & Ferland (2006) Osterbrock, D. E. & Ferland, G. J. 2006, Astrophysics of gaseous nebulae and active galactic nuclei (University Science Books)
  * Perna et al. (2020) Perna, M., Arribas, S., Catalán-Torrecilla, C., et al. 2020, A&A, 643, A139
  * Perna et al. (2022) Perna, M., Arribas, S., Colina, L., et al. 2022, A&A, 662, A94
  * Perna et al. (2015a) Perna, M., Brusa, M., Cresci, G., et al. 2015a, A&A, 574, A82
  * Perna et al. (2015b) Perna, M., Brusa, M., Salvato, M., et al. 2015b, A&A, 583, A72
  * Perna et al. (2019) Perna, M., Cresci, G., Brusa, M., et al. 2019, A&A, 623, A171
  * Perna et al. (2017) Perna, M., Lanzuisi, G., Brusa, M., Cresci, G., & Mignoli, M. 2017, A&A, 606, A96
  * Rigby et al. (2022) Rigby, J., Perrin, M., McElwain, M., et al. 2022 [arXiv:2207.05632]
  * Robertson et al. (2022) Robertson, B. E., Tacchella, S., Johnson, B. D., et al. 2022, arXiv e-prints, arXiv:2212.04480
  * Schwarz (1978) Schwarz, U. J. 1978, A&A, 65, 345
  * Seaton (1978) Seaton, M. J. 1978, MNRAS, 185, 5P
  * Shen (2016) Shen, Y. 2016, ApJ, 817, 55
  * Shirazi & Brinchmann (2012) Shirazi, M. & Brinchmann, J. 2012, MNRAS, 421, 1043
  * Smith et al. (2007) Smith, J. D. T., Armus, L., Dale, D. A., et al. 2007, PASP, 119, 1133
  * Strom et al. (2017) Strom, A. L., Steidel, C. C., Rudie, G. C., et al. 2017, ApJ, 836, 164
  * Tacchella et al. (2022) Tacchella, S., Johnson, B. D., Robertson, B. E., et al. 2022, arXiv e-prints, arXiv:2208.03281
  * Tozzi et al. (2021) Tozzi, G., Cresci, G., Marasco, A., et al. 2021, A&A, 648, A99
  * Tozzi et al. (2023) Tozzi, G., Maiolino, R., Cresci, G., et al. 2023, MNRAS, 521, 1264
  * Trakhtenbrot et al. (2017) Trakhtenbrot, B., Lira, P., Netzer, H., et al. 2017, ApJ, 836, 8
  * Trefoloni et al. (2023) Trefoloni, B., Lusso, E., Nardini, E., et al. 2023, arXiv e-prints, arXiv:2305.07699
  * Übler et al. (2023) Übler, H., Maiolino, R., Curtis-Lake, E., et al. 2023, arXiv e-prints, arXiv:2302.06647
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, CiSE, 13, 22
  * van Dokkum (2001) van Dokkum, P. G. 2001, PASP, 113, 1420
  * Van Wassenhove et al. (2012) Van Wassenhove, S., Volonteri, M., Mayer, L., et al. 2012, ApJ, 748, L7
  * Vayner et al. (2023) Vayner, A., Zakamska, N. L., Ishikawa, Y., et al. 2023, arXiv e-prints, arXiv:2303.06970
  * Veilleux & Osterbrock (1987) Veilleux, S. & Osterbrock, D. E. 1987, ApJS, 63, 295
  * Venemans et al. (2020) Venemans, B. P., Walter, F., Neeleman, M., et al. 2020, ApJ, 904, 130
  * Vietri et al. (2020) Vietri, G., Mainieri, V., Kakkad, D., et al. 2020, A&A, 644, A175
  * Villar Martín et al. (2020) Villar Martín, M., Perna, M., Humphrey, A., et al. 2020, A&A, 634, A116
  * Volonteri et al. (2021) Volonteri, M., Habouzit, M., & Colpi, M. 2021, Nature Reviews Physics, 3, 732
  * Wisnioski et al. (2019) Wisnioski, E., Förster Schreiber, N. M., Fossati, M., et al. 2019, ApJ, 886, 124
  * Woo et al. (2016) Woo, J.-H., Bae, H.-J., Son, D., & Karouzos, M. 2016, ApJ, 817, 108
  * Worseck et al. (2021) Worseck, G., Khrykin, I. S., Hennawi, J. F., Prochaska, J. X., & Farina, E. P. 2021, MNRAS, 505, 5084
  * Wylezalek et al. (2022) Wylezalek, D., Vayner, A., Rupke, D. S. N., et al. 2022, ApJ, 940, L7
  * Zana et al. (2022) Zana, T., Gallerani, S., Carniani, S., et al. 2022, MNRAS, 513, 2118
  * Zuo et al. (2015) Zuo, W., Wu, X.-B., Fan, X., et al. 2015, ApJ, 799, 189
  * Zuo et al. (2020) Zuo, W., Wu, X.-B., Fan, X., et al. 2020, ApJ, 896, 40

##  Appendix A Nuclear spectra

Figure 15 shows the simultaneously fit of four spectra extracted from circular
regions with radius of 0.2″(4 spaxels) and centred at different positions
within a few spaxels from the peak emission of the QSO, used to reduce the
degeneracy between BLR and NLR. All spectra are normalised so that the BLR
wings of the H$\alpha$ and H$\beta$ have the same fluxes, and can be fitted
with the same broken power-law functions. During the fit, BLR profiles are
therefore tied, assuming that these emission components originate from the
same unresolved region. All other components are free to vary as they
originate from more extended and likely resolved regions. The small aperture
radius is required to observe significant variations in the H$\alpha$-[N ii]
complex (e.g. with respect to the integrated spectrum in Fig. 3).

Figure 15:  Integrated spectra extracted from circular regions with a radius
of 0.2″and centred at different positions within a few spaxels of the peak
emission of the QSO. The best-fit models shown here were obtained by fitting
the four spectra with the same BLR profiles, as explained in Sect. 4.

##  Appendix B Sinusoidal-type patterns

In Sect. 5.1, we proved that our approach is capable of modelling and
correcting for the spurious wiggles in the single-spaxel spectra. However, it
has important limitations. For instance, some residual wiggles are present in
a few spaxels (see e.g. Fig. 18). Moreover, in the very innermost nuclear
regions, the H$\alpha$ BLR emission covers a significant number of wavelength
channels, preventing a proper modellisation of the underlying wiggles in the
vicinity of the H$\alpha$ line (see Fig. 6). This can affect the
reconstruction of the H$\alpha$ kinematics.

Another aspect is related to the emission line fluxes: an improper correction
of the wiggles implies an incorrect reconstruction of the emission line
profile and, as a consequence, an incorrect measurement of its integrated
flux. For LBQS 0302$-$0019, we check that the [O iii]$\lambda\lambda
4959,5007$ line ratio is preserved at 1:3, which is consistent with theory
(Osterbrock & Ferland 2006). Figure 20 shows the nuclear spectra extracted
from different areas (integrating over circular regions with radius from 1 to
5 spaxels), from the original cube (top panel) and the one corrected for
wiggles (bottom). All spectra are continuum-subtracted and normalised to the
[O iii] peak; the inset in the bottom panel shows that the [O iii]$\lambda
4959$ peaks at $\sim 0.33$, consistent with the expectations. We note however
that significant deviations (up to $\sim 50\%$) are observed in individual
spaxels, both in the original and in the corrected spectra, although the
corrected ones have line ratios closer to the theoretical 1:3 ratio. We also
checked that our corrections preserve the shape of the spectrum and integrated
fluxes, as shown in Fig. 21.

Therefore, we caution that the presence of wiggles might affect both the
kinematics and flux ratio measurements; a proper modellisation and subtraction
of the wiggles is required to mitigate their effects. In fact, as shown in
Fig. 20, off-centred integrated spectra are always affected by these wiggles.

Figure 16: Sinusoidal-type patterns in single-spaxel spectra extracted from
the emsm data cube with a spaxel size of 0.1″. Top panel: LBQS 0302$-$0019
spectrum integrated over an aperture of $r=0.5$″(orange curve), in comparison
with the spectrum of the brightest spaxel (blue curve). Both spectra are
normalised to their maximum values for visualisation purposes. The wiggles
affecting the single-spaxel spectrum are reported in grey and are obtained as
the difference between the blue and orange curves (see Fig. 5 for details).
Bottom panel: Wiggles obtained from the eight pixels closest to the brightest
one.  Figure 17: Wiggle-corrected spectra extracted from the emsm data cube
with a spaxel size of 0.1″. Top panel: LBQS 0302$-$0019 spectrum integrated
over an aperture of $r=0.5$″(orange curve), in comparison with the spectrum of
the brightest spaxel, after the wiggle subtraction (blue curve). Both spectra
are normalised to 1 for visualisation purposes. The residuals are reported in
grey and are obtained as the difference between the blue and orange curves
(see Fig. 5 for details). Bottom panel: Residuals obtained from the eight
spaxels closest to the brightest one. The most significant residuals are found
at the position of the brightest emission lines: they are not due to the
wiggles, but to the line profile variations.  Figure 18: Same as Fig. 17, but
for the drizzle data cube, with spaxels of 0.05″.  Figure 19: Same as Fig. 6,
but for the drizzle data cube of the QSO VDES J0020-3656 (Marshall et al.
2023), with spaxels of 0.05″.

Figure 20: Integrated spectra extracted from circular regions containing 1 to
49 spaxels (corresponding to radii of 1 to 5 spaxels), centred at 0.3″east of
the LBQS 0302$-$0019 nucleus (from the drizzle cubes with spaxels of 0.05″).
The top panel shows the original spectra, while the bottom panel shows the
same spectra after the correction for the wiggles at the spaxel level (Sect.
5.1). All spectra are continuum-subtracted and are normalised to the peak of
[O iii]; for those extracted from regions with radii $<5$ spaxels, we added
vertical offsets to ease the visualisation. The insets show a zoomed-in view
of the vicinity of the [O iii] and H$\beta$ lines, without any vertical
offset; these spectra show that the [O iii]$\lambda 4959$ peaks at $\sim 0.33$
(indicated by the horizontal dashed line), consistent with theoretical
expectations.  Figure 21:  Integrated spectra extracted from circular regions
containing 5 to 49 spaxels (corresponding to radii of 2 to 5 spaxels), centred
on the LBQS 0302$-$0019 nucleus (from the drizzle cubes with spaxels of
0.05″). The solid lines show the spectra after the wiggle subtraction, while
the dotted lines show the original spectra. All spectra are normalised to the
peak of [O iii]; for those extracted from regions with radii $<5$ spaxels, we
added vertical offsets to ease the visualisation. The figure proves that our
correction preserves the integrated fluxes and the shape of the spectrum.

##  Appendix C 3D-Barolo fit

Figures 22, 23, and 24 show the 3D-Barolo best-fit modellisation for the three
targets that display broadly regular velocity gradients. We caution that the
significant residuals, likely due to the superposition of different kinematic
components associated with distinct clumps (or targets) on the same LOS, call
into question the reliability of the inferred best-fit parameters.

Figure 22: LBQS 0302$-$0019 host galaxy disk kinematic best fit of Moment 0,
1, and 2 (first to third rows). These best fits are inferred from the analysis
of the narrow [O iii] component obtained from our multi-component Gaussian fit
(i.e. all components with FWHM $<600$ km s-1in Fig. 9). The black and green
lines identify the major axis and the zero-velocity curve, respectively; the
red cross identifies the QSO position.  Figure 23: Jil2 Moment 0, 1, and 2 and
the 3D-Barolo disk kinematic best fit. See Fig. 22 for details.  Figure 24:
Jil3 Moment 0, 1, and 2 and the 3D-Barolo disk kinematic best fit. See Fig. 22
for details.