The ultra-dense, interacting environment of a dual AGN at z 3.3 revealed by JWST/NIRSpec IFS

: The ultra-dense, interacting environment of a dual AGN at z 3.3 revealed by JWST/NIRSpec IFS


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β, [O iii] λλ4959,5007, Hα, and [N ii] λλ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 ∼ 2−3; e.g.Förster Schreiber & Wuyts 2020) to the epoch of re-ionisation (z 7; e.g.Robertson et al. 2023;Curtis-Lake et al. 2023).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 2023; Bunker et al. 2023;Cameron et al. 2023;Cresci et al. 2023;Curti et al. 2023;Kocevski et al. 2023;Tacchella et al. 2023;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 highluminosity, 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. 2022Chen et al. , 2023;;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. 2018a).
The QSO LBQS 0302−0019 (RA 3 h 4 m 49.93 s , Dec −0 • 8 13.10 , J2000) at z ∼ 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α 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. 2018a) report the detection of a Lyα 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α, C ivλ1549, He iiλ1640, and C iiiλ1909 ultraviolet emission-line diagnostics.The He ii line luminosity, L(He ii) ∼ 1.7 × 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 AGN ∼ 10 45 erg s −1 ), if located within the compact region emitting He ii.
Follow-up ground-based K s-band imaging and near-infrared spectroscopy are presented in Husemann et al. (2018b), who successfully detected Jil's host galaxy emission, with an estimated stellar mass of ∼10 11 M , and the optical [O iii]λ5007 line ([O iii] hereinafter), with L([O iii]) ∼ 2.5×10 42 erg s −1 .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 ∼1 (∼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 +222 −109 Myr and log(M * /M ) = 11.2 +0.3 −0.1 , and Jil2, associated with the compact He ii emission, with t * = 19 +74 −14 Myr and log(M * /M ) = 9.4 +0.9 −0.4 .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 1000 km s −1 .A velocity offset of C iv relative to the centroid of the Hβ 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.NIR-Spec 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 Sect.7.
Throughout, we adopt a Chabrier (2003) initial mass function (0.1−100 M ) and a flat Λ cold dark matter cosmology with H 0 = 70 km s −1 Mpc −1 , Ω Λ = 0.7, and Ω 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.

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 A89, page 2 of 23 NIRSpec's IFS mode, which provides spatially resolved spectroscopy over a contiguous 3.1 × 3.2 sky area, with a sampling of 0.1 /spaxel and a spatial resolution from ∼0.04 (at ∼1 µm) to ∼0.15 (at ∼5 µm; see Böker et al. 2022;Rigby et al. 2023).The IFS observations were taken with the grating/filter pair G235H/F170LP.This results in a data cube with spectral resolution R ∼ 2700 over the wavelength range 1.7−3.1 µ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 s.

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 detectorlevel 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: -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 countrate 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 σ-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 (PI 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 quality 1 .We therefore decided to follow two different approaches to remove the spikes: the first one uses an algorithm similar to lacosmic The second approach consists of a post-processing correction, and is done applying a σ 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 /0.1 ) 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.0 2 ; we also patched the photom script, applying the corrections implemented in the same version 3 , which allows more reasonable flux densities to be inferred (i.e. a factor of ∼100 smaller with respect to those obtained with standard pipeline 1.8.2).

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 ∆RA = −0.492and ∆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.pmap 4.
1 At the time of this writing, the newest version of the pipeline, v1.9.4, and the latest context file, jwst_1063.pmapare still affected by these issues. 2 jwst_1063.pmapcorrects for a ∼4 pixels systematic offset associated with the coordinate transformation between the "OTEIP" and the world systems, but not for a smaller offset (∼0.2−0.4 pixels) between the "GWA" and the "virtual slit" frame (see Dorner et al. 2016).

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β 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 ∼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 calibration5 .(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α 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.
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 σ-clipping method.

Spectral analysis of the integrated
LBQS 0302-0019 spectrum

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α and Hβ lines, the [O iii] λλ4959,5007, [N ii] λλ6548,83, and [S ii] λλ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α and Hβ 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 Å.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β-[O iii] and Hα-[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. A.1) 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 Hβ.All these previous measurements are within ≈±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).

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 ± 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.

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α and Hβ BLR region following the single-epoch calibrations from Dalla Bontà et al. ( 2020), with an intrinsic scatter of ∼0.3 dex, and from Greene & Ho (2006): (2) with larger intrinsic scatters of ∼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α/Hβ flux ratio of 3.88 +0.19 −0.12 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 ∼500 low-z Seyfert 1 and QSOs with minimal dust extinction effects, Hα/Hβ = 3.06 ± 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 M BH .
The Balmer line luminosities and widths are measured from our best-fit BLR profiles shown in Fig. 3 (i.e. the broken powerlaw components); we obtain estimates of the black hole mass of the order of ∼2 × 10 9 M .These values are broadly consistent with those previously reported in the literature and are based on Hβ and C iv BLR measurements (with the latter being slightly larger, as commonly reported in the literature; e.g.Coatman et al. 2017).
Using our black hole mass estimate from the Hβ BLR (Eq.( 1), which has smaller scatter than Eqs.( 2) or ( 3)), we find an Eddington ratio of λ Edd = 0.9 ± 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.
−0.11 × 10 9 log(L bol /(erg s −1 )) DB20 47.22 ± 0.01 λ Edd 0.9 ± 0.1 log(L [OIII] /(erg s −1 )) 45.71 ± 0.03    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β profile and the wings of the [O iii] lines.See used in the literature: V10 ∼ −760 km s −1 , the velocity at the 10th percentile of the overall emission-line profile, W80 ∼ 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 ∼ 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.

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 1D 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 (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 ).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 ∼ 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. B.1), although the amplitude of their wiggles is ≈2−3 times smaller than in the drizzle cubes.Moreover, the use of emsm implies a decrease in spatial resolution, down to ∼0.2 (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 download6 .

Modelling of the wiggles
Figures 5 and B.1 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 × 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π f w w + φ) + B, where A is the amplitude, f w is the frequency in 1/µm, w is the wavelength, φ is the phase shift, and B is the continuum level; we repeated the process in small portions of the wavelength range (∼0.1 µ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 signifi-cantly 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 ∼ 40 µm −1 at shortest and longest wavelengths, and f w ∼ 5 µ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 B.2 and B.3 show the same residuals presented in Figs. 5 and B.1, 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 with 160 • ).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.

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(Husemann et al. , 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α line, taking as a reference the Hα BLR emission, and one for those in the vicinity of Hβ, taking as a reference the Hβ 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 (Sect.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.
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α and Hβ.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 highvelocity 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.

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 A89, page 9 of 23 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 (∆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.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 ∼±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α velocity field is nois-  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 rot /σ 0 ratio, where v rot is the intrinsic maximum rotation velocity (corrected for inclination, v rot = v LOS / sin(i)) and σ 0 is the intrinsic velocity dispersion of the rotating disk, related to its thickness.In this work, we define σ 0 as the measured line width in the outer parts of the galaxy, corrected for the instrumental spectral resolution (e.g. more affected by noise).The rotation-to-random motion ratio v rot /σ 0 ≈ 2 indicates that this galaxy is associated with a dynamically warm disk, consistent with z ∼ 2 galaxies presented in Förster Schreiber et al. (2018), with v rot /σ 0 spanning the range from 0.97 to 13 (with a median of 3.2), as inferred from Hα gas kinematics (see also e.g.Wisnioski et al. 2019).

QSO host disk
The 3D-Barolo best-fit velocity maps also show significant residuals in the receding part, at ∼0.15 south-west of the nucleus, with velocities ≈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α line, if due to Hα line, would correspond to L(Hα) ≈ 10 43 erg s −1 , 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): M dyn = (14 ± 6) × 10 10 M , within a radius of 2.4 ± 0.6 kpc (corrected for the PSF, and containing 85% of the [O iii] total flux, as inferred from the QSOsubtracted cube).Combining this measurement with the M BH derived in Sect.4, we obtained a M BH /M dyn ≈ 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).

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β.We used the equation from Cresci et al. (2015), where L 41 (Hβ) is the Hβ luminosity associated with the outflow component in units of 10 41 erg s −1 , n e is the electron density, v out is the outflow velocity, and R 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 ∼3 higher mass rate would be obtained for instance using the electron density measured in the outflowing gas of the QSO XID2028 at z ∼ 1.5 (i.e.360 ± 180 cm −3 ), as measured from recent JWST/NIRSpec IFS observations (Cresci et al. 2023).
Here we consider the L 41 (Hβ) to be the luminosity of the Hβ 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β)/[erg s −1 ]) = 45.17 +0.09 −0.13 .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.07  −0.15 .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 (∼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 out = 930 +60 −110 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 out = t d × v out 9 kpc.This estimate is compatible with the extension of ionised outflows observed in other QSOs at high z, in the range ≈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 Ṁout (Hβ) ∼ 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 A89, page 11 of 23   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.

Target z (∆v [km s
Integrated [O iii] and Hα 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β detection, we measured the log([O iii]/Hβ) lower limit assuming that the Hβ upper limit is three times smaller than Hα.The non-parametric velocity W80 refers to the [O iii] line profile.momentum powers are Ėout = 1/2 Ṁout v 2 out ∼ 4 × 10 45 erg s −1 and Ṗout = Ṁout v out ∼ 8 × 10 37 dyne, respectively.Hence, the kinetic power is ∼2% of the radiative luminosity of the AGN, while the momentum rate is in excess of ∼15 times the radiative momentum flux (L bol /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).

Further considerations of the outflow extension
We report here two further arguments to better justify the assumed outflow extension (≈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 α out of a few degrees) would be required to explain the presence of a spatially unresolved ( 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 α out ≈ 10−60 • ; 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.Ṁout ∼ 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.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 ∼±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 wavelength (Å) Fig. 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.
can therefore speculate that this companion is contributing to the feeding of the QSO host.In this case, we detect a lower limit for log([O iii]/Hβ) > 0.7, consistent with the presence of high ionisation.The line ratio diagram also shows that Jil2 and Jil3 galaxies are associated with very stringent upper limits for the log([N ii]/Hα), of the order of −1.This may indicate that they are metal-poor AGN or galaxies, consistent with model predictions (Z 0.5 Z ; 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. 2018a).On the other hand, the LBQS 0302−0019 host might be associated with a higher metallicity (Z ≈ Z ; according to the same grid models), because of the higher [N ii]/Hα.Kauffmann et al. (2003) between starforming 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 ∼ 2.
The relative proximity of Jil5, Jil6, Jil7, and Jil8 to the QSO is a likely explanation for the high [O iii]/Hβ in such targets.On the other hand, the high [O iii]/Hβ in Jil1, Jil2, Jil3, and Jil9 can be explained by the presence of an AGN in Jil2, inferred by Husemann et al. (2018a) on the basis of ultraviolet diagnostics.In support of this scenario, we used the He ii λ4686 diagnostics (Shirazi & Brinchmann 2012;Nakajima & Maiolino 2022;Übler et al. 2023;Tozzi et al. 2023).Since the He ii λ4686 is undetected in NIRSpec, we used the ratio He ii λ1640/He ii λ4686 = 7.2, expected for recombination (Seaton 1978), to infer the He ii λ4686 flux in Jil2 (correcting for extinction).This gives A89, page 14 of 23 for Jil2 a log(He ii λ4686/Hβ) = −0.22,consistent with AGN ionisation (see Fig. 7 in Übler et al. 2023).We stress that the detection of He iiλ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 ∼2σ) 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σ detection in Jil2).
We infer for Jil2 an AGN bolometric luminosity log(L bol /[erg s −1 ]) ∼ 45.8, from the narrow Hβ luminosity (corrected for extinction; see Table 2), following (Netzer 2019).This result is consistent with the predictions reported in Husemann et al. (2018aHusemann et al. ( , 2021)), to explain the presence of He ii λ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 ∼ 3.3.

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 ∼20 kpc, all within a velocity range of ∼±250 km s −1 from the QSO host systemic velocity.A blank field at z ∼ 3 is expected to have a space density of ∼0.01 [O iii] emitters (with L([O iii]) >10 41 erg s −1 ) per Mpc −3 (Khostovan et al. 2015;Hirschmann et al. 2023); this corresponds to 5 × 10 −4 expected galaxies within a ∼3 × 3 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.
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 (∼0.8 kpc in a FOV of 25 × 25 kpc 2 at z ∼ 3), is revealing tidal bridges and tails at kiloparsec scales connecting such companions, hence allowing the study of galaxy interactions at such high redshifts.

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 ∼3 × 3 (23 × 23 kpc 2 ), which allowed us to map the extension of the QSO host as well as characterise its environment with a spatial sampling of ∼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β and Hα broad lines: M BH ≈ 2 × 10 9 M .With a bolometric luminosity of log(L bol /[erg s −1 ]) ∼ 47.2, this QSO is accreting material close to the Eddington limit (λ Edd = 0.9 ± 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 rot /σ 0 ≈ 2, as inferred from 3D-Barolo modelling.We also derived a tentative dynamical mass for the host, M dyn = (14 ± 6) × 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 out ∼ 1000 km s −1 and a mass rate Ṁout ∼ 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. 2018a(Husemann et al. , 2021)), for a total of nine companions within 30 kpc of the QSO.All of these companions are within ±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 ∼20 kpc from LBQS 0302−0019 on the basis of [O iii]/Hβ, [S ii]/Hα, and He ii/Hβ 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 A89, page 15 of 23 in the early Universe.With a total exposure time of ∼1 h, 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.Figure A.1 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α and Hβ 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α-[N ii] complex (e.g. with respect to the integrated spectrum in Fig. 3).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.    .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 A: Nuclear spectra
A89, page 22 of 23

(
van Dokkum 2001) to remove outliers in individual exposures (at the end of Stage 2): because our sources are undersampled 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× the noise, whichever was highest) and we rejected the 95th percentile of the resulting distribution (see D'Eugenio et al. 2023 for details).

Fig. 1 .
Fig. 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α emission before (top) and after (bottom) the QSO PSF subtraction, from Husemann et al. (2018a), with contours from the PSF-subtracted HST image.The right panels show the [O iii] λ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.

Fig. 2 .
Fig. 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β and [O iii] lines.

Fig. 3 .
Fig. 3. Multi-component simultaneous best-fit results for the continuum-subtracted spectrum of LBQS 0302−0019, around the Hβ-[O iii] (left) and Hα-[N ii] regions (right; integrated over a circular region with r = 0.5 ).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 Fig. A.1 show additional peaks and/or inflection points in the Hα-[N ii] complex, due to the presence of strong [N ii] emission line compo-nents; 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.3is likely due to the more prominent BLR emission, and the stronger degeneracy between BLR and NLR kinematic components.4.2.Systemic redshiftWe 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 ± 0.0003, which is in agreement withZuo et al. (2015Zuo et al. ( , 2020) ) but at odds with other redshift measurements from the literature.Husemann et al. (2018a) reported values in the range 3.2882−3.2887(for different ultraviolet lines);Coatman et al. (2019) reported z = 3.2856 ± 0.0002 for the [O iii], and z = 3.2868 ± 0.0012 for the

Fig. 4 .
Fig. 4. Comparison of the BLR profiles of the C iv (from Shen 2016) and Hα and Hβ (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.
Fig. B.1 for analogous effects in the emsm cube.

Fig. 6 .
Fig.6.Modelling of the wiggles in single-spaxel spectra.Top panel: integrated LBQS 0302−0019 spectrum (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.

Fig. 7 .
Fig. 7. Reconstructed PSF from the spatially unresolved BLR emission.Left: PSFs measured from the drizzle cube for the Hα and Hβ BLR emission, respectively, as described in Sect.5.3.The reconstructed Hβ PSF is less extended than Hα, 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β (bottom) and Hα (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.

Fig. 8 .Fig. 9 .
Fig. 8. PSFs measured from the drizzle cube for the nuclear emission around Hα (left) and Hβ (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.

Figure 10
Figure 10 shows an overview of the flux distribution and kinematics of the narrow component in the LBQS 0302−0019 host ier than the [O iii] one, because of the BLR subtraction step, and the degeneracy between Hα and [N ii] lines.The [O iii] and Hα 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α and [O iii] maps might be present.As the Hα 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. ( Fig. 10.From left to right: Hα 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.

Fig. 11 .
Fig. 11.From left to right: Hα 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.
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α (bottom).The flux distribution shows multiple clumps in the north-east regions, as well as plumes and irregular structures within ≈1 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 ∼±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.12we show a few narrow-band 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.
For this companion we measure a log([O iii]/Hβ) > 0.7, consistent with flux ratios measured in the QSO host, and hence likely ionised by the QSO radiation.Jil6 is located at ∼10 kpc south-east of the QSO, and is detected in [O iii] and Hα (and in Hβ and [N ii] at S /N ∼ 2−3).Both log([O iii]/Hβ) ∼ 0.7 and log([N ii]/Hα) ∼ 0.1 suggest a QSO ionisation.Jil7 is located at ∼10 kpc south-east of the QSO, and is detected in [O iii] and Hβ.It shows a prominent blue wing in the [O iii] (V10 ∼ −550 ± 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β) ∼ 0.52).Jil8 is located at ∼8 kpc east of the QSO nucleus, with an extension of ∼2 kpc.It is detected in [O iii], Hα, and Hβ.The broad components in the emission lines are due to PSF artefacts.For this companion, log([O iii]/Hβ) ∼ 0.7 suggests a QSO ionisation.Jil9 is located at ∼23 kpc north-east of the blue QSO, and is detected in [O iii] and Hα.It shows a velocity offset of ∼300 km s −1 from Jil3, and narrow line profiles (W80 ∼ 170 km s −1 ).
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α 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β ratios.

Fig. 14 .
Fig. 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β ratios are reported outside of the BPT, as Hα 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 lowmetallicity 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 byKauffmann et al. (2003) between starforming galaxies (left) and AGN (right) at low z; the solid line fromKewley et al. (2001) includes more extreme starbursts and composite objects among the star-forming galaxies at low z; the dot-dashed grey line fromStrom et al. (2017) shows the locus of star-forming galaxies at z ∼ 2.

Fig
Fig.A.1.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.
Fig. B.3.Same as Fig. B.2, but for the drizzle data cube, with spaxels of 0.05 .

Fig. B. 5 .
Fig. B.5. 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β lines, without any vertical offset; these spectra show that the [O iii]λ4959 peaks at ∼0.33 (indicated by the horizontal dashed line), consistent with theoretical expectations.

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).

Table 2 .
Properties of the companion sources in the LBQS 0302−0019 environment.