© The Authors 2023

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β, [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 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. 2022, 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^h4^m49.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⁴² erg s⁻¹, 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⁴⁵ erg s⁻¹), 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. (2018b), who successfully detected Jil’s host galaxy emission, with an estimated stellar mass of ∼10¹¹ M_⊙, and the optical [O III]λ5007 line ([O III] hereinafter), with L([O III]) ∼ 2.5 × 10⁴² 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 ∼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_{\ast }={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_{\ast }={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.

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.

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⁻¹. 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⁻¹ 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 Sect. 7.

Throughout, we adopt a Chabrier (2003) initial mass function (0.1 − 100 M_⊙) and a flat Λ cold dark matter cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_Λ = 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.

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″ × 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.

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:

• 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 σ-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¹. 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× 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 σ 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″)² 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²; we also patched the photom script, applying the corrections implemented in the same version³, 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).

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 #x0394;RA = −0.492″ and #x0394;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⁴.

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β and [O III] lines.

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.

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 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 σ-clipping method.

4. Spectral analysis of the integrated LBQS 0302–0019 spectrum

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

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

We note that the integrated spectra reported in 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 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  ±  0.0003, which is in agreement with Zuo et al. (2015, 2020) but at odds with other 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 Hβ. All these previous measurements are within ≈ ± 100 km s⁻¹ 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 ± 60 km s⁻¹. 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α and Hβ BLR region following the single-epoch calibrations from Dalla Bontà et al. (2020),

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

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

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

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

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.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 ∼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 power-law components); we obtain estimates of the black hole mass of the order of ∼2 × 10⁹ 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).

We calculated the bolometric luminosity of LBQS 0302−0019 following Dalla Bontà et al. (2020), hence using the Hβ BLR luminosity (their Eq. (25)): log(L_bol/[erg s⁻¹]) ∼47.2 (consistent with the value inferred from the continuum luminosity at 5100 Å, and using the (Netzer 2019) bolometric correction, ∼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(L_bol/[erg s⁻¹]) ∼ 46.7.

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.

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λ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⁻¹, 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).

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.

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 ∼ −760 km s⁻¹, the velocity at the 10th percentile of the overall emission-line profile, W80 ∼ 1080 km s⁻¹, 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⁻¹, 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.

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

Fig. 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 Fig. B.1 for analogous effects in the emsm cube.

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 download⁶.

5.2. 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) = Asin(2πf_ww + ϕ)+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.

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.

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 ∼ 40 μm⁻¹ at shortest and longest wavelengths, and f_w ∼ 5 μm⁻¹ 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 Fig. B.4 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° and on October 16, 2022, 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.

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

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.

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

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.

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⁻¹) associated with significant blueshifts (#x0394;v < −100 km s⁻¹), 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⁻¹). 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.

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

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α (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.

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.

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⁻¹, 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 noisier 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. (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/σ₀ ratio, where v_rot is the intrinsic maximum rotation velocity (corrected for inclination, v_rot = v_LOS/sin(i)) and σ₀ is the intrinsic velocity dispersion of the rotating disk, related to its thickness. In this work, we define σ₀ 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. C.1. From them, we infer an inclination i = 10 ± 7°, a v_rot = 360 ± 80 km s⁻¹, and a σ₀ = 190 ± 45 km s⁻¹ (at ∼0.3″, i.e. ∼2.4 kpc from the nucleus, as the more external regions are more affected by noise). The rotation-to-random motion ratio v_rot/σ₀ ≈ 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/σ₀ 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).

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⁻¹; 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⁴³ 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): M_dyn = (14 ± 6)×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 QSO-subtracted 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).

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α) = −0.3 ± 0.2 and log([O III]/Hβ) = 0.7 ± 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; Sect. 6.4.1).

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β. We used the equation

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{out}}(\mathrm{H}\beta )=8.6\times \frac{L_{41}(\mathrm{H}\beta )v_{\mathrm{out}}}{n_{\mathrm{e}}{R}_{\mathrm{out}}}{M}_{\odot }{\mathrm{yr}}^{-1}\end{array}$$(4)

