GA-NIFS: Black hole and host galaxy properties of two z$\simeq$6.8 quasars from the NIRSpec IFU

Integral Field Spectroscopy (IFS) with JWST NIRSpec will significantly improve our understanding of the first quasars, by providing spatially resolved, infrared spectroscopic capabilities which cover key rest-frame optical emission lines that have been previously unobservable. Here we present our results from the first two z>6 quasars observed as a part of the Galaxy Assembly with NIRSpec IFS (GA-NIFS) GTO program, DELS J0411-0907 at z=6.82 and VDES J0020-3653 at z=6.86. By observing the H$\beta$, [OIII], and H$\alpha$ emission lines in these high-z quasars for the first time, we measure accurate black hole masses, $M_{\rm{BH}}=1.85e9$ and $2.9e9$M$_\odot$, corresponding to Eddington ratios of $\lambda_{\rm{Edd}}=0.8$ and 0.4 for DELS J0411-0907 and VDES J0020-3653 respectively. These provide a key comparison for existing estimates from the more uncertain MgII line. We perform quasar-host decomposition using models of the quasars' broad lines to measure the underlying host galaxies. We also discover multiple emission line regions surrounding each of the host galaxies, which are likely companion galaxies undergoing mergers with these hosts. We measure the star formation rates, excitation mechanisms, and dynamical masses of the hosts and companions, measuring the $M_{\rm{BH}}/M_{\rm{dyn}}$ ratios at high-z using these estimators for the first time. DELS J0411-0907 and VDES J0020-3653 both lie above the local black hole--host mass relation, and are consistent with the existing observations of $z\gtrsim6$ quasar host galaxies with ALMA. We detect ionized outflows in [OIII] and H$\beta$ from both quasars, with mass outflow rates of 58 and 525 M$_{\odot}$/yr for DELS J0411-0907 and VDES J0020-3653, much larger than their host star formation rates of<33 and<54 M$_\odot$/yr. This work highlights the exceptional capabilities of the JWST NIRSpec IFU for observing quasars in the early Universe.


Introduction
Observational surveys over the last two decades have revealed several hundreds of high-quasars ( ≳ 6; e.g. Fan et al. 2000Fan et al. , 2001Fan et al. , 2003Willott et al. 2009Willott et al. , 2010Kashikawa et al. 2015;Bañados et al. 2016Bañados et al. , 2018Banados et al. 2022;Matsuoka et al. 2018;Wang et al. 2019;Yang et al. 2023a), less than a billion years into the Universe's history. The black holes powering these quasars have been measured to be as massive as ∼ 10 10 M ⊙ (Wu Article number, page 1 of 27 arXiv:2302.04795v2 [astro-ph.GA] 2 Aug 2023 A&A proofs: manuscript no. manuscript et al. 2015), raising questions about the nature of the first black hole seeds and their extreme growth mechanisms. Theories to explain the growth of these first quasars generally require massive black hole seeds, sustained accretion at the Eddington limit, or phases of accretion at super-Eddington rates (e.g. Volonteri & Rees 2005;Volonteri et al. 2015;Sĳacki et al. 2009). Indeed, some quasars have been measured to have accretion rates above the Eddington limit (e.g. Bañados et al. 2018;Pons et al. 2019;Farina et al. 2022), although when considering the uncertain black hole mass measurements and large scatter in the adopted scaling relations, Farina et al. (2022) found no clear evidence for a significant population of super-Eddington high-quasars. Thus accurate measurements of the black hole mass and accretion rate for high-quasars provide vital constraints for these theories.
Current high-quasar black hole mass measurements rely on calibrations that can be quite uncertain. The H and H lines are generally considered the most reliable single-epoch indicators of black hole mass, as they are used in many reverberation mapping studies (e.g. Wandel et al. 1999;Kaspi et al. 2000). However, these reverberation mapping studies are performed for low-, low-luminosity active galactic nuclei (AGN). At higher-, and for luminous quasars, extrapolations are required, and so at higheven these best estimates are uncertain. More critically, for highquasars these lines are redshifted into infrared wavelengths > 2.5 m, beyond the range observable from the ground (for sources at ≳ 4 for H and ≳ 3 for H ). Thus black hole mass measurements for high-quasars instead have typically relied on the Mg ii and C iv lines. Reverberation mapping studies based on Mg ii and C iv are more uncertain (e.g. Trevese et al. 2014;Kaspi et al. 2021), and the resulting scaling relations required to obtain single-epoch black hole mass estimates are affected by luminosity-dependent biases and provide highly uncertain mass measurements, for example due to significant blueshifts caused by outflows (e.g. Shen et al. 2008;Shen & Kelly 2012;Peterson 2009;Coatman et al. 2016;Farina et al. 2022, and references within). With the James Webb Space Telescope (JWST; Gardner et al. 2006Gardner et al. , 2023McElwain et al. 2023) and the Near-Infrared Spectrograph (NIRSpec; Jakobsen et al. 2022;Böker et al. 2023), H is now observable up to = 9.8, and H up to = 7.0 (and beyond with the Mid-Infrared Instrument, MIRI), allowing us to more accurately measure the black hole masses and Eddington ratios of high-quasars with the same estimator at all redshifts (e.g. Eilers et al. 2022;Yang et al. 2023b). This will improve our understanding of the first black holes and their growth rates by answering questions such as: are these black holes truly so massive, and are they accreting above the Eddington limit?
The JWST not only allows us to access the near infrared spectral range with unprecedented sensitivity, but also provides high angular resolution and a stable point-spread function (PSF). With these improvements over HST, simulations predict that the host galaxies of these high-quasars will be detectable for the first time (Marshall et al. 2021), and indeed the first signs of host emission have been detected with the Near-Infrared Camera (NIRCam; Ding et al. 2022). Thus JWST and the NIRSpec Integral Field Unit (IFU; Böker et al. 2022) will permit the characterization of the spatially resolved spectral properties of the hosts, measuring their stellar properties, interstellar medium (ISM) structure and kinematics, excitation mechanisms, and stellar-based SFRs. This will significantly advance our understanding of quasars, their host galaxies, and the rapid growth of early supermassive black holes.
Based on JWST and NIRSpec IFU observations, in this paper we use the standard H and H single-epoch black hole mass estimators to derive accurate black hole mass measurements for two > 6 quasars, DELS J0411−0907 and VDES J0020−3653. DELS J0411−0907 (VHS J0411-0907) at = 6.82 was discovered independently by Pons et al. (2019) and Wang et al. (2019). DELS J0411−0907 is the most distant of the Pan-STARRS quasars, and is the brightest quasar in the J-band above = 6.7 (Pons et al. 2019). With a relatively low Mg ii black hole mass of (6.13 ± 0.51) × 10 8 M ⊙ , and a bolometric luminosity of (1.89 ± 0.07) × 10 47 erg s −1 , Pons et al. (2019) measured an Eddington ratio of 2.37 ± 0.22, implying that this quasar could be a highly super-Eddington quasar undergoing extreme accretion. The quasar VDES J0020−3653 at = 6.86 was discovered by Reed et al. (2019). VDES J0020−3653 has a larger Mg ii black hole mass than DELS J0411−0907, measured as (16.7 ± 3.2) × 10 8 M ⊙ by Reed et al. (2019), and rederived as 25.1 × 10 8 M ⊙ by Farina et al. (2022) using the Shen et al. (2011) Mg ii-black hole mass relation. With a bolometric luminosity of Bol = (1.35 ± 0.03) × 10 47 erg s −1 , these correspond to Eddington ratios of 0.62 ± 0.12 and 0.43 respectively, suggesting that this quasar is powered by a larger supermassive black hole with less extreme accretion than DELS J0411−0907. However, these black hole mass and accretion measurements are based on the Mg ii line; robust H and H mass estimates will better constrain the accretion of these quasars. This paper is outlined as follows. In Section 2 we describe the JWST NIRSpec observations, and our data reduction and analysis techniques. In Section 3 we present our results, with the black hole properties in Section 3.1, the host galaxy properties in Section 3.2, and the outflow properties in Section 3.3. We present a discussion in Section 4, before concluding with a summary of our findings in Section 5.
2.87-5.27 m. The observations were taken with a NRSIRS2 readout pattern with 25 groups, using a 6-point medium cycling dither pattern, resulting in a total exposure time of 11029 seconds. We constrained the complete position angle (PA) window available for observing these targets with JWST to a more limited PA range, to minimise the leakage of light through the microshutter assembly (MSA) from bright sources. Additionally, these quasars have been observed with NIRSpec fixed-slit spectroscopy at shorter wavelengths (with G140H /F070LP, G235H/F170LP), which will be analysed in a separate paper.
DELS J0411−0907 was observed on Oct 1, 2022, at a PA of 61.893 degrees. VDES J0020−3653 was observed on Oct 16, 2022, at a PA of 160.169 degrees. Unfortunately, only three out of the six planned dithers were observed due to a spacecraft guiding failure. The final three dithers were observed on Nov 17, 2022, at a position angle of 188.177. With half of the VDES J0020−3653 observations observed at a different PA, combining these exposures and their subsequent data analysis was more challenging than the DELS J0411−0907 observations. Due to an error in the reference files of the pipeline (see Perna et al. 2023, for details), we performed a bona fide astrometric registration matching the peak pixel locations of the H BLR collapsed images. We applied a constant offset in RA and Dec to the pipeline Stage 1 products of the second set of observations, which allowed the pipeline to correctly combine all six dithers in a final data cube.
Due to an uncertainty in the specified astrometry of our observations, we assume that the quasar is centred on the known positions from the literature: 04h 11m 28.63s, −09d07m49.8s for DELS J0411−0907 (Wang et al. 2019) and 00h 20m 31.472s, −36d53m41.82s for VDES J0020−3653 (Reed et al. 2019). These were the targeted positions for our IFU observations, which have been found to have a typical pointing accuracy of ∼ 0.1 ′′ compared to the requested pointing (Rigby et al. 2022).