from Cresci et al. (2015), where L₄₁(Hβ) is the Hβ luminosity associated with the outflow component in units of 10⁴¹ erg s⁻¹, 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⁻³, 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⁻³), as measured from recent JWST/NIRSpec IFS observations (Cresci et al. 2023).

Here we consider the L₄₁(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⁻¹]) = $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 (∼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_{\mathrm{out}}={930}_{-110}^{+60}$ km s⁻¹.

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⁴ M_⊙ yr⁻¹. 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 ∼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β clouds in the NLR of an AGN. The kinetic and momentum powers are ${\dot{E}}_{\mathrm{out}}=1/2{\dot{M}}_{\mathrm{out}}{v}_{\mathrm{out}}^2\sim 4\times {10}^{45}$ erg s⁻¹ and Ṗ_out = Ṁ_outv_out ∼ 8 × 10³⁷ 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).

6.3. 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⁵ M_⊙ yr⁻¹). 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α (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⁻¹ 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⁻¹. The velocity width (Moment 2) in these companions is significantly smaller than the ones in the QSO host.

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.

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.

Fig. 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σ). For the remaining sources, we considered 3 × 3 spaxel integration (Jil5, Jil6, Jil7, and Jil9), or 5 × 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.

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.

Jil1 is detected in [O III] at ∼3σ. 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α, and Hβ, but not in [N II] or [S II] lines. It shows a clumpy morphology, with an extension over ∼4 kpc in projection. A velocity gradient with amplitude of ∼ ± 100 km s⁻¹ 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. C.2. 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_rot/σ₀ ≈ 0.7, and question the presence of a rotating disk. Nevertheless, we infer a tentative dynamical mass for this system, log(M_dyn) = $(8_{-6}^{+16})\times {10}^9{M}_{\odot }$, considering a circular velocity of 75 km s⁻¹ (corrected for an inclination i = 80° ±15°), 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 ≈3× the stellar mass estimate inferred by Husemann et al. (2021, from spectral energy distribution analysis), but still consistent within the errors.

The presence of an obscured AGN in Jil2, initially proposed by Husemann et al. (2018a) 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β) ∼ 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α, and Hβ, but not in [N II] or [S II] lines. It shows an elongated morphology, with a clear and regular gradient with amplitude ±200 km s⁻¹ over ∼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. C.3. We infer a tentative dynamical mass log(M_dyn) = (1.3 ± 0.5)×10¹⁰ M_⊙, assuming a circular velocity of 70 km s⁻¹ (corrected for an inclination i = 80 ± 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β) ∼ 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 ∼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 ∼200 km s⁻¹ in 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β) > 0.7, consistent with flux ratios measured in the QSO host, and hence likely ionised by the QSO radiation.

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.

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⁻¹), 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⁻¹ from Jil3, and narrow line profiles (W80 ∼ 170 km s⁻¹). In this case, we detect a lower limit for log([O III]/Hβ) > 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α 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. 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 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 ∼ 2.

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

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

For Jil2, we also report a ∼4σ detection of the [S II] doublet, and hence a log([S II]/Hα) = −0.52 ± 0.07. This value places Jil2 in the Seyfert-like region of the line ratio diagnostic diagram [O III]/Hβ versus [S II]/Hα (Veilleux & Osterbrock 1987).

We infer for Jil2 an AGN bolometric luminosity log(L_bol/[erg s⁻¹]) ∼ 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. (2018a, 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.

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 ∼20 kpc, all within a velocity range of ∼ ± 250 km s⁻¹ 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⁴¹ erg s⁻¹) per Mpc⁻³ (Khostovan et al. 2015; Hirschmann et al. 2023); this corresponds to 5 × 10⁻⁴ 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.

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α 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 (∼0.8 kpc in a FOV of 25 × 25 kpc² 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.

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 ∼3″ × 3″ (23 × 23 kpc²), 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⁹ M_⊙. With a bolometric luminosity of log(L_bol/[erg s⁻¹]) ∼ 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/σ₀ ≈ 2, as inferred from 3D-BAROLO modelling. We also derived a tentative dynamical mass for the host, M_dyn = (14 ± 6)×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⁻¹ and a mass rate Ṁ_out ∼ 10⁴ M_⊙ yr⁻¹. 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, 2021), for a total of nine companions within 30 kpc of the QSO. All of these companions are within ±250 km s⁻¹ 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¹⁰ 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 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.

Acknowledgments

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. M.P., S.A., and B.R.P. acknowledge support from the research project PID2021-127718NB-I00 of the Spanish Ministry of Science and Innovation/State Agency of Research (MICIN/AEI). M.P. also acknowledges support from the Programa Atracción de Talento de la Comunidad de Madrid via grant 2018-T2/TIC-11715. M.A.M. 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. R.M., J.S. and F.D.E. acknowledge support by the Science and Technology Facilities Council (STFC), from the ERC Advanced Grant 695671 “QUENCH”. R.M. and J.S. also acknowledge funding from a research professorship from the Royal Society. G.C. 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. A.J.B., G.C.J. and A.J.C. 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). S.C. acknowledges support from the European Union (ERC, WINGS, 101040227). P.G.P.-G. acknowledges support from Spanish Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033 through grant PGC2018-093499-B-I00. I.L. 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 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

Appendix A: Nuclear spectra

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

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.

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. B.3). Moreover, in the very innermost nuclear regions, the Hα BLR emission covers a significant number of wavelength channels, preventing a proper modellisation of the underlying wiggles in the vicinity of the Hα line (see Fig. 6). This can affect the reconstruction of the Hα kinematics.

Fig. B.1. 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.

Fig. B.2. 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.

Fig. B.3. Same as Fig. B.2, but for the drizzle data cube, with spaxels of 0.05″.

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]λλ4959, 5007 line ratio is preserved at 1:3, which is consistent with theory (Osterbrock & Ferland 2006). Figure B.5 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]λ4959 peaks at ∼0.33, consistent with the expectations. We note however that significant deviations (up to ∼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. B.6.

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. B.5, off-centred integrated spectra are always affected by these wiggles.

Fig. B.4. 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″.

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.

Fig. B.6. 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 C.1, C.2, and C.3 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.

Fig. C.1. 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−1 in Fig. 9). The black and green lines identify the major axis and the zero-velocity curve, respectively; the red cross identifies the QSO position.

Fig. C.2. Jil2 Moment 0, 1, and 2 and the 3D-BAROLO disk kinematic best fit. See Fig. C.1 for details.

Fig. C.3. Jil3 Moment 0, 1, and 2 and the 3D-BAROLO disk kinematic best fit. See Fig. C.1 for details.

All Tables

All Figures

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.

In the text

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.

In the text

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

In the text

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.

In the text

Fig. 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 Fig. B.1 for analogous effects in the emsm cube.

In the text

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.

In the text

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.

In the text

	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.
In the text

	Fig. 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.
In the text

	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.
In the text

	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.
In the text

Fig. 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σ). For the remaining sources, we considered 3 × 3 spaxel integration (Jil5, Jil6, Jil7, and Jil9), or 5 × 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.

In the text

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.

In the text

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 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 ∼ 2.

In the text

	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.
In the text

Fig. B.1. 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.

In the text

Fig. B.2. 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.

In the text

	Fig. B.3. Same as Fig. B.2, but for the drizzle data cube, with spaxels of 0.05″.
In the text

	Fig. B.4. 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″.
In the text

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.

In the text

Fig. B.6. 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.

In the text

Fig. C.1. 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−1 in Fig. 9). The black and green lines identify the major axis and the zero-velocity curve, respectively; the red cross identifies the QSO position.

In the text

	Fig. C.2. Jil2 Moment 0, 1, and 2 and the 3D-BAROLO disk kinematic best fit. See Fig. C.1 for details.
In the text

	Fig. C.3. Jil3 Moment 0, 1, and 2 and the 3D-BAROLO disk kinematic best fit. See Fig. C.1 for details.
In the text