Data Reduction
The data reduction was made with the JWST calibration pipeline version 1.8.2 (Bushouse et al. 2022), using the context file jwst_1023.pmap. All of the individual raw images were first processed for detector-level corrections using the Detec-tor1Pipeline module of the pipeline (Stage 1 hereinafter). The individual products (count-rate images) were then calibrated through Calwebb_spec2 (Stage 2 hereinafter), where WCScalibration, flat-fielding, and the flux-calibrations are applied to convert the data from units of count rate to flux density. The individual Stage 2 images were then resampled and coadded onto a final data cube through the Calwebb_spec3 processing (Stage 3 hereinafter). A number of additional steps and corrections in the pipeline code have been applied to improve the data reduction quality. In particular, -The individual count-rate frames were further processed at the end of Stage 1, to correct for a different zero-level in the dithered frames. For each image, we therefore subtracted the median value, computed considering the entire image, to get a base-level consistent with zero counts per second. -We further processed these count-rates to subtract the 1/ noise (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 sigmaclipping. The modelled 1/ noise is then subtracted before proceeding with the Stage 2 of the pipeline. -The outlier_detection step of the third and last step of the pipeline is required to identify and flag all remaining cosmic-rays and other artefacts that result in a significant number of spikes in the reduced data. In pipeline version 1.8.2, this step tends to identify too many false positives that compromises the data quality significantly. We therefore implement a modified version of the lacosmic algorithm (van Dokkum 2001) in the JWST pipeline, which is applied at the end of Stage 2 (details in D' Eugenio et al. in prep). -We applied the cube_build step to produce combined data cubes with a spaxel size of 0.1 ′′ , using the Exponential Modified Shepard Method (EMSM) weighting. This provides high signal-to-noise at the spaxel level yet slightly lower spatial resolution with respect to other weightings; the loss in resolution is compensated by the fact that these cubes are less affected by oscillations in the extracted spectra due to resampling effects, for example (see Law et al. 2023;Perna et al. 2023, and Appendix C).
We modified the cube_build script to fix an error affecting the drizzle algorithm, as implemented in the pipeline version 1.8.5 which was not released at the time of data reduction1. Finally, we also modified the photom script, applying the corrections implemented in the same unreleased pipeline version2, which allows us to infer more reasonable flux densities (i.e. a factor ∼ 100 smaller with respect to those obtained with standard pipeline 1.8.2). To test this calibration we also reduced the data using an external flux calibration (details in Perna et al. 2023), and obtained consistent results.
A complete description of the data reduction process and the tests we have performed to ensure its robustness can be found in Perna et al. (2023).

Quasar and Host Integrated Line Fitting
To obtain our integrated quasar spectrum, we integrate our data cube across an aperture of radii 0.35" for the wavelengths blueward of the detector gap, and 0.45" red-ward of the gap. These apertures were chosen by examining the PSF structure around H and H . These aperture radii contain the central peak of the PSF, but not the ring of reduced flux that occurs prior to the second radial PSF peak. This allows us to maximise our signal to noise, while containing the majority of the quasar flux.
We compare the flux in this aperture to the total flux in the image (to radius 1.5") at the peak H and H wavelengths, finding an average fraction for both quasars of 81.4% at H , and 86.6% at H . We use these percentages as our flux correction, multiplying the integrated aperture flux blue-ward of the detector gap by 1.228 and the integrated flux red-ward of the detector gap by 1.155.
To fit our quasar spectra we use the Markov chain Monte Carlo (MCMC)-based technique of QubeSpec3. Our quasar spectra contain H and H emission from the BLR, which are blended with the narrow galaxy H , [O iii] 4959,5007,and H lines (galaxy narrow: 'GN'), and also with broader lines from the galaxy associated with a kinematically distinct component (e.g. an outflow; galaxy broad: 'GB'). We note that we see no evidence for any [N ii] 6548, 6583 emission in our cubes, from either the quasar or the host galaxies and companions, and so we do For the H and H lines, we find that the BLR signal cannot be well fit with a single Gaussian. We instead fit the BLR H and H lines with a broken power law, that is convolved with a Gaussian with standard deviation (as in e.g. Nagao et al. 2006;Cresci et al. 2015). In our spectral fit we assume that the galaxy lines are composed of a narrow component (GN), as well as a broad component (GB), which is typically blue-shifted and therefore generates a wing in the total profile, which could be attributed to ionized gas outflows. and 2 , the initial guess is the respective H fit property, with normally distributed priors surrounding that value. This is again because the lines likely arise from the same physical region with similar kinematics. We also choose an initial guess for the H GN, GB, and BLR peak values to be 3.5 times that of the respective H values, with normally distributed priors around that value, which helps the fit to converge on a H model that has realistic flux ratios relative to H . Despite the H -[O iii] 4959, 5007 and the H complexes being fit independently, the agreement among the various kinematic properties was excellent, particularly for the host components (see Table 1). However, for VDES J0020−3653 we note that the H BLR full width at half maximum (FWHM) is 3092 +39 −1 km/s, which is significantly lower than the H FWHM of 4461 +233 −182 km/s. This discrepancy may be due to the H line being inaccurately measured as it lies close to the edge of the wavelength coverage of the detector; we are more confident in the H line measurement. There are also fit degeneracies which could be affecting this measurement.
We show the full integrated quasar spectrum for DELS J0411−0907 and VDES J0020−3653 in Figure 1. In Figure 2 we show the H and H regions of these integrated spectra, alongside the best-fitting model. The velocity offset and FWHM for each of the model components are given in Table 1, and the resulting BLR line fluxes and luminosities are given in Table 2.

Continuum Subtraction
To further analyse our quasar data cubes, we first subtract the continuum emission from the spectra in each spaxel, by considering the flux in the continuum windows at rest frame 3790-3810Å, 4210-4230Å, 5080-5100Å, 5600-5630Å and 5970-6000Å, following Kuraszkiewicz et al. (2002) and Kovacevic et al. (2010); these windows are selected to be uncontaminated by emission lines and blended iron emission. We also added an additional continuum window at the red side of H at the edge of the spectral range, at 5.260-5.266 m in the observed frame, as well as a window near the blue edge of the spectral range, at 2.975-2.990 m in the observed frame, to better constrain the continuum. We calculate the signal-to-noise ratio (SNR) at each wavelength in each spaxel, taking the noise to be the 'ERR' cube produced by the pipeline. For each spaxel in which the median SNR within these continuum windows is greater than 2, we interpolate between these continuum windows using scipy's interp1d function using a spline interpolation of second order; this is found to be more robust than fitting a power law curve in each spaxel. In each of these spaxels the interpolation is subtracted from the spectra to create our continuum-subtracted cube. Figure 1 shows the original and continuum-subtracted spectrum of the brightest quasar spaxel.
We note that for VDES J0020−3653 there may be some broad H flux in this reddest continuum window in the brightest quasar spaxels, as at this redshift the IFU spectral coverage does not extend to the pure continuum region red-ward of H . However, we find that this window is necessary for finding a reasonable fit to the continuum around H . To reduce any potential effect of this on our results, we use the full, non-continuum subtracted cube in Section 2.3.1 when we fit our integrated quasar model to investigate the broad line properties. With our QubeSpec modelling we fit the emission lines as well as the continuum in the form of a power law. As the integrated quasar spectrum is well described by a power law, this approach does not require the continuum window strategy used here. In contrast, within individual spaxels the continuum is not well described by a power law due to the spatial variation of the PSF with wavelength, and so this continuum window approach is more reliable. For our spatially resolved analysis we are primarily interested in the narrow line emission, which is further from the edge of the detector, and so this should be less biased by this issue. However, we do caution that at this high redshift, the H measurement of VDES J0020−3653 could be generally susceptible to edge effects from the detector and an imperfect continuum subtraction.

Quasar Subtraction
In order to fully understand our quasar-host systems, we decompose the observed emission into its various components. The aim is to subtract the bright unresolved flux from the quasar, revealing the extended narrow line emission from the host galaxy.
The full cube contains broad line emission from the quasar's broad-line region (BLR) surrounding the black hole, for example from the H and H lines. The BLR is spatially unresolved at these redshifts, and thus this broad line emission is spread spatially according to the PSF of the instrument. In our BLR Quasar Subtraction method (see below), we subtract this BLR component from the cube. The remaining flux in the cube is narrow line emission, for example from [O iii] 4959, 5007 as well as the H and H lines which also have a strong broad component. This narrow emission is emitted both by the quasar's narrow line region (NLR), photoionized by the AGN and shocks, as well as gas throughout the extended host galaxy (and any companion galaxies) which is photoionized by star formation. While the NLR can be spatially extended, significant NLR emission is produced by the spatially unresolved central region surrounding the AGN, and so this component will clearly trace out the PSF shape. To understand the true spatial extent and kinematics of the extended narrow line emission, and not be biased by the bright unresolved NLR emission, we subtract this narrow unresolved flux from our cube. This is our Narrow & Broad (N&B) Quasar Subtraction, as detailed below. To separate the host and quasar emission, we use the code QDeblend3D (Husemann et al. 2013(Husemann et al. , 2014, which uses the relative strength of quasar broad lines in each spaxel to map out the spatial PSF. This strategy allows us to use the specific quasar observation to create an accurate PSF model, instead of relying on an external observation of a known point source (i.e. star) or a theoretical PSF model. The spatial PSF varies with the precise position of the source on the detectors, the general position on the sky, and time, so matching the PSF of our specific observations with a PSF observation or model is less precise than the broad line model we can create here. This methodology has successfully been used for IFU observations of lower-redshift quasar host galaxies (e.g. Husemann et al. 2013Husemann et al. , 2014Vayner et al. 2016;Kakkad et al. 2020).
As the NIRSpec PSF varies with wavelength, we opt to perform the quasar subtraction at two different wavelength regimes. For the spectra blueward of the detector gap, we use the H line for the subtraction, and for the spectra redward of the detector gap we use the H line for subtraction. This allows for the best subtraction for the wavelength region around each of these two lines.
Step 1, Spatial PSF model: We determine the spatial PSF shape by measuring the relative flux of the quasar BLR at each spatial location. The BLR is spatially unresolved and thus this emission will trace the PSF shape of the instrument. To measure this relative flux we use the continuum-subtracted cube. For H , we take the mean of the flux between rest-frame 6500-6523Å and 6592-6615Å, either side of the peak of H . These broad spectral windows are free from any narrow lines, which would bias the measurement; they are marked in Figure 3, which shows an example of this quasar subtraction process. For H , the equivalent is performed for emission between rest frame 4800-4830Å and 4880-4910Å. Once this broad line flux is measured in each spatial pixel, it is normalised to that of the central brightest quasar spaxel, giving a 2D fractional map of the relative flux of the spatially unresolved BLR, that is, the 2D PSF.
To identify spaxels with significant quasar flux from the PSF model, we calculate the SNR in each individual wavelength channel across the selected broad spectral window for each spaxel, using the 'ERR' extension as the noise. For our final PSF shape, we select spaxels in which at least three consecutive wavelength channels have SNR > 10, and then expand this region spatially by selecting any adjacent spaxels with SNR > 2 over two consecutive wavelength channels, following the find_signal_in_cube algorithm of Sun (2020) (see also Sun et al. 2018). To reduce the effect of any artefacts and companion galaxies on the measured PSF shape, we cut the PSF at 10 , where = FWHM/(2 √︁ 2 log 2) and the FWHM is the diffraction limit of the telescope at the observed wavelength of H and H , 0.146" and 0.201" respectively. This ensures we are only performing the quasar subtraction in regions that have a broad line signal, while ensuring we extract as much of the quasar as possible. The resulting masked PSF shapes are shown in Figure  A.1, for both the H and H subtractions for both quasars. Measuring the FWHM of the resulting H PSFs gives 0.18" for VDES J0020−3653 and 0.16" for DELS J0411−0907, and 0.23" for H for both VDES J0020−3653 and DELS J0411−0907, consistent with expectations for the PSF size.
Step 2, Quasar spectral model: To perform the quasar subtraction we require a 3D quasar model cube to subtract from the initial quasar and host cube. To create a 3D cube for the quasar emission, alongside the PSF shape (spatial x,y dimensions) we also need a model for the quasar spectrum (spectral dimension), for which we consider two approaches. The first approach is to use the observed spectrum of the brightest spaxel, the centre of the quasar emission, under the assumption that the host line emission is negligible. We refer to this as the 'N&B Quasar Subtraction', as this quasar cube includes both the narrow (N) and broad (B) components of the central spectrum, and so both components will be subtracted. Our second approach is to use a model of the BLR component of the quasar spectrum. We take this model from our fit to the integrated quasar spectrum (the red curve in Figure 2, with values quoted in Table 1). This quasar model only includes the BLR component of the fit, and so we refer to this as the 'BLR Quasar Subtraction' as only the BLR emission is subtracted from the initial cube. Our process for integrating and fitting the quasar spectrum is described in detail in Section 2.3.1.
In our quasar subtraction, the spatial variation of the quasar's BLR flux is used to determine the PSF shape. However, the N&B subtraction subtracts both the broad and narrow emission that are present in the central quasar spaxel. These emission lines are emitted from the central, unresolved NLR, which is spread according to the PSF shape. This N&B subtraction results in a clear view of the narrow line flux in the underlying, spatially extended host galaxy, without contamination from the bright lines from the central NLR region that have been smeared by the PSF shape. This is ideal for measurements of the spatial distribution and kinematics of emission line regions, as well as measurements of the SFRs which would otherwise be significantly biased by the large amount of AGN-photoionized flux in the centre. However, this method uses the spectrum from the individual brightest quasar spaxel, and so can be subject to observational artefacts that would be removed when integrating over a larger region, for example oscillations in the spectra (see Perna et al. 2023). We chose to use the individual spaxel spectrum and not an integrated spectrum for the model as we found that this resulted in the best quasar subtraction. This model also assumes that the host contribution to the line flux in this central spaxel is negligible, while there may be some small contribution at the position of the optical lines analysed in this work.  , showing the spectrum from the brightest spaxel (navy), alongside the integrated spectrum (orange). The integrated spectrum is integrated over a radius of 0.35" for the wavelengths blue-ward of the detector gap, and 0.45" red-ward of the gap, corresponding to the apertures selected for our quasar fitting (Section 2.3.1). Our flux corrections of 81.4% and 86.6% have been applied to the integrated spectrum. The spectrum of the brightest spaxel is increased by a constant scaling factor, 8.230 for DELS J0411−0907 and 9.354 for VDES J0020−3653, such that the peak of H is equal for both spectra in the plots. The lower panel in each shows the spectra after continuum subtraction. For the brightest spaxel spectrum, the continuum model subtracted from the spectrum is shown in black in the top panel. For the continuum-subtracted integrated spectrum in the lower panel, the integration is performed on the continuum-subtracted cube (i.e. the continuum subtraction was first performed on each of the spaxels individually). The dotted grey vertical lines depict lines of interest at the redshift of the quasar ( = 6.854 for VDES J0020−3653 and = 6.818 for DELS J0411−0907), and the grey regions show the continuum windows. The detector gap from 4.06-4.17 m has been masked. We note that some minor oscillations are present in the single spaxel spectrum (navy), which can be particularly noticed below 3.3 m in the single spaxel spectrum. These are caused by spatial undersampling of the PSF, which is a known issue when considering spectra in individual NIRSpec IFU spaxels; we describe this in detail in Perna et al. (2023). When measuring the quasar properties in Section 3.1, we use the integrated spectrum (orange), thus removing these single pixel effects. For the quasar subtraction, we find the best performance using the single spaxel spectrum (navy), however these oscillations do not significantly affect our results which directly focus on the regions around H and H .
In contrast, the BLR quasar subtraction method allows us to view the complete narrow line emission maps, including the contribution of the central unresolved NLR. Our BLR subtraction technique subtracts only the BLR quasar flux, with all of the narrow flux from the central core remaining. Much of this remaining narrow line flux is spread according to the PSF shape, which significantly contaminates any spatial analysis of the extended host region. This approach is also model dependent, varying based on the BLR model fit to the integrated quasar spectrum. With their various advantages and disadvantages, we use both approaches throughout this work.
Step 3, Subtracting the 3D PSF model: For each of these choices for the model quasar spectrum, we scale the quasar spectrum by the 2D PSF to create a 3D quasar cube. We then subtract this quasar cube from the original cube, revealing our estimate of the host galaxy's 3D spectrum. We visually inspect the quasar model fit and residual for a range of spaxels to confirm our subtraction is robust; an example for one spaxel is shown in Figure  3.

Host Line Maps
To estimate the noise for our host galaxy maps, we consider that the quasar subtraction introduces additional noise to the resulting host cube that should be accounted for. Following a Monte Carlo process, we generate 200 different realisations of the continuum-subtracted cube with normally-distributed noise, where is taken from the 'ERR' extension from the original cube. We then perform the H and H quasar subtraction on each of these 200 cube realizations. To create our modified noise cube, we calculate the standard deviation of these 200 realisations of the quasar subtracted cube at each wavelength in each spaxel.
For the H -[O iii] 4959, 5007 line complex the values in these noise maps are up to 34% higher than the original noise maps in the central spatial regions of the quasar for the N&B quasar subtraction, and 4% higher for the BLR subtraction. For the H line, the values in the noise map are up to 47% higher than the original noise map for the N&B quasar subtraction, and 10% higher for the BLR subtraction. The N&B subtraction introduces more noise than the BLR subtraction, as the N&B quasar model spectrum contains noise while the BLR quasar spectrum is a smooth model and the model uncertainties are not considered in this calculation. The quasar subtraction adds additional uncertainty, particularly where more quasar flux has to be subtracted.
To create maps of the H , H , and [O iii] 4959, 5007 emission of the host galaxy, we fit the spectrum in each spaxel as a series of Gaussians, using AstroPy's Levenberg-Marquardt algorithm. We assume that the narrower galaxy lines are com-Article number, page 7 of 27 A&A proofs: manuscript no. manuscript  This shows the spectrum in spaxel (23,24), with the quasar peak located at spaxel (26, 24), in image coordinates. The orange spectrum is the original continuum-subtracted spectrum in the spaxel, containing line flux from both the quasar and host galaxy. The black line is the quasar model spectrum. The top panels show the choice that the quasar model spectrum is that of the central brightest quasar spaxel, which contains both narrow and broad emission-the N&B subtraction. The bottom panel shows the choice that the quasar model spectrum is the BLR component from our best-fit model to the quasar spectrum ( Figure 2)-the BLR subtraction. Using the broad spectral windows marked in grey, the quasar model (black) is scaled to fit the combined spectrum (orange). The purple spectrum shows the residual to that fit, which is the host spectrum. posed of a narrow component (GN) as well as a broader wing (GB; see Section 2.3.1). For the GN and GB components, we fit each of the four lines as a Gaussian with mean velocity and line widths equal in velocity space. We tie the amplitude of [O iii] 4959 to be 3 times less than that of [O iii] 5007, for both the GN and GB components. These theoretical flux ratio constraints provide a good fit for the host galaxy [O iii] 4959, 5007 lines for both the N&B and BLR quasar subtraction, which gives further confidence in our subtraction techniques. As noted above, we see no evidence for any [N ii] 6548, 6583 emission in our cubes, and do not include those lines in our fit. When the spectrum can be well fit with two Gaussians, the Gaussian with the smallest velocity dispersion is assigned as the GN component, while the largest is assigned to the GB component. When only one Gaussian is preferred, if it has ≤ 33Å, or approximately 600 km/s, it is assigned as a GN component, otherwise as a GB component.
To determine the flux of each of the GN and GB line components we integrate the individual Gaussian model. The SNR is calculated by integrating the noise cube, calculated from the quasar subtraction algorithm, across the same spectral region.
Using these noise maps we search for spaxels with SNR > 5, and then expand this region spatially by selecting any adjacent spaxels with SNR > 1.5, following the find_signal_in_cube algorithm of Sun (2020). In the maps we only display spaxels meeting this criterion. The resulting maps are shown in Figures  4 and 6 for the N&B subtraction for DELS J0411−0907 and VDES J0020−3653 respectively, with the BLR subtraction maps for both quasars shown in Figure B.1 to show the differences between the two subtraction methods.
To create velocity maps, we calculate the median velocity 50 , a non-parametric measure for the line centre. For the velocity dispersion maps, we calculate 80 , the width containing 80% of the line flux. This is a non-parametric measure, which is approximately equal to the FWHM for a Gaussian line ( 80 = 1.088 FWHM; see e.g. Zakamska & Greene 2014). The velocity dispersion maps are derived after removing the instrumental width in quadrature at each spaxel, where for NIRSpec in the G395H/F290LP configuration, FWHM inst,H ≃ 115 km/s, and FWHM inst,H ≃ 85 km/s for the redshift of these quasars (Jakobsen et al. 2022). These velocity and velocity dispersion Table 2. The fluxes and luminosities measured for the 5100Å continuum, H and H BLR lines, alongside the derived black hole masses and Eddington ratios using these properties. The errors for the black hole masses and Eddington ratios include the fit uncertainties as well as the 0.43 dex uncertainty due to the intrinsic scatter in the scaling relations. We note that all flux uncertainties are measurement uncertainties and do not include any intrinsic uncertainty in the flux level from the data calibration.
maps are also shown in Figures 4 and 6 for the N&B subtraction, with the BLR subtraction maps also presented in Figure B.1.

Redshift
We show the integrated quasar spectrum for DELS J0411−0907 and VDES J0020−3653 in Figure 2, with the parameters from our spectral fit given in Table 1. We derive the quasar redshift based on the galaxy narrow (GN) [O iii] 5007 emission. For DELS J0411−0907, we find a redshift of = 6.818 ± 0.001. This is broadly in agreement with the Lyman-redshift of DELS J0411−0907, Ly = 6.81 ± 0.03 (Pons et al. 2019;Wang et al. 2019), and the Mg ii redshift, Mgii = 6.824 (Pons et al. 2019

Black Hole Masses
Black hole masses can be estimated using the virial relation: where is the size of the BLR, Δ is the width of the emission line, is the gravitational constant and is a scale factor. To robustly determine the BLR size , reverberation mapping studies measure the time delay of variability in the emission line to the continuum region, with = , where is the speed of light. Empirical correlations have been found between and the continuum luminosity of the quasar (e.g. Kaspi et al. 2000), with ∝ , allowing the black hole mass to be estimated from single-epoch measurements.
However, the continuum luminosity can be difficult to measure due to contamination by Fe ii lines and, for lower luminosity quasars, the host galaxy itself. To circumvent this, Greene & Ho (2005) and Vestergaard & Peterson (2006) Using our measured H and H broad line properties, we calculate the black hole masses using Equations 4 and 5 for both the Greene & Ho (2005) and Vestergaard & Peterson (2006) correlations. These are also listed in Table 2. As these empirical black hole mass relations have large scatter, we take the error in these relations to be 0.43 dex, the quoted scatter from Vestergaard & Peterson (2006), and add this in quadrature with our measurement uncertainty from the MCMC fit.

DELS J0411−0907:
For DELS J0411−0907, we find that the various scaling relations (Equations 3-5) give a black hole mass between 0.9 and 2.5 × 10 9 M ⊙ (Table 2). These measures are all consistent within the uncertainties. We note that while the H -L 5100 black hole mass estimated by the Greene & Ho (2005) and Vestergaard & Peterson (2006) relations are similar, at 2.0 and 1.7×10 9 M ⊙ respectively, the pure H black hole mass estimates are quite different between the two versions of the relation. The Greene & Ho (2005) H relation gives a low black hole mass of 1.3 × 10 9 M ⊙ , while the Vestergaard & Peterson (2006) relation gives 2.5 × 10 9 M ⊙ , twice as large. This variation between the two scaling relations could be due to the brightness of this quasar, which is brighter than the quasars considered in those studies. We also note that the Greene & Ho (2005) H black hole mass for DELS J0411−0907 is 0.9 × 10 9 M ⊙ , the lowest of these estimates, yet this is consistent with the other estimates within the uncertainties. As the pure H and H relations involve a further calibration from 5100 to H /H , we opt to take the average of the two BH,H ,5100 measurements as our best black hole mass estimate: 1.85 +2 −0.8 ×10 9 M ⊙ . We reiterate that all of these black hole mass estimates are consistent with each other within the uncertainties.
In comparison, from the FWHM of the Mg ii line and the restframe 3000Å luminosity measured from Magellan FIRE spectra (Pons et al. 2019), the DELS J0411−0907 BH mass derived using the Vestergaard & Osmer (2009) relation is BH,Mgii = (0.61 ± 0.05) × 10 9 M ⊙ . This is much lower than our measurements. This quoted error considers only measurement uncertainties, and does not include the additional ∼ 0.55 dex uncertainty from the scatter in the scaling relation (Vestergaard & Osmer 2009), which would be BH,Mgii = 0.6 +1.6 −0.4 × 10 9 M ⊙ . Alternatively, applying the Shen et al. (2011)  −0.7 × 10 9 M ⊙ . While these Mg ii black hole masses are all lower than our best estimates, when including the large uncertainty from the scatter in the scaling relations, all of the measurements are consistent.
In summary, the best estimate for the black hole mass of DELS J0411−0907 from NIRSpec/IFU, BH = 1.85 +2 −0.8 × 10 9 M ⊙ , is larger than the existing measurements from Mg ii, BH,Mgii = 0.6-1 × 10 9 M ⊙ . Our measurements use the H -5100Å black hole mass relation with scatter of ∼ 0.43 dex, compared with the larger ∼ 0.55 dex scatter in the Mg ii relation, highlighting the benefit of using these Balmer lines as black hole mass estimators.

VDES J0020−3653:
For VDES J0020−3653, we find that the various scaling relations give a wide range of black hole mass estimates, between 1.1 and 4.5 × 10 9 M ⊙ . These estimates are all consistent within the uncertainties. We note that the lowest measurement is the H mass measurement of 1.1 × 10 9 M ⊙ . For VDES J0020−3653, the H line falls at the edge of the wavelength coverage of the detector, which could cause issues with this measurement. Thus the H measurements of 2.3-4.5 × 10 9 M ⊙ are likely more reliable, although all measurements are consistent. As the pure H and H relations involve a further calibration from 5100 to H /H , we opt to take the average of the two BH,H ,5100 measurements as our best black hole mass estimate: 2.9 +3.5 −1.3 × 10 9 M ⊙ . We note again that the Greene & Ho (2005) and Vestergaard & Peterson (2006) black hole masses are quite different, with the Greene & Ho (2005) estimate larger for the H -L 5100 relation, while the Vestergaard & Peterson (2006) mass is twice as large for the pure H relation. This variation could again be due to the brightness of this quasar, which is brighter than the quasars considered in those studies. We reiterate that all of these black hole mass estimates are consistent within the uncertainties.
In comparison, Reed et al. (2019) measured a BH mass from the Mg ii line and the rest-frame 3000Å luminosity of BH,Mgii = (1.67 ± 0.32) × 10 9 M ⊙ , using the Vestergaard & Osmer (2009) relation. Accounting for the additional ∼ 0.55 dex scatter in the Mg ii black hole mass relation, this is BH,Mgii = 1.7 4.3 −1.2 × 10 9 M ⊙ . Applying the Shen et al. (2011) relation to the same Mg ii line FWHM and rest-frame 3000Å luminosity from Reed et al. (2019), Farina et al. (2022) derived a larger BH,Mgii = 2.51 × 10 9 M ⊙ . These estimates are in agreement with our estimated black hole mass of 2.9 +3.5 −1.3 × 10 9 M ⊙ . For our two quasars, we find that VDES J0020−3653 has a H -L 5100 black hole mass that is consistent with the existing Mg ii black hole mass estimates, while the H -L 5100 black hole mass for DELS J0411−0907 is larger than the Mg ii estimates, albeit consistent within the large uncertainties. Shen et al. (2008) found that for ∼ 8000 low-SDSS quasars, log( BH,H / BH,MgII ) follows a Gaussian distribution with mean of 0.034 and dispersion of 0.22 dex, generally indicating a good agreement yet with some quasars found to have significantly different values of log( BH,H / BH,MgII ) up to ±1.5. While high-quasars may follow a different distribution, our quasars are consistent with this Shen et al. (2008) scenario where Mg ii and H mass estimates can be consistent, although for some quasars these estimates can vary significantly. This highlights the need for obtaining H and H masses of a larger sample of highquasars with JWST, to more accurately measure black hole mass in the early Universe (see also Yang et al. 2023a).

Eddington Ratios
We now consider the Eddington Ratio Edd = Bol / Edd , the ratio of the quasars' bolometric to Eddington luminosity.
The quasars' bolometric luminosities have been estimated from their rest-frame 3000Å luminosities using the conver- The Eddington luminosity, the theoretical maximum luminosity an object could have while balancing radiation pressure and gravity, is typically assumed to be where is the gravitational constant, the proton mass, the speed of light, and the Thomson scattering cross-section. Using the various black hole mass estimates from Section 3.1.2, we calculate the corresponding Eddington ratios which are listed in Table 2.
From our best estimate of the black hole mass for DELS J0411−0907, BH = 1.85 +2 −0.8 × 10 9 M ⊙ , we estimate an Eddington ratio of Edd = 0.8 +0.7 −0.4 , with our full range of estimates from 0.6-1.6. In comparison, the Pons et al. (2019)  VDES J0020−3653 has a lower bolometric luminosity than DELS J0411−0907, and a larger black hole mass, resulting in a lower Eddington ratio for this quasar. From our best estimate of the black hole mass for VDES J0020−3653, BH = 2.9 +3.5 −1.3 × 10 9 M ⊙ , we estimate an Eddington ratio of Edd = 0.4 +0.3 −0.2 , with our full range of estimates from 0.2-0.5, except for the large H estimate of 0.9 +0.6 −0.4 , which could be affected by its proximity to the edge of the detector. In comparison, the estimate of Edd based on the Mg ii BH mass was Edd,Mg ii = 0.62 ± 0.12 for VDES J0020−3653 (Reed et al. 2019), or alternatively Edd,Mg ii = 0.43 as re-calculated by Farina et al. (2022). Our derived Eddington ratios are consistent with both values, which all suggest a sub-Eddington accretion rate for this quasar.

ISM Structure and Kinematics
In this section we study the spatial structure of the emission line regions present within our data cubes. To verify the spatial structures that we see, we produce a second set of data cubes by sub-sampling to a pixel scale of 0.05"/pix, and combining the various dithers with a different weighting scheme (see Appendix C), giving higher spatial resolution data cubes. Unfortunately, due to PSF effects the quasar subtraction does not produce successful results using these cubes. As a result, we only use these cubes to qualitatively examine the spatial morphologies in greater detail, with all quantitative results instead using the 0.1"/pix cubes. We show the [O iii] 5007 maps, as measured from these 0.05"/pix data cubes in Figure C.1, and discuss this further in Appendix C.

DELS J0411−0907:
Our quasar-subtracted cube for DELS J0411−0907 reveals three extended narrow-line emission regions, as shown in Figure  5 and labelled as Regions 1, 2, and 3. These three structures can be seen in both the 0.1"/pix and 0.05"/pix cubes. The integrated spectra for each of the regions are also shown in Figure 5. Our kinematic maps for the subtracted cube, showing the flux of the lines, their relative velocity, and linewidth in each spaxel are shown in Figure 4 for the N&B quasar subtraction and in Figure B.1 for the BLR quasar subtraction. We fit the emission lines in the integrated spectra of each of the three regions with either one or two Gaussians, depending on whether a broader wing component is seen, using the spectral fitting method used to construct our kinematic maps described in Section 2.3.4. For Region 2, we exclude the central 9 pixels surrounding the quasar peak, as these are highly corrupted by the quasar subtraction and introduce significant noise and artefacts. The extracted flux values and velocities are given in Table 3 and the flux ratios and our derived physical properties are listed in Table 4.
Region 1 is a companion galaxy to the southeast of the quasar. The ellipse defining this region in Figure 5 is offset 0.78 ′′ from the quasar position. In WMAP9 cosmology at = 6.818, 1 ′′ = 5.446 kpc, and so this offset is 4.26 kpc, quasar-to-centre. For reference, the region ellipse, which we fit by eye, has = 0.814 ′′ = 4.4 kpc, = 0.341 ′′ = 1.9 kpc and a PA of 54.6 degrees. By fitting a 2D Gaussian to the spatial distribution of the [O iii] 5007 flux in this region, we find this region can be best described with a = (0.173 ± 0.007) ′′ = 0.94 ± 0.04 kpc, an axis ratio of 0.70 ± 0.04, and a PA of 54 ± 4 degrees. We note that this is an unconvolved Gaussian fit to the observed shape, we have not considered the PSF extent or shape to determine an intrinsic size.
From the integrated spectrum, this region has a relative velocity shift of +65 ± 51 km/s with respect to the quasar's peak [O iii] 5007 emission at = 6.  Tables 3  and 4. This region is spatially offset from the quasar such that narrow flux from the central spaxel, smeared by the PSF, does not contribute strongly. However, some quasar flux is present, hence the slightly different results for the different quasar subtractions.
Region 2 is the host galaxy of the quasar. The circle defining this region in Figure 5 has a radius = 0.458 ′′ or 2.5 kpc, for reference. By fitting a 2D Gaussian to the spatial distribution of the [O iii] 5007 flux in this region, we find this can be best described with a = 0.14 ± 0.02 ′′ = 0.75 ± 0.09 kpc, an axis ratio of 0.81 ± 0.14, and a PA of 83 ± 20 degrees. The GN H , [O iii] 4959, 5007, and H lines are detected with a SNR ≥ 17, while the brighter [O iii] 4959, 5007 lines also exhibit a broader (GB) wing, detected with SNR ≥ 9. The GN line is offset from the quasar at = 6.818 by +163 ± 51 km/s, with a line of 314 ± 51 km/s, while the broader wing is offset by −783 ± 51 km/s with a line of 526 ± 51 km/s. These line properties for Region 2 vary significantly when considering the BLR quasar subtracted cube, that has removed only the BLR model of the quasar, instead of the N&B quasar subtracted cube, in which unresolved NLR flux from the centre of the galaxy is also removed. The line properties for the BLR subtracted cube are presented in Table 3, alongside the N&B subtracted properties. The integrated spectrum of this region is also shown for the BLR-subtracted cube in Figure 5, and clearly contains significantly more flux than for the N&B quasar subtracted cube. This central region is coincident with the quasar, and so this BLR subtracted cube contains the unresolved narrow line flux from the NLR surrounding the quasar, spread across the image according to the PSF shape. This can also be seen in the flux maps for this BLR-subtracted cube, presented in Figure B.1, which show the bright central narrow line emission. In the BLRsubtracted cube there is significant GB flux, offset by −374 ± 51 km/s from the quasar's GN [O iii] 5007 peak, with a line of 807 ± 51 km/s. Region 3 sits to the northwest of the quasar. The ellipse we have fit by eye to define this region in Figure 5 is offset 0.90 ′′ = 4.9 kpc from the quasar centre-to-centre, with = 0.509 ′′ = 2.8 kpc, = 0.  narrow-line flux from the quasar has negligible contribution.

VDES J0020−3653:
Our quasar-subtracted cube for VDES J0020−3653 also reveals three extended narrow-line emission regions, with their [O iii] 5007 flux and integrated spectra shown in Figure 7. The [O iii] 5007 maps made using the higher spatial resolution 0.05"/pix cubes (Figure C.1) expose a fourth emission line region, separated by only ∼ 0.37 ′′ or 2 kpc (quasar peak to companion peak) to the south-southwest of the quasar. This is so close to the quasar host that in the 0.1"/pix cubes this appears to be one elongated host galaxy. Using the spatial structure from the higher resolution cube, we define four region ellipses from which we extract the spectra in each region. Thus, while the 0.1"/pix cubes can barely distinguish the host and this nearby companion galaxy, using the additional spatial information from the higher resolution cubes we are able to more accurately measure the host galaxy properties from the 0.1"/pix cubes, excluding the additional flux from the nearby companion. The extracted flux values from each of the four region ellipses are given in Table 3, with their flux ratios and derived physical properties in Table 4. Our kinematic maps for the subtracted cube, showing the flux of the lines, their relative velocity, and linewidth in each spaxel are shown in Figure 6 for the N&B quasar subtraction and Figure  B.1 for the BLR quasar subtraction.
Region 1 is a companion galaxy to the south-southwest of the quasar, with such a small separation that it is barely resolvable in the 0.1"/pix cube. The ellipse defining this region in Figure 7 is offset 0.52 ′′ from the quasar position. In WMAP9 cosmology at = 6.854, 1 ′′ = 5.429 kpc, and so this offset is 2.8 kpc, quasar-tocentre. The ellipse has = 0.338 ′′ = 1.8 kpc, = 0.228 ′′ = 1.2 kpc and a PA of 109 degrees. As this region is blended with the host emission at this small separation, the emission does not follow a 2D Gaussian profile and so we do not quote the best fitting parameters as for the other regions. From the integrated spectrum, this region has a relative velocity shift of −30 ± 51 km/s with respect to the quasar's peak [O iii] 5007 emission at = 6.854. The GN H , [O iii] 4959, 5007, and H lines are all detected with a SNR ≥ 11, with a line of 319 ± 51 km/s. These lines are well-fit by a single Gaussian component, with no secondary GB flux detected. When using the BLR subtracted cube, the flux is measured to be slightly larger due to contamination by the narrow-line flux from the quasar at this separation. The flux seen around [N ii] 6583 in Region 1 is likely remaining broad H flux from the quasar.
Region 2 is the host galaxy of the quasar. For reference, the circle defining this region in Figure 7 has a radius = 0.385 ′′ = 2.1 kpc. By fitting a 2D Gaussian to the spatial distribution of the [O iii] 5007 flux in this region, we find this can be best described with a = (0.18 ± 0.02) ′′ = (1.0 ± 0.1) kpc, an axis ratio of 0.66 ± 0.12, and a PA of 26 ± 10 degrees.  The integrated spectrum for the three regions (coloured lines), along with our best fit model in black. The middle panel, for Region 2, also shows the integrated spectra from the BLR quasar subtraction, in cyan, with the best fit model in grey. This region is the host of the quasar (position marked as a cyan cross in the top panel), and so the residual spectrum in this region depends heavily on the type of quasar subtraction. For Regions 1 and 3 the spectrum for both subtraction methods look equivalent, and so the BLR-subtracted spectra are not shown. For Region 2, we exclude the central 9 pixels surrounding the quasar peak, as these are highly corrupted by the quasar subtraction and introduce significant noise and artefacts. This means that we will slightly underestimate the total flux in this region, however the fluxes will be significantly more reliable than if these 9 most corrupted spaxels were included. In Region 1, the spectral lines can be described by a single narrow Gaussian (GN), whereas in Regions 2 and 3 there is also a broader component fit with a second Gaussian (GB); these separate components can be seen as the dashed black lines for the N&B quasar subtraction, and dotted grey lines for the BLR quasar subtraction.
As for the host region in DELS J0411−0907, the VDES J0020−3653 line properties for Region 2 vary significantly when considering the BLR quasar subtracted cube. The line properties for the BLR subtracted cube are presented in Table 3, alongside the N&B subtracted properties, and the integrated spectrum for Region 2 for the BLR subtraction is also shown in Figure 7. As for DELS J0411−0907, this central region is coincident with the quasar, and so this BLR subtracted cube contains the unresolved narrow line flux from the NLR surrounding the quasar, spread across the image according to the PSF shape. This can also be seen in the kinematic maps for this cube, presented in Figures 6  and B.1, which show the bright central emission.
Region 3 is likely a second companion galaxy to the northwest of the quasar. The ellipse defining this region in Figure 7 is offset 0.77 ′′ = 4.17 kpc from the quasar centre-to-centre, with = 0.417 ′′ = 2.3 kpc, = 0.271 ′′ = 1.5 kpc and a PA of 11 degrees. By fitting a 2D Gaussian to the spatial distribution of the [O iii] 5007 flux in this region, we find that this can be best described with a = (0.4 ± 0.1) ′′ = (2.0 ± 0.5) kpc, an axis ratio of 0.6 ± 0.2, and a PA of 193 ± 5 degrees. In this region there are significant detections of GN [O iii] 4959, 5007, H and H components, with SNR ≥ 8.9. The GN line component is offset from the quasar at = 6.854 by +232 ± 51 km/s, with a line of 250 ± 51 km/s. No broader GB component is detected. The properties of Region 3 are relatively unaffected when considering instead the BLR quasar subtracted cube; these are also given in Tables 3 and 4. This region is spatially offset from the quasar such that no significant narrow flux from the central region contributes.
Region 4 is a companion galaxy to the southwest of the quasar. The ellipse defining this region in Figure 7 is offset 1.03 ′′ = 5.61 Table 3. Fluxes and velocities for the H , H and [O iii] 4959, 5007 galaxy line components, from the integrated spectra over the regions shown in Figure 5 and 7 for DELS J0411−0907 (upper panels) and VDES J0020−3653 (lower panels). The values in braces give the integrated SNR of each line. We present the results from both the N&B quasar subtraction and the BLR quasar subtraction. We fit the narrow line spectra with two components, a true narrow line (GN) as well as a galaxy 'broad' wing (GB) component likely associated with outflows; note that this 'broad' component is not from the quasar BLR, which has been subtracted from these spectra. In regions where no broader GB line component is detected, we fit only a single GN component. Upper limits for undetected lines are 3 limits on the line flux, taken from the ERR array. Uncertainties on the velocity and velocity dispersion assume an uncertainty of ±1 wavelength element, ±6.65Å or 51 km/s. The velocity dispersions have been corrected for an instrumental broadening of FWHM inst = 115 km/s. kpc from the quasar position, quasar-to-centre. For reference, the region ellipse that we have fit by eye has = 0.439 ′′ = 2.4 kpc, = 0.299 ′′ = 1.6 kpc and a PA of 0 degrees. By fitting a 2D Gaussian to the spatial distribution of the [O iii] 5007 flux in this region, we find that this region can be best described with a = 0.43 ± 0.08 ′′ = 2.3 ± 0.4 kpc, an axis ratio of 0.47 ± 0.13, and a PA of −4 ± 2 degrees. We note again that this is an unconvolved Gaussian fit to the observed shape. From the integrated spectrum, this region has a relative velocity shift of +367 ± 51 km/s with respect to the quasar's peak [O iii] 5007 emission at = 6.854. The GN H , [O iii] 4959, 5007, and H lines are all detected with a SNR ≥ 13, with a line of 189 ± 51 km/s. These lines are well-fit by a single Gaussian component, with no GB flux detected. As for Region 3, the BLR subtracted cube has consistent results to the N&B quasar-subtracted cube, as Region 4 is spatially offset enough that the narrow-line flux from the quasar has negligible contribution.

Limits on Star Formation Rates and Excitation Mechanisms
Using our integrated spectrum for each region, we estimate the SFR using the GN H line flux. Assuming no extinction, a solar abundance, and a Salpeter IMF, the SFR can be estimated as: This equation assumes that all of the H emission is from star formation; for our host galaxies, however, we note that even for our N&B quasar subtracted cubes we expect some AGNphotoionized narrow line flux to contribute. The N&B quasar subtraction removes all of the GN component that is emitted from the central unresolved region surrounding the quasar, although spatially resolved GN flux from both the NLR and star-formation remains. Thus, while the N&B quasar subtracted cubes provide the most accurate measurement of the SFRs of our emission-line regions, in comparison to the BLR subtracted cubes, for our host components these are upper limits. For comparison, we also include the SFRs calculated using the BLR quasar subtracted cube in Table 4, which includes all quasar NLR emission, showing the significant improvement that the narrow-line quasar subtraction can have on being able to measure host galaxy properties. For DELS J0411−0907, the host SFR measured from the N&B subtraction is < 33 M ⊙ /yr, whereas the BLR subtraction can only constrain the SFR to < 91 M ⊙ /yr, a much larger upper limit. The host galaxy of VDES J0020−3653 has SFR constraints of < 55 M ⊙ /yr, consistent for both the BLR and N&B subtractions, likely due to the noisy subtraction around H . In the absence of dust extinction, theoretically the ratio of H to H flux is H / H ≃ 2.86, estimated using a temperature = 10 4 K and an electron density = 10 2 cm −3 for Case B recombination (Osterbrock 1989;Domínguez et al. 2013). Hence, the amount of dust can be estimated from these two lines via ( − ) = 1.97 log 10 [( H / H )/2.86] following Domínguez et al. (2013), and so H = (3.33 ± 0.80) × ( − ) assuming the Calzetti et al. (2000) reddening law. For companion regions where both H and H are significantly detected and H / H > 2.86, we calculate ( − ), H and a dustcorrected SFR, which are given in Table 4. As the host regions contain NLR quasar-photoionized flux which can have significantly higher H /H ratios than the 2.86 assumed here, we do not calculate dust-corrected SFR limits for the host galaxies.
For DELS J0411−0907, the companion galaxy Region 1 has a SFR of 27 ± 1 M ⊙ /yr, as measured from the N&B quasar subtracted spectrum. As the H / H ratio is < 2.86, we do not estimate the dust corrected values. These flux ratios are unlikely to be truly < 2.86, with this more likely being an issue with the quasar subtraction technique and its removal of flux, as the N&B subtraction has H / H = 2.4 ± 0.1 while the BLR subtraction has H / H = 2.8 ± 0.1. The BLR quasar subtraction results in slightly more H flux, corresponding to a SFR of 29 ± 1 M ⊙ /yr, as well as a lower H flux which also contributes to the larger H / H . For this region only, we therefore consider the BLR subtraction to provide the best estimate of the SFR, 29 ± 1 M ⊙ /yr, whereas for all other regions the N&B subtraction provides the best estimate.
For VDES J0020−3653, the companion galaxy Region 1 has a SFR of 24±2 M ⊙ /yr as measured from the N&B quasar subtracted spectrum, or 51 +50 −22 M ⊙ /yr when corrected for dust attenuation. Using the BLR subtraction, the SFR is similarly 26 ± 2 M ⊙ /yr. However, this subtraction results in a much lower H /H ratio of 2.91 ± 0.3, due to contamination from the unresolved quasar emission, and so in this case there is minimal change in the SFR when corrected for dust attenuation, 27 +15 −8 M ⊙ /yr . The companion galaxy Region 3 has a SFR of 18±2 M ⊙ /yr, as measured from the N&B quasar subtracted spectrum, or 70 +118 −40 M ⊙ /yr when corrected for dust attenuation, due to its large H /H ratio. Region 4 has a SFR of 17±1 M ⊙ /yr, as measured from both the N&B quasar subtracted and BLR quasar subtracted spectrum. When corrected for dust-attenuation, this SFR increases to 42 +42 −18 M ⊙ /yr. The subtraction method makes little difference for the companion regions 3 and 4, which are sufficiently spatially separated from the quasar such that the unresolved quasar emission has little contribution (see Table 4).
We note that our quasar hosts and companion galaxies have no detected [N ii] 6548, 6583 emission, implying low metallicities. For our SFR calculation we use the typically used conversion from H in Equation 7, which is calculated for solar abundances. showing the non-parametric central velocity of the line ( 50 ; middle) and the line width ( 80 ; right). As these are non-parametric, this combines both the GN and GB components. We describe our method for creating these maps in Section 2.3.4. The black ellipses in the lower middle panel depict the regions as in Figure 7. synthesis models. Assuming a case B recombination with nebular temperature of 10 4 K, a density of 100 cm −3 , and a Salpeter IMF, Lee et al. (2009) found that the ratio SFR/ H in Equation 7 decreases from 7.9 × 10 −42 at = ⊙ , to 6.2 × 10 −42 at = 0.2 ⊙ and 5.2 × 10 −42 at = 0.02 ⊙ . This would reduce our calculated SFRs by 22% and 34%, respectively.
We do not detect any [N ii] 6548, 6583 emission from our quasars, their hosts or their companions. To place upper limits on the amount of [N ii] 6583 present in the different emission-line regions, we integrate the ERR array at the expected location of the line, assuming a line width equal to that of [O iii] 5007. Our [N ii] 6583 limits are given in Table 3, with the resulting [N ii] 6583/H flux ratio limits in Table 4. Using our measured flux values and limits, we are able to place our emission line regions on a BPT diagram (Baldwin et al. 1981), which we show in Figure 8. We only consider regions where both the H and H lines are significantly detected (see Table 4). This diagnostic aims to characterise the main photoionization source based on the narrow line flux. For this we use the BLR quasar subtracted cube, as this gives the purest measure of the full narrow line flux, whereas in the N&B subtraction a large portion of the quasar NLR flux has been removed.
We find that for both quasars, the companion emission regions For VDES J0020−3653, the quasar host lies within the lowtransition region or star-forming regime. Theoretical models predict that high-AGN will reside in the left side of the BPT diagram, in the AGN regime that is sparsely populated at = 0 as well in the star-forming regime, mainly due to their lower metallicities (see e.g. Feltre et al. 2016;Strom et al. 2017;Nakajima & Maiolino 2022;Hirschmann et al. 2022). This is consistent with our measurements for both quasar host galaxies. For further discussion on this topic, see Übler et al. (2023).
To accurately determine the primary photoionization mechanism within these various regions, we would require either deeper observations for a [N ii] 6583 detection or tighter limit, or alternative line-ratio diagnostic methods. We could use a [ ii] 6717, 6731-based diagnostic (e.g. Veilleux & Osterbrock 1987), however [ ii] 6717, 6731 emission was not detected for these quasars. For the integrated quasar spectrum, by integrating the corresponding ERR array across the expected location of the [ ii] 6717 line, assuming the same line width as the narrow H component (Table 1)   The integrated spectrum for the four regions (coloured lines), along with our best fit model in black. The middle panel, for Region 2, also shows the integrated spectra from the BLR quasar subtraction, in cyan, with the best fit model in grey. This region is the host of the quasar (position marked as a cyan cross in the top panel), and so the residual spectrum in this region depends heavily on the type of quasar subtraction. We note that to account for the overlapping area between Regions 1 and 2, we assign the flux in those spaxels to Region 1 only. For Region 2, we exclude the central 9 pixels surrounding the quasar peak, as these are highly corrupted by the quasar subtraction and introduce significant noise and artefacts. This means that we will slightly underestimate the total flux in this region, however the fluxes will be significantly more reliable than if these 9 most corrupted spaxels were included. In Regions 1, 3, and 4, the spectral lines can be described by a single narrow Gaussian (GN), whereas in Region 2 there is also a broader component which we fit with a second Gaussian (GB); these separate components can be seen as the dashed black lines for the N&B quasar subtraction, and dotted grey lines for the BLR quasar subtraction. We note that the emission around 3.8 m in Regions 3 and 4 is likely an artefact in the cube. We do not see any significant detection of [N ii] 6548, 6583. correction, we find [ II] 6717 < 2.8 × 10 −18 erg/s/cm 2 (3 ) for DELS J0411−0907. For VDES J0020−3653, this line is beyond the wavelength coverage of the detector. These lines are either very close to or just above the detectable wavelength range for this G395H/F290LP NIRSpec configuration, making accurate measurements difficult if a line was detected. We have also not detected any He ii 4686 emission, which is particularly useful for distinguishing star formation and AGN dominated regions in the high-Universe where metallicities are low (see e.g. Nakajima & Maiolino 2022;Übler et al. 2023). For the integrated quasar spectrum, following the same process for [ ii] 6717 but instead assuming the same line width as the narrow [O iii] 5007 component (Table 1), we find HeII 4686 < 2.3 × 10 −18 erg/s/cm 2 (3 ) for DELS J0411−0907 and HeII 4686 < 1.6 × 10 −18 erg/s/cm 2 (3 ) for VDES J0020−3653.

Dynamical Masses and Black Hole-Dynamical Mass Ratios
To calculate the dynamical masses of the various emission regions, we follow the method as in Decarli et al. (2018) for ALMA quasar host studies. Firstly, we assume each region is dispersiondominated, with: where is the radius of the region, the integrated line width, and the gravitational constant. As an alternative, we assume that the regions are rotation-dominated thin discs, with: where 0.75 FWHM scales the line FWHM to the width of the line at 20% of the peak as in Willott et al. (2015), and is the inclination angle of the thin disc. Finally, we also consider the virial approximation dyn,virial ≈ 5 2 .
This commonly-used scaling factor of 5 was found to be the optimal scaling factor when calibrated to a sample of 25 galaxies by Cappellari et al. (2006) (see e.g. van der Wel et al. 2022, for further discussion). This estimate will be a factor of 10/3 larger than dyn,dispersion by construction. As it is most common for high-quasar host studies to assume a rotating disc geometry (e.g. Wang et al. 2013;Willott et al. 2015;Decarli et al. 2018), we choose our best dynamical mass estimate to be dyn,rotation .
To measure the radii, we fit a 2D Gaussian to the spatial distribution of the [O iii] 5007 emission from each region, and take to be the of the semi-major axis; these values are described in Section 3.2.1. From this Gaussian fit we also measure the axis-ratio, / , from which we derive the inclination angle, , via cos( ) = / . As above, we note that and / have not been deconvolved by the PSF. For the emission line widths we take the GN [O iii] 5007 component, which is the brightest of the lines we have measured. We note that this line width represents the total dynamics of the gas, including rotation as well as non-rotational components such as dispersion and winds (see e.g. Wisnioski et al. 2018). We consider only the N&B quasar subtracted cube for this calculation, since otherwise the measurements could be skewed by the unresolved NLR regions around the quasar. The kinematics of this central NLR are not dominated by the larger-scale galaxy potential, thus this emission does not trace the dynamical mass. Considering only the N&B quasar subtracted cube ensures that we are measuring the spatially extended host regions and not the PSF shape, which is particularly important for the host regions (Region 2). The values used in these calculations and the resulting dynamical mass estimates are listed in Table 5.
For DELS J0411−0907, for Region 1 we find a dynamical mass between 0.6-2.5×10 10 M ⊙ , with a best estimate of dyn,rotation = (2.5 +3.1 −1.6 ) × 10 10 M ⊙ , from an inclination angle = 45 ± 3 degrees. For the host Region 2, we find that the rotational measure for the dynamical mass, dyn,rotation = (16 +59 −9 ) × 10 10 M ⊙ , is larger than the two other estimates, dyn,dispersion = (2.6 +1.3 −1.0 ) × 10 10 M ⊙ and dyn,virial = (8.6 +4.4 −3.3 ) × 10 10 M ⊙ , albeit with large uncertainties particularly due to the inclination angle, = 36 +12 −17 degrees. This difference is due to our measured inclination angle for this region being relatively low, and so the inclination correction is relatively large. For their ALMA quasar host galaxies that are marginally resolved, and thus where their axis ratio cannot be accurately measured, Willott et al. (2015) and Decarli et al. (2018) assume an inclination of 55 degrees, the median inclination of the resolved sources in Willott et al. (2015) and Wang et al. (2013). As our quasar hosts are difficult to measure, requiring the quasar subtraction, and they are interacting systems which complicates their structure, we are not confident in our inclination measurements. If we instead assume the median = 55 degrees, the rotational dynamical mass is (8.0 +4.1 −3.0 ) × 10 10 M ⊙ , similar to the virial estimate. For Region 2 of DELS J0411−0907, we therefore assume a best dynamical mass estimate of dyn,virial = (8.6 +4.4 −3.3 ) × 10 10 M ⊙ . Region 3 cannot be well described by a 2D Gaussian, and so we cannot perform our fit as for the other regions. VDES J0020−3653 Region 1 also cannot be well described by a 2D Gaussian, due to contamination by the quasar host, and so we do not estimate its dynamical mass. For the VDES J0020−3653 host Region 2, we find a dynamical mass between 4.6-17 ×10 10 M ⊙ , with a best estimate of dyn,rotation = (17 +18 −8 ) × 10 10 M ⊙ , from an inclination of = 49 +8 −10 degrees. This is similar to the median inclination of the resolved sources in Willott et al. (2015) and Wang et al. (2013), = 55 degrees, which would reduce the value slightly to dyn,rotation = (14 +7 −5 ) × 10 10 M ⊙ . Region 3 has dynamical mass estimates of 4.7-15×10 10 M ⊙ , with our best estimate dyn,rotation = (15 +39 −10 ) × 10 10 M ⊙ from an inclination of 54 +15 −19 degrees. For VDES J0020−3653 Region 4 we find a dynamical mass between 3.1-10×10 10 M ⊙ , with a best estimate from dyn,rotation of (8.1 +11.2 −5.1 ) × 10 10 M ⊙ with an inclination of = 62 ± 9 degrees. Table 5. Dynamical masses, and values used in their calculation, for the emission-line regions in DELS J0411−0907 (upper rows) and VDES J0020−3653 (lower rows), using the N&B quasar subtraction. DELS J0411−0907 Region 3 and VDES J0020−3653 Region 1 cannot be well described by a 2D Gaussian, and so we do not measure their dynamical mass. We plot the dynamical mass-black hole mass relation in Figure 9 for these two quasars. We find that both of our quasars lie above the local Kormendy & Ho (2013) black hole-stellar mass relation, although these are consistent within the scatter and uncertainties. Previous ALMA observations have found that luminous quasars generally lie above the local black hole-host relation, while lower luminosity quasars lie around or below the relation (see e.g. Willott et al. 2017;Izumi et al. 2019). Our luminous quasars lie above the relation, which is consistent with this existing picture, however a larger sample of high-quasars is required to test this scenario with NIRSpec. For a true comparison to the local black hole-stellar mass relation, we require stellar and not dynamical masses.
We note that in the current observations the faint host galaxy continuum emission is not detectable, and so we cannot measure their stellar masses. Our high spectral resolution ( ≃ 2700) observations were designed to study the emission lines in detail. This sacrifices the sensitivity that is required to detect the host continuum, which is expected to be achievable in similar exposure times at lower spectral resolution (either ∼ 100 IFU spectra, or alternatively images as in e.g. Ding et al. 2022;Harikane et al. 2023). Future work within the GA-NIFS GTO program will use the ∼ 100 prism mode to search for the continuum emission of several high-quasar host galaxies, which would make stellar mass measurements possible.
In Figure 9 we also find that both of our quasars have dynamical masses and black hole-dynamical mass ratios consistent with ALMA observations of ≳ 6 quasars. However, we note that DELS J0411−0907 and VDES J0020−3653 do not have published ALMA dynamical mass estimates that could be used for a direct comparison of these methods for our targets. A detailed comparison of ALMA and JWST NIRSpec dynamical mass and JWST stellar mass estimates will be a key aim of future quasar studies, which will allow for a much deeper understanding of the high-black hole-host mass relations.

Outflow properties
The clear detection of an outflow component, as traced by the broad kinematic component (GB components in Figure 2) via the analysis of the [O iii] 5007 velocity and spatial distribution, opens the possibility of investigating the ionized outflow properties in these early quasars. We note that this is the first time that it is possible to properly characterise the ionized outflows in quasars in the Epoch of Reionisation. Previous studies have relied on spatially unresolved information, and in particular on the detection of UV broad and blueshifted absorption features of UV resonant lines (e.g. Wang et al. 2018;Ginolfi et al. 2018;Onoue et al. 9.5 10 2019; Bischetti et al. 2022Bischetti et al. , 2023. These features trace primarily the nuclear outflow, which is generally a small fraction of the global, galactic scale outflow, and the lack of spatial information generally prevents a reliable estimate of the outflow parameters such as the outflow rate, kinetic power and momentum rate. We also note that NIRSpec fixed-slit spectroscopic observations of these two quasars covering the rest-frame UV do not reveal any broad absorption features (Willott, priv.comm.), indicating that the presence of ionized outflow is not necessarily traced by broad absorption features in the UV part of the spectrum (see e.g. Rankine et al. 2020).
A&A proofs: manuscript no. manuscript Table 6. The properties of the quasar outflows as measured in Section 3.3, for DELS J0411−0907 and VDES J0020−3653. The luminosities and velocities are taken from the H and [O iii] 5007 GB components from the full integrated quasar spectrum ( Figure 2 and  Figure 2). In this section we use the flux, velocity and spatial information of the H and [O iii] 5007 broad galaxy (GB) components to infer some of the basic properties of the ionized outflow. We use the GB fluxes and velocities measured from our full integrated spectral fit for the quasar and host galaxy, prior to quasar subtraction, as described in Section 2.3.1 and shown in Figure 2. This gives the most reliable measurements of the integrated quantities, as these are not affected by the quasar subtraction technique. To measure the spatial extent of the outflow, we use the quasar BLR-subtracted cube, from which we can measure the [O iii] 5007 line width in each spaxel without contamination from the quasar's BLR. For this we consider the full cube, without splitting the emission into regions as in Section 3.2 above. The broad flux emanates from the quasar, approximately covering Region 2 for both quasars, as well as extending into Region 3 for DELS J0411−0907-this can be seen in the 80 maps in Figure from Carniani et al. (2015), assuming a density of the gas in the ionized outflow of ⟨ ⟩ = 500cm −3 , as in Carniani et al. (2015), and where is a factor taking into account the clumpiness of the medium ( = ⟨ ⟩ 2 /⟨ 2 ⟩), which we assume to be unity. Here we consider the H to be the luminosity of the GB H component as measured from our full integrated spectral fit for the quasar and host galaxy, as described in Section 2.3.1 and shown in Figure 2. We note that as in our calculations above, we do not correct for extinction. We perform the equivalent calculation using the and assuming a solar oxygen abundance [O/H] = [O/H] ⊙ . These luminosities and calculated outflow masses are given in Table 6. By using the luminosities of the "outflow" component of H we infer outflow masses of ionized gas of (17 +13 −11 ) × 10 7 M ⊙ and (97 +14 −17 )×10 7 M ⊙ for DELS J0411−0907 and VDES J0020−3653, respectively. We note that the masses using the broad, blueshifted component of [O iii] 5007 are lower than the H masses, by a factor of approximately 2 for DELS J0411−0907 and an order of magnitude for VDES J0020−3653. This is a well known issue and it is understood to be a consequence of the ionization structure of the clouds in the NLR of AGNs, whereby the [O iii] 5007 only traces a small part of the mass of each cloud (e.g. Perna et al. 2015).
Once the outflow mass is inferred, the ionized outflow rate can be calculated under the thin shell approximation (which generally provides a better description of galactic outflow rates, relative to the uniformly filled outflows, as discussed in Lutz et al. 2020) by using the equation where out is size of the outflow and out is its velocity. To estimate out , we measure the maximum extent from the quasar centre of the [O iii] 5007 emission with 80 > 600 km/s, in the BLRsubtracted cube. This 80 > 600 km/s is a conservative criterion that has been used in previous studies to classify outflows (e.g. Harrison et al. 2015;Kakkad et al. 2020). We note that for DELS J0411−0907 this outflow region extends from the host galaxy in Region 2 into Region 3, which is seen to have an outflow GB component in Figure 5. For VDES J0020−3653, when we measure the maximum extent we exclude the south-southwest direction, avoiding the nearby Region 1 which may be artificially contributing to the large 80 values as its velocity is offset from the quasar. This out is a projected radius and it reproduces the true radius only if the directions of the outflow expansion include the plane of the sky. For the outflow velocity we adopt the same approach as Fluetsch et al. (2019) by defining out = peak + FWHM /2, where peak is the peak velocity of the GB component of [O iii] 5007 (and equivalently H as they are kinematically tied in our spectral fits) and FWHM its FWHM. These velocities are measured from our full integrated spectral fit for the quasar and host galaxy, which are given in Table 1. Our measured out and out are given in Table 6. Based on these measured values we obtain ionized outflow rates from H of 58 +44 −37 M ⊙ yr −1 and 525 +75 −92 M ⊙ yr −1 for DELS J0411−0907 and VDES J0020−3653, respectively; the [O iii] 5007 rates are also quoted in Table 6.
The implications of these findings are discussed in Section 4.2.

Quasar companions
Both of our quasars show additional emission-line regions nearby their host galaxies, which are consistent with being companion galaxies.
For DELS J0411−0907, Region 1 to the southeast of the quasar has a projected offset of ≃ 4.3 kpc from the quasar, with a widest extent of ≃ 8.8 kpc. From the integrated spectrum we find that this region has a velocity offset relative to the quasar host of −98 ± 72 km/s, indicating minimal line-of-sight offset ( Figure 5 and Table 3). From the emission maps in Figures 4 and  B.1 we see that the velocity spans about 100 km/s either side of the host velocity. This companion galaxy is likely undergoing a merger with the quasar host. Region 3 exhibits some broader [O iii] 4959, 5007 emission, indicating an ionized gas outflow in this region, and is likely a tidal tail.
For VDES J0020−3653, Region 1 to the south-southwest of the quasar has a projected offset of ≃ 2.8 kpc from the quasar, with a widest extent of ≃ 3.6 kpc. This region has a line-of-sight velocity offset by −55±72 km/s from the host (Figures 6 and B.1). Region 3 to the northwest of VDES J0020−3653 has a projected offset of ≃ 4.2 kpc from the quasar, with a widest extent of ≃ 4.6 kpc. This region has a line-of-sight velocity offset by +83 ± 72 km/s from the host, with the bulk of the gas at a similar offset. Region 4 to the southwest of the quasar has a projected offset of ≃ 5.6 kpc, and a widest extent of ≃ 4.8 kpc. This region has a line-of-sight velocty offset by +218 ± 72 km/s from the host, again with the bulk of this gas at a similar velocity.
All of these regions satisfy the common merger condition of having a projected distance of < 20/ℎ kpc, where ℎ is the dimensionless Hubble constant ℎ = 0 /(100 km/s/Mpc), and a velocity difference < 500 km/s (e.g. Patton et al. 2000;Conselice et al. 2009).
The presence of quasar companions, or lack thereof, provides important insights into quasar growth mechanisms. For example, theory predicts that kpc-scale interactions trigger AGN activity (e.g. Sanders et al. 1988;Hopkins et al. 2006), however mergers and companion galaxies are not ubiquitously observed around low-quasars (Cisternas et al. 2011;Kocevski et al. 2012;Mechtley et al. 2016;Marian et al. 2019). However, these early quasars have more intense accretion, and higher merger rates are likely in the early Universe, suggesting that these mergers could be a key driver of this early black hole growth. ALMA observations have detected companion galaxies around high-quasars at ≃ 8-60 kpc separations which are interpreted as major galaxy interactions (e.g. Wagg et al. 2012;Decarli et al. 2017), at rates of up to 50% in some samples (Trakhtenbrot et al. 2017;Nguyen et al. 2020). These companions are less frequently observed in the rest-frame UV (Willott et al. 2005), although some potential companions have been discovered with HST (McGreer et al. 2014;Marshall et al. 2020). The lack of UV-discovered companions may be due to dust obscuration (Trakhtenbrot et al. 2017), however our observations which reveal such close neighbours may indicate that resolution has previously limited our ability to identify merging quasars in the early Universe. It is indeed striking that both of our quasars have these companion galaxies, and it will be insightful to see whether other upcoming JWST observations also reveal similar structures around other highquasars, building up a larger statistical sample. Discovering these close companions highlights the ability of JWST to understand these early quasars in significantly greater detail, through detailed measurements of their environments and thus potential quasar triggering mechanisms (see also Perna et al. 2023).

Quasar driven outflows
In Figure 10 we show the ionized outflow rate as a function of AGN luminosity for our two quasars at ∼ 6.8, and for AGN whose ionized outflows had been measured at lower redshift (from the compilation of Fiore et al. 2017). Clearly, our two high-quasars follow the same relation as lower redshift quasars.
The inferred outflow rates (∼ 58M ⊙ yr −1 for DELS J0411−0907 and ∼ 525M ⊙ yr −1 for VDES J0020−3653) are significantly larger than the upper limits we have measured for the host SFRs (< 33M ⊙ yr −1 for DELS J0411−0907 and < 54M ⊙ yr −1 for VDES J0020−3653; see Section 3.2.2). These outflow rates are much larger than the SFRs, especially if we consider that the measured outflow rates are only for the ionized component and that the other phases, in particular the molecular phase, may also contribute. This result implies that, if prolonged, the quasar-driven outflow may be the main agent responsible for depleting the galaxy of its gas content. However, one has to take into account that the quasar activity is likely short lived and bursty, hence it may not continue for a time long enough to compete with star formation in depleting the gas mass.
More quantitative consideration on the potential feedback effect of the quasar driven outflow can be obtained by considering the kinetic power of the outflows. Efficient feedback, with high coupling with the ISM, in many scenarios requires that the kinetic power of the outflow should be about 5-7% of the AGN bolometric luminosity. The right panel of Figure 10 shows the kinetic power of ionized outflows as a function of AGN luminosity, where locii of constant P /L AGN ratio are also indicated with dashed lines. Clearly, the ionized outflows of two quasars at z∼6.8 have a kinetic power that is less than 1% of L AGN (in the case of DELS J0411−0907 actually <0.1%), indicating that the outflow may not be an effective feedback mechanism.
Moreover, the momentum rate ( out out ) of the two quasars is 0.07 L AGN / and 1.28 L AGN / for DELS J0411−0907 and VDES J0020−3653, respectively. These values are much lower than expected by the (very feedback-effective) energy-driven outflows, which expect a momentum boost by a factor of ∼20 relative to the radiation momentum rate L AGN /c. Based on these values, these outflows seem to be more consistent with momentum-driven outflows, possibly resulting from radiation pressure on the dusty clouds (Fabian 2012;Costa et al. 2014). However, once again, we caution that with our NIRSpec data we can only probe the ionized phase of these outflows. The other phases of the outflows, and in particular the molecular component, may be much more massive and energetic (Fluetsch et al. 2019), although Fiore et al. (2017) suggest that at such high luminosities the ionized outflows may be comparable or stronger than the other phases. Future deep observations with ALMA may provide constraints on the other phases of these outflows and, therefore, provide additional information on the feedback role of AGN-driven outflows in these two early quasars.
By analysing their integrated quasar spectrum, we measured black hole masses from the H and H broad lines, which are the most reliable single-epoch mass estimators. These are some of the first observations of Balmer broad emission lines in high-quasars, making such a measurement possible for the first time. We derive a best estimate for the black hole mass in DELS J0411−0907 of BH = 1.85 +2 −0.8 × 10 9 M ⊙ , implying an Eddington ratio of Edd = 0.8 +0.7 −0.4 . This is more massive than the Mg ii mass estimates of BH,MgII = (0.6-1) × 10 9 M ⊙ which implied more extreme Eddington ratios of Edd = 1.3-2.4. By measuring a more accurate black hole mass, we have revised our  Fig. 10. Left: Ionized gas outflow mass rate as a function of AGN luminosity, for the two quasars at ∼ 6.8 presented in this paper, and for various other AGN at lower redshift (0 < < 3; from Fiore et al. 2017). The solid line is a linear fit (in logarithmic scales) of the latter data. The ionized outflows of the two quasars at ∼ 7 follow the same relation as quasars at lower redshifts. Right: Kinetic power (P ) of the ionized outflows as a function of AGN luminosity (same symbols as in the left plot). The dashed lines show constant ratios of P and L AGN as labelled. We note that the Fiore et al. (2017) galaxies have been shifted by log(5/2) for consistency with our density assumption of 500cm −3 , as they assume 200cm −3 .
understanding of the accretion physics of this early quasar. For VDES J0020−3653 our best estimate of the black hole mass is BH = 2.9 +3.5 −1.3 × 10 9 M ⊙ , corresponding to an Eddington ratio of Edd = 0.4 +0.3 −0.2 . This is more consistent with the Mg ii measurements of Reed et al. (2019). Together, these quasars show that the Mg ii black hole mass can be significantly biased relative to the H estimate for some high-quasars, while for others these two techniques may give consistent black hole mass measurements, highlighting the benefit of measuring these more reliable H masses.
We perform quasar subtraction on our IFU observations, by modelling the PSF using the variation in the flux of the quasar broad H and H lines across the spatial domain. From this subtraction, we reveal that the host galaxy of DELS J0411−0907 shows emission extending to a radius of 2.5 kpc, described by a 2D Gaussian with semi-major axis = 0.75 kpc. This host has narrow [O iii] 4959, 5007, H and H emission, alongside a broader wing emission from an outflow. We estimate that the host has a SFR of < 33 M ⊙ /yr, with a dynamical mass in the range of 2.6-16×10 10 M ⊙ with a best estimate of 8.6 +4.4 −3.3 × 10 10 M ⊙ . This host galaxy lies within the AGN-dominated regime of the low-BPT diagram.
The host galaxy of VDES J0020−3653 has emission extending to a radius of 2.1 kpc, described by a 2D Gaussian with semi-major axis = 1.0 kpc. Both narrow [O iii] 4959, 5007, H and H and broad-wing emission is detected. We estimate a SFR of < 54 M ⊙ /yr for this host, with a dynamical mass between 4.6-17×10 10 M ⊙ with a best estimate of 17 +18 −8 × 10 10 M ⊙ . The host of VDES J0020−3653 lies within the star-forming or transition region of the low-BPT diagram, however this is likely due to the low metallicity of these quasar hosts, with low-metallicity, high-AGN predicted to lie in this region.
We find that from our measured black hole and dynamical masses, DELS J0411−0907 and VDES J0020−3653 both lie slightly above the local black hole-host mass relation (Kormendy & Ho 2013) and are consistent with the existing observations of ≳ 6 quasar host galaxies with ALMA.
Both quasars are found to have additional emission-line regions near their host galaxies, with two regions for DELS J0411−0907 and three for VDES J0020−3653. We study these companion regions, measuring their SFRs, excitation mechanisms, and kinematics. We find that four companion regions are likely to be either galaxies merging with the quasar host galaxy, and one an extended tidal tail. It is striking that both of our highquasars are undergoing galaxy mergers, and it will be insightful to see the fraction of high-quasars that are found to also have companions in future NIRSpec IFU observations.
We detect ionized outflows in [O iii] 4959, 5007 and H from both VDES J0020−3653 and DELS J0411−0907. We measure the rates of the ionized outflow from H to be 58 +44 −37 M ⊙ yr −1 for DELS J0411−0907 and 525 +75 −92 M ⊙ yr −1 for VDES J0020−3653. These outflow rates are much larger than the SFRs of the host galaxies, implying that if this is prolonged, the quasar-driven outflow may be the main driver for gas depletion within the galaxies. We find that these outflows are likely momentum-driven, and are not likely to be effective feedback mechanisms.
This work highlights the exquisite capabilities of the JWST NIRSpec IFU for observing high-quasars. With the highresolution spectral coverage in the infrared, we are able to study key rest-frame optical emission lines that have been previously unobservable, such as H , [O iii] 4959, 5007 and H . This allows us to measure the black hole masses and Eddington ratios of these quasars more accurately than has been previously possible. Obtaining spatially resolved IFU spectra also allows us to study the host galaxies, local environments, and outflow properties of these quasars in great detail. Overall, this work shows the possibilities of the NIRSpec IFU, and we look forward to the insights that this instrument will give about the formation and growth of quasars and their host galaxies in the early Universe. , from the N&B quasar subtraction, as described in Section 2.3.3. These are created from the best fit quasar cube, which is the peak quasar spectrum scaled in flux to match the observed broad line flux in each spaxel. The relative intensity here depicts this scaling factor, which is the amount of flux in each spaxel relative to the brightest spaxel. We create two PSFs for each quasar, one from each of H and H , as the PSF varies with wavelength. and the line width ( 80 ; right). As these are non-parametric, this combines both the GN and GB components. We describe our method for creating these maps in Section 2.3.4. The black ellipses in the lower middle panels depict the regions as in Figures 5 and 7.