GA-NIFS: Early-stage feedback in a heavily obscured AGN at z = 4 . 76

Dust-obscured galaxies are thought to represent an early evolutionary phase of massive galaxies in which the active galactic nucleus (AGN) is still deeply buried in significant amounts of dusty material and its emission is strongly suppressed. The unprecedented sensitivity of the James Webb Space Telescope (JWST) enabled us for the first time to detect the rest-frame optical emission of heavily obscured AGN and unveil the properties of the hidden accreting super-massive black holes (BHs). In this work, we present the JWST / NIRSpec integral field spectroscopy (IFS) data of ALESS073.1, a massive (log(M ⋆ / M ⊙ ) = 10 . 98) dusty, star-forming galaxy at z = 4 . 755 hosting an AGN at its center. The detection of a very broad ( > 9000 km s − 1 ) H α emission associated with the broad line region (BLR) confirms the presence of a BH (log( M BH / M ⊙ ) > 8 . 7) accreting at less than 18% of its Eddington limit. The identification of the BLR classifies the target as a type 1 AGN despite the observed high column density of N H ∼ 10 24 cm − 2 . The rest-frame optical emission lines also reveal a fast ( ∼ 1700 km s − 1 ) ionized gas outflow marginally resolved in the galaxy center. The high sensitivity of NIRSpec allowed us to perform the kinematic analysis of the narrow H α component, which indicates that the warm ionized gas velocity field is consistent with disk rotation. Interestingly, we find that in the innermost nuclear regions ( < 1 . 5 kpc), the intrinsic velocity dispersion of the disk reaches ∼ 150 km s − 1 , which is ∼ 2 − 3 times higher than the velocity dispersion inferred from the [C ii ]158 µ m line tracing mostly cold gas. Since at large radii the velocity dispersion of the warm and cold gas are comparable, we conclude that the outflows are injecting turbulence in the warm ionized gas in the central region, but they are not su ffi ciently powerful to disrupt the dense gas and quench star formation. These findings support the scenario that dust-obscured galaxies represent the evolutionary stage preceding the unobscured quasar when all gas and dust are removed from the host.


Introduction
Supermassive black holes (BHs) are thought to reside in the center of the majority of local massive galaxies (Hopkins et al. 2008).A symbiotic connection between the growth of BHs and their hosts is suggested by the observed relations between the BH masses and the galaxy properties (e.g., Kormendy & Ho 2013;Heckman & Best 2014).Specifically, observations of local galaxies reveal tight relations between the BH mass and the stellar velocity dispersion and the mass and luminosity of the galactic bulge.These relations hold up throughout several orders of magnitude in black hole masses and galaxy properties and up to high redshift with different normalizations (Carraro et al. 2020;Suh et al. 2020).
During the accretion phase, BHs are revealed as active galactic nuclei (AGNs) due to the radiation emitted from radio to Xray wavelengths by the accretion disk (Padovani et al. 2017).X-ray surveys show that the most luminous and massive AGNs were most numerous at z > 1 and, in particular, observations find the cosmic black hole accretion rate density peaks at z ∼ 2, such as the cosmic star formation rate (SFR) density (Shankar et al. 2009;Delvecchio et al. 2014;Madau & Dickinson 2014;Aird et al. 2015;Ananna et al. 2019;Brandt & Yang 2022).Both the cosmic star formation rate and BH accretion rate are driven by the availability of cold gas in the system (Hopkins et al. 2008).
⋆ E-mail: eleonora.parlanti@sns.it The tight relations between the BH and host galaxy properties and the similar evolution of activity with redshift suggests that BHs and the galaxies they inhabit have undergone a common evolutionary process (Kormendy & Ho 2013).
The origin and the mechanisms that regulate the coevolution of BHs with their host galaxy are still unclear.Massive galactic outflows driven by the radiation emitted by the most luminous AGNs are considered a fundamental physical process in the evolution of galaxies.They are believed to regulate star formation (SF, e.g., Fabian 2012;Zubovas & King 2014;Muratov et al. 2015;Yuan et al. 2018;Yoon et al. 2018;Nelson et al. 2019) and reduce the number of galaxies at the high-mass end of the stellar mass function (Benson et al. 2003;Puchwein & Springel 2013).Such fast outflows can potentially accelerate a substantial mass of gas beyond the escape velocity of the local gravitational potential, inject turbulence in the interstellar medium (ISM), and/or heat the gas in the galaxy, in this way damping or even halting SF in their host galaxies ("negative feedback" ; e.g., Fabian 2012;Harrison et al. 2017).On the other hand, AGN feedback has been observed to enhance star formation ("positive feedback"; e.g., Shin et al. 2019), with stars actively forming in outflowing material (e.g., Maiolino et al. 2017;Gallagher et al. 2019).It is thus fundamental to investigate the feedback mechanism over the various phases of the evolution of the BHs and galaxies to understand their role in the coevolution process.
Submillimeter galaxies (SMGs) are a class of high-redshift (z > 0.1) galaxies mostly characterized by a high luminosity in the far infrared (FIR) continuum emission (L FIR > 10 11 L ⊙ ), high SFRs (SFR ∼ 10 3 M ⊙ yr −1 ), and high dust content (M dust = 10 8−10 ; Santini et al. 2010).Although SMGs are a rare cosmological population of galaxies, they account for ∼ 20% of the SFR cosmic density at 1 < z < 5 (Swinbank et al. 2014).These galaxies are also thought to be the precursors of local massive and quiescent early-type galaxies (ETGs) that host BHs with masses M BH > 10 8 M ⊙ in their centers (Sanders et al. 1988;Swinbank et al. 2004;Hopkins et al. 2008;Toft et al. 2014).Quiescent galaxies are already common at (z ∼ 2-3), and are observed up to z ∼ 4.6 (Carnall et al. 2023).The evolutionary path that connects SMGs and quiescent galaxies is thought to be driven by the interaction between the galaxy and the BH (Sanders et al. 1988;Hopkins et al. 2008).The current theoretical scenarios expect that the feeding of galaxies with gas from the cosmic web or via mergers triggers both episodes of intense SF and BH accretion.This process is thought to be selfregulated.In particular, the feedback from the accreting BHs should be able to balance the SF in the host galaxy and the accretion of gas on the BH itself due to powerful radiation-driven outflows that are able to sweep out the host galaxy gas reservoir (Debuhr et al. 2012) and halt the accretion of gas (Peng et al. 2015).These outflows are thought to become important when the accretion rate on the BH reaches the Eddington rate (King & Pounds 2015).The dust in the system is also swept out, allowing the radiation coming from the BH accretion zone to be detected and identified as a luminous quasar (QSO) (Sanders et al. 1988;Hopkins et al. 2008).The removal and heating of gas due to the AGN feedback halts SF in the system, turning the host into a "red and dead" galaxy.This is consistent with what is observed in local ETGs, which are dominated by old stellar populations, with formation redshifts of z > 2 (Thomas et al. 2005;McDermid et al. 2015).Several studies also find galaxies whose gas content has been depleted on short timescales (Bezanson et al. 2019;Williams et al. 2021).Moreover, recent observations of quiescent galaxies at z > 3 indicate that fast quenching processes are already in place in the first 2 Gyr of the Universe (Glazebrook et al. 2017;Valentino et al. 2020;Carnall et al. 2023).In conclusion, it is crucial to study the SMG population to asses whether they really represent the evolutionary stage preceding the active QSO phase and, thus, whether they are the progenitors of local ETGs or not.Moreover, SMGs enabled us to study the first phases of the BH feedback process outlined above.
The ISM properties of the high-redshift SMG population have been mainly studied through millimeter observations, which do not allow us to determine whether a BH is hidden at the center of the galaxies.Up until the advent of the James Webb Space Telescope (JWST), the rest-frame optical emission lines from galaxies at z > 4 were extremely difficult to access from ground-based telescopes because the lines are redshifted to wavelength which are outside the atmospheric windows.Thanks to its unique sensitivity, the NIRSpec instrument (Jakobsen et al. 2022) on board JWST has already proven its capabilities to detect faint z > 4 galaxies and identify serendipitous AGNs by observing the emission that arises from the broad line regions (BLRs) surrounding super-massive BHs (Kocevski et al. 2023;Maiolino et al. 2023;Matthee et al. 2023;Übler et al. 2023) By using the Integral Field Spectrograph (IFS) mode of NIR-Spec (Böker et al. 2022) we can also exploit the rest-frame optical lines at high redshift to spatially resolve the emission from the ISM and stellar population and determine the impact of the feedback mechanism on the host galaxy (Cresci et al. 2023;Marshall et al. 2023;Perna et al. 2023;Übler et al. 2023).
This work is structured as follows.In section 2, we describe the target and the new JWST IFS observations.In section 3, we present the analysis of the spatially integrated and singlespaxel spectra.In section 4, we investigate which is the primary excitation mechanism of the gas.In Section 5 we compute the black hole mass and compare the position of ALESS073.1onthe M ⋆ − M BH plane with other low and high redshift AGNs and QSOs.In section 6, we study the properties and energetics of the outflow.In section 7, we perform a detailed kinematic analysis of the host galaxy.We discuss our results in section 8, and we draw our conclusions in section 9.In this work, we adopt the cosmological parameters from Planck Collaboration et al. ( 2016): H 0 = 67.7 km s −1 Mpc −1 , Ω m = 0.307, and Ω Λ = 0.691, 0.1 ′′ = 0.66 kpc at z =4.755.

Target
ALESS073.1 is part of the Extended Chandra Deep Field South (Lehmer et al. 2005) (RA: 03h32m29.290s,DEC:-27d56m19.30s).In Figure 1 we show the RGB image of the source.It was identified as a strong submillimeter (Coppin et al. 2009;De Breuck et al. 2011) and X-ray source (Gilli et al. 2011(Gilli et al. , 2014) ) as well as Lyα emitter (Vanzella et al. 2006(Vanzella et al. , 2009)).The presence of a narrow Lyα, as well as a broader emission of Nv λ1240 Å (FWHM ∼ 2000 km s −1 ) (Vanzella et al. 2006(Vanzella et al. , 2009;;Coppin et al. 2009) identify the target as an AGN.This was furthermore confirmed by the detection of X-ray emission which resulted in an estimation of the column density of N H = 17.0 +11.7 −6.8 × 10 23 cm −2 (Circosta et al. 2019) that implies the presence of a Compton thick AGN (Gilli et al. 2011;Vito et al. 2013;Gilli et al. 2014;Circosta et al. 2019).The estimated intrinsic luminosity in the 2-10 keV band is 1.3× 10 44 erg s −1 (Luo et al. 2017) with the AGN bolometric luminosity of 2.66 ± 0.80 × 10 12 L ⊙ estimated via SED fitting (Circosta et al. 2019) classifying the target as a low-luminosity obscured QSO.The kinematic of the host galaxy was extensively studied through high-angular resolution observations of the [Cii]158µm emission line, revealing a massive, dusty disk with ordered rotation and low levels of turbulence (random, noncircular motion of the gas) (De Breuck et al. 2014;Lelli et al. 2021).The kinematics indicate the presence of a stellar bulge, which is a sign of an already-evolved galaxy.The fast evolution of the galaxy is also supported by metallicity measurements computed by exploiting the ratio between the FIR lines [Cii]158µm and [Nii]205µm that report an estimated gas-phase metallicity close to solar but we note that that the measurements have large uncertainties and this tracer of metallicity is less reliable than other optical diagnostics (Nagao et al. 2012;De Breuck et al. 2014).Gilli et al. (2014) suggest the presence of an outflow because of the observed ∼ 400 km s −1 velocity shift between the Lyα emission and other observed lines at submillimeter wavelengths ([Cii], CO(2-1), [CI], CO(7-6), [Nii]λ205 µm ;De Breuck et al. 2011;Nagao et al. 2012;D'Amato et al. 2020).

JWST observations and data reduction
The target was observed on 12 September 2022, as part of the NIRSpec GTO "Galaxy Assembly with NIRSpec IFS survey (GA-NIFS)," under the program 1216 "Integral Field Spectroscopy in GOODS-S" (PI: Nora Lützgendorf).The observations were carried out using a medium cycling pattern between four dithers to achieve a total exposure time of ∼ 1h with both G235H/F170LP and G395H/F290LP gratings/filters combinations to target the emission of the galaxy from 1.7 to 5.2 µm at high spectral resolution (on average R∼2700).The raw data were downloaded from the MAST archive and then processed with a modified version of the JWST Science Calibration Pipeline version 1.8.2 with the CRDS context "jwst_1068.pmap."The count rate maps were created by correcting at the detector level by using the module Detector1Pipeline of the pipeline.The calibra-tion was performed by applying the Calwebb_spec2 stage of the pipeline.Finally, the cube was created by adding the individual calibrated images with a drizzle weighting and a spaxel size of 0.05 ′′ using the Calwebb_spec3 step of the pipeline.The field of view of NIRSpec IFS is shown as white dashed square in Figure 1.Several corrections were made to the pipeline steps to allow for better data quality and to correct known bugs in the pipeline.These corrections are presented in detail in Perna et al. (2023), but here we mention only the major changes.The 1/ f correlated noise was subtracted from the count-rate images.The rejection of outliers was performed by using an algorithm similar to lacosmic (van Dokkum 2001) that removes the outliers on individual 2-d exposures before combining them to create the final data cube.

ALMA observation and data reduction
In this work, we also use the Atacama Large Millimeter Array (ALMA) high-resolution observation of the [Cii] emission line.ALESS073.1wasobserved with ALMA Band-7, with a total integration time of ∼ 1.5 h, to target the [Cii] emission line.The maximum and minimum baselines are 2516 and 14 m, respectively.We retrieve the [Cii] raw data from the ALMA archive (2017.1.01471.S, PI: Lelli), and then we use the pipeline scripts included in the datasets to calibrate the visibilities with the Common Astronomy Software Application CASA (McMullin et al. 2007).
Using the CASA task tclean, we perform the cleaning on the calibrated visibilities with a natural weighing scale and a spaxel scale of 0.05 ′′ to create the final datacube.The resulting datacube has a beam size of 0.17 ′′ × 0.14 ′′ .

Spectral fit of the central region
We first analyzed the emission from the nuclear region in which we expect to find emission features associated with the hidden AGN. Figure 2 illustrates the G395H/F290LP spectrum extracted from a circular aperture of radius 0.15 ′′ and centered on the central region.The errors on the spectrum were initially computed by summing in quadrature the noise from the error extension present in the data cube ("ERR") in the spaxels of the selected region.To take into account the spatial correlations of the noise due to the PSF size, we scaled the errors to match the standard deviation in line-free regions of the spectrum (see also Übler et al. 2023).In the spectrum, we clearly see the narrow (FWHM < 100 Å) emission lines of Hα and [Nii] doublets.In addition to the narrow components, it is evident the presence of a broad (FWHM > 1000 Å) line associated only with the permitted line of Hα that is consistent with being emission from the Broad Line Region (BLR).This confirms the presence of an AGN in ALESS073.1 and identifies the target as a Type 1 AGN, in contrast to what was expected due to the high value of N H .
The mismatch between X-ray and optical classification, despite not being so common, has been observed in both low and high-z AGNs (Merloni et al. 2014;Ordovás-Pascual et al. 2017;Circosta et al. 2018;Shimizu et al. 2018;Kamraj et al. 2019) and it was already proposed for this target by Gilli et al. (2014).In particular, ∼ 10-23% of optically classified Type-1 AGNs show X-ray absorbed spectra with a high value of N H (Scott et al. 2012;Ordovás-Pascual et al. 2017).This can be ascribed to different reasons: i) massive, metal-rich, high redshift galaxies can have a high ISM column density that can be large enough to Fig. 2. Spectrum in the central region from the G395H/F290LP cube.In the upper panel, we report in blue the spectrum extracted from a circular aperture of radius 0.15 ′′ centered in the central region with the associated error (gray-shaded region).We show the wavelengths around the Hα complex.The solid black line is the best-fit model resulting in the sum of the dashed lines.Each dashed line represents the best-fit result of each Gaussian component or the best-fit continuum emission.In red, the emission lines associated with the narrow component tracing the host-galaxy, in green the broader component tracing the outflow, in blue the BLR, and in dark blue the best-fit polynomial continuum.The solid vertical lines on the top represent the expected position of the [Nii] Hα and [Sii] lines.In the lower panel, we report in as a solid gray line the residuals of the fit and as a gray-shaded region the errors associated with the data.absorb X-ray radiation (Gilli et al. 2014(Gilli et al. , 2022)).If the dust-togas ratio is low or the dust grains are large (>0.03µm), the ISM dust attenuation is not sufficiently high to absorb UV and optical radiation resulting in an optically classified Type 1 AGN, with high X-ray obscuration (Maiolino et al. 2001a,b); ii) the line of sight of our observation allows us to observe the BLR unobscured from the torus, while the corona, being smaller in size than the BLR, is obscured (Shimizu et al. 2018;Kamraj et al. 2019).
To model the emission in the central region, we performed a least-square fitting by modeling the emission-line spectrum as a sum of Gaussian profiles.The Hα line profile required the addition of a broad (FWHM > 5000 km s −1 ) Gaussian component to reproduce the BLR emission.The continuum emission was fitted with a power-law function.The narrow components of both Hα and [Nii] were modeled with two Gaussian profiles each to take into account the presence of ionized outflows because using a singular Gaussian component was not sufficient to reproduce their asymmetric profiles.To disentangle between the outflow (broad) and galaxy (narrow) components of Hα and [Nii] doublets, we allowed the width of the narrow component to vary between 0 < σ < 250 km s −1 , while the outflow line width was free to vary between 250 km s −1 < σ < 1000 km s −1 .
For each Gaussian component, we tied the centroid and line width of the [Nii] doublet to those of Hα.The two emission lines of the [Nii] doublet, originating from the same upper level, were fitted with an intensity ratio I(6584)/I(6548) fixed at ∼2.94 (Storey & Zeippen 2000).Finally, the model spectrum, obtained by combining all components, was convolved with a Gaussian kernel with a dispersion of 49 km s −1 to reproduce the line spread function of the instrument at the Hα wavelength.
In Figure 2 we show the best-fit profiles for each component with dashed lines, while the best-fit values for the fluxes and the FWHMs are reported in Table 1.
We also analyzed the spectrum from the same aperture in the G235H/F170LP cube, which covers the rest-frame wavelengths from ∼ 2900Å to ∼ 5500Å.The spectrum around the Hβ -[Oiii] complex is shown in Figure 3.We did not identify any clear emission line, except for a tentative detection of the [Oiii]λ5007Å with a S /N = 2, the low S/N of the line is possibly due to high dust extinction at bluer wavelengths.We performed a single Gaussian fitting to reproduce the [Oiii]λ5007Å emission because the low signal-to-noise ratio (S/N) of the spectrum does not allow us to perform a multiple Gaussian fitting as performed for the Hα complex.We left the width of the line free to vary between 20 and 500 km s −1 to allow for the possible presence  1 together with the upper limit on the Hβ flux derived by assuming an FWHM as large as that of Hα.
By using the ratio between the Hα flux and the upper limit on the Hβ flux, we measure a lower limit on the Balmer decrement of F Hα,narrow+outflow /F Hβ > 3.6.This value is higher than the theoretical value for star-forming galaxies assuming case B recombination F Hα /F Hβ = 2.86 (Osterbrock & Ferland 2006), implying dust absorption as expected for SMG population.Assuming a Calzetti et al. ( 2000) curve we estimate a lower limit on the extinction A V > 0.77 , but we expect much higher extinction as found for other SMG (A V ∼ 4, Álvarez-Márquez et al. 2023).
To obtain a value of A V , and estimate the properties of the galaxy, we performed a new SED fitting to disentangle between galaxy and AGN contribution.The SED fitting of this target was already presented in Circosta et al. (2019) and Gilli et al. (2014), but thanks to the JWST rest-frame optical observation we are now able to characterize the AGN as a type 1, hence using the correct AGN templates.The analysis will be presented in Circosta et.al. in preparation, but here we summarize in Appendix A the methods and report the results relevant to this work.From the SED fitting we obtain a color excess E(B − V) = 0.55 ± 0.07.By using the Calzetti et al. (2000) reddening curve, we thus calculate a dust extinction A V = 2.2 ± 0.3 that is consistent with the lower limit A V > 0.77 determined by the Balmer decrement measurement.

Spatially resolved emission
We performed a spaxel-by-spaxel fitting of the G395H/F290LP data cube by exploiting the model adopted to reproduce the spectrum from the central region.We allowed the spectral components to vary except for the BLR Hα, which is spatially unresolved, and thus its centroid and FWHM were fixed to the bestfit results obtained from the analysis of the circular aperture presented in the previous section.
For each spaxel, two alternative models were adopted for the narrow Hα and [Nii] profiles: one with the outflow component and the second without it.We then selected the most suited model for each spaxel based on the Bayesian information criterion1 (BIC) test (Liddle 2007).For each spaxel, we estimated the BIC, and in those cases where the difference between BIC without outflow and BIC with outflow was larger than 2, we selected the model with two components as having a ∆BIC> 2 is consid-    1) and the flux map of the BLR component.The [Nii]λ6584Å-host emission is predominant in the central region reaching flux values comparable with the Hα-host line (see also Sec. 4).The Hα emission extends to a larger distance from the center compared to the [Nii].Based on the BIC test, the additional second Gaussian profile is necessary only for the central region of the galaxy whose size is comparable to the PSF FWHM (see Appendix B).This indicates that the outflow emission is only marginally resolved by JWST and the region directly affected is limited to the central 1 kpc.

Morphology of the host galaxy
We estimated the size of the Hα emission directly from the flux map obtained by collapsing the data cube in the wavelength range 3.775 -3.781 µm that covers the FWHM of the Hα narrow emission (Figure 2) as the flux map of the narrow component created with the pixel-by-pixel fitting (see Section 3.2) has large uncertainties due to the low S/N.
We thus performed a two-dimensional multicomponent photometric decomposition of the map.In particular, we used a 2D Gaussian profile to reproduce the emission from the unresolved BLR, the marginally resolved outflows (see Figure 4) and also taking into account the possibility of an unresolved bulge (Lelli et al. 2021 find that the bulge size is less than 300 pc, hence it is unresolved in our observations) and a 2D Sérsic profile (Sérsic 1963) with index equal to 1 (exponential disk) to describe the emission arising from the galactic disk.We also added a 2D constant to account for a possible residual of background emission.The combination of the three models was then convolved with a Gaussian point spread function of FWHM ≃ 0.202 ′′ × 0.167 ′′ obtained from the BLR flux map (see Appendix B for a detailed analysis).
The fit was carried out by using Dynesty (Speagle 2020), a Dynamic Nested Sampling Python code that allows us to estimate the Bayesian evidence and the posterior distribution of the free parameters.For the 2D Gaussian model we assumed as free parameters the centroid position, the standard deviation along the RA and Dec directions, and the amplitude.For the 2D Sérsic model the free parameters are the normalization constant, the position of the center, the effective radius, the ellipticity, and the position angle.
The data and best-fitting model are shown in Figure 5, while the corner plots for the posterior distributions of the free parameters and their best-fitting values are shown in Figure C.1.The majority of the observed flux is coming from the AGN emission (BLR+outflows) with a 2D Gaussian size of ∼ 0.01 ′′ that is five times smaller than the size of one spaxel and thus consistent with the emission from a point-like source.The effective radius of the disk is r e = 0.460±0.012′′ corresponding to an exponential scale radius of r D = 1.8 ± 0.5 kpc2 .The Hα disk scale radius is comparable within the errors with the disk radius measured with the [Cii] emission line (r D = 1.4 − 1.2 kpc De Breuck et al. 2014;Parlanti et al. 2023).We discuss the similarities and differences between the two tracers and the origin of such discrepancies in Section 8.
The ellipticity of the galaxy is 0.09±0.02corresponding to an inclination angle of 24 ± 3 deg3 assuming an infinitely thin disk.We note that the centroids of the two components have different positions (see best fit results in Figure C.1), with the Sérsic component being shifted northward by 0.9 kpc with respect to the BLR+outflow flux map centroid.This offset is probably caused by a central concentration of dust in the galaxy that absorbs Hα emission from the core (see Lelli et al. 2021, Fig. 1a).

Excitation mechanism
The rest-frame optical emission lines can be used to characterize the primary source of excitation of the gas in the ISM of galaxies.We thus exploit the Baldwin-Phillips-Terlevich (BPT) diagram (Baldwin et al. 1981), [Oiii]λ5007Å/Hβ versus [Nii]λ6584Å/Hα, to determine the dominant source of ionizing radiation and distinguish the regions mainly excited by young stars from those where the ionization mechanism is dominated by AGN radiation.
Figure 6 shows the line ratios of the narrow component for the integrated nuclear 1D spectrum, whose line fluxes are reported in Table 1.Since we do not detect the Hβ line, we can only report a lower limit on [Oiii]λ5007Å/Hβ line ratio of 3.50.This is sufficient to conclude that the gas in the nuclear region of the galaxy is excited by the AGN radiation (Kauffmann et al. 2003;Kewley et al. 2013).For the outflow component, we have upper limits for both [Oiii] and Hβ and so we can report only a vertical line at the location of [Nii]λ6584Å/Hα line ratio in the BPT diagram.The line ratio indicates that the photoionization by the central AGN dominates the gas ionization with no obvious sign of a major contribution from young stars.This supports the fact that the outflowing gas is likely accelerated by AGN radiation.The kinematics of this gas will be discussed in detail in Section 6.
We also investigate the spatially resolved excitation properties of the ionized gas but, since we do not detect Hβ and [Oiii] in the individual spaxels of the data cube, we can only probe the [Nii]/Hα line ratio.In the upper left panel of Figure 7 we report the spaxel-by-spaxel value of log([Nii]/Hα).We observe a gradient from the central regions to the outskirts of the galaxy where the values of log([Nii]/Hα) decrease at increasing radii.Since we only have this diagnostic to infer the excitation mechanism of the ISM in the galaxy, we define the three following possible ranges of log([Nii]/Hα) based on the demarcation lines suggested by Kewley et al. (2001Kewley et al. ( , 2013) ) and Kauffmann et al. (2003): region likely dominated by SF excitation (log([Nii]/Hα) ≤ −0.25), region photoionized by AGN (log([Nii]/Hα) > 0.25), and composite region −0.25 < log([Nii]/Hα) ≤ 0.25).In the BPT diagram in Figure 6 we show the number of spaxels as a function of log([Nii]/Hα) of the host galaxy component color-coded according to the aforementioned three categories, while in the lower left panel of Figure 7 we report their spatial distribution.
The majority of the spaxels show a line ratio consistent with "composite" excitation and only a few spaxels have high enough [Nii]/Hα flux to end up in the AGN region, but they reside in the outer region of the galaxy where the S/N is lower, also they are nonadjacent, so consistent with being due to S/N fluctuations shifting the category from composite to AGN.At large radii from the center, there is a number of spaxels with low (< 0.6) [Nii]/Hα suggesting that in these regions the excitation mechanism is likely dominated by star-formation activity.In conclusion, the spatially resolved BPT diagram indicates that the excitation mechanism is due to both an AGN and a young stellar population.Based on the results obtained from the nuclearintegrated 1D spectrum, we speculate the central part of the galaxy is mainly AGN-dominated while the excitation by young stars dominates at large radii.On the other hand, if we compute [Nii]/Hα spaxel-by-spaxel for the outflow component, we find that most of the spaxels have log([Nii]/Hα)> 0.25 (upper and lower right panels of Figure 7) indicating that the outflows are likely driven by AGN activity.

Black hole properties
Assuming that the gas in the BLR is virialized, we can estimate the BH mass by using the calibration by Greene & Ho (2005): L Hα 10 42 erg s −1 (0.55±0.02)

×
FWHM Hα 10 3 km s −1 (2.06±0.06)where L Hα and FWHM Hα are the dust-corrected luminosity and the FWHM of the broad Hα line associated with the BLR.
We note that the lack of detection of the Hβ line does not allow us to correct the Hα for the dust extinction of the galaxy and AGN torus.Hence the inferred luminosity for the BLR component is only a lower limit and consequently, we can only determine a lower limit on the BH mass: log (M BH /M ⊙ ) > 8.7.We can also compute a limit on the Eddington luminosity by using: where m p is the proton mass and σ T is the Thompson scattering cross-section.We obtain an Eddington luminosity of L Edd >  Appendix A), we obtain an Eddington ratio of λ Edd < 0.18 implying that the BH is accreting at a much lower rate than the Eddington limit.This is consistent with an evolutionary sequence where obscured AGN, like ALESS073.1, are in an early phase of QSO evolution that lasts until the Eddington ratio reaches values close to unity and the AGN radiation is able to sweep away gas and dust from the galaxy, revealing the emission of the bright unobscured AGNs (i.e, blue QSOs).According to this evolutionary path, we might conclude that the outflow in ALESS073. 1 is not yet energetic enough to affect the star-formation activity in the galaxy.We further investigate this possibility in Section 8.

Outflow
In this Section, we study the mass outflow rate and the energetics of the warm (T ∼ 10 4 K) ionized gas traced by the "outflow" component of Hα and [Nii] identified as the broader component associated with each line in the Gaussian fit. Figure 9 shows the kinematic maps of the outflow.The maps reveal a pattern that is not compatible with a rotating disk.Most of the spaxels in the velocity map (left panel) show negative values that are consistent with gas approaching along the line of sight.As often reported for other AGN-host galaxies (Fischer et al. 2013;Bae & Woo 2014;Perna et al. 2017), a corresponding redshifted com-ponent of the outflow is missing in ALESS073.1,probably due to dust obscuration of the receding side of the outflow.The central panel of Figure 9 shows the v 10 map, the velocity at the 10th percentile of the outflow component in each spaxel, which is usually adopted to trace the highest-velocity blueshifted gas in the outflows.We find regions in which the gas reaches a velocity as high as v 10 = −700 km s −1 supporting the fact that this gas cannot be associated with the rotation of the disk, given that the maximum velocity of the rotating disk is 400 km s −1 (Lelli et al. 2021;Parlanti et al. 2023).
The mass of gas expelled by the outflow can be estimated by following Cresci et al. (2023) as where L Hα,outflow is the extinction-corrected Hα luminosity of the outflow and n e is the electron density of the outflow.To determine the intrinsic Hα luminosity we correct the value in Tab. 1 for the value of A V estimated from the SED fitting assuming a Calzetti et al. (2000) extinction law obtaining an extinctioncorrected Hα luminosity of L Hα,outflow = 1.6 × 10 42 erg s −1 .For the electron density, since the density-diagnostic [Sii]λλ6716,31 line doublet is not detected in our observation, we have assumed the fiducial value of 1000 cm  cm −3 , based on the outflow densities of densities measured in high-redshift galaxies (Isobe et al. 2023).We thus obtain a mass of the ionized outflow of M out = 5.1 +20 −2.5 × 10 6 M ⊙ .The ionized outflow rate is calculated as follows assuming time-averaged thin expelled shells (Lutz et al. 2020) where v out and R out are the outflow velocity and radius, respectively.We adopt the prescription by Genzel et al. (2011) to estimate the velocity of the outflowing gas that takes into consideration that the emission line from the outflowing outflow is spectrally broadened due to projection effects and the velocity of the line wing traces the velocity component of the outflow directed along the line of sight, hence tracing the outflow intrinsic velocity.Thus, we obtain: v out = |∆v host galaxy,outflow | + 2σ out = 1710 km s −1 where for the values of the velocity shift and the outflow velocity dispersion we use the values obtained from the fit of the spatially integrated spectrum in Sec 3.1.For the outflow extent, we use the half width at half maximum (HWHM) of the JWST PSF at 3.78 µm given that the surface brightness emission of the outflow component is marginally resolved.We infer a mass outflow rate of Ṁout = 11 +57 −5.5 M ⊙ yr −1 .We note that the outflow radius might be smaller than the angular resolution, therefore the reported mass outflow rate may be considered a lower limit.
We also compute the lower limits on the kinetic and momentum rate of the ionized outflow as Ėout = 1 2 Ṁout v 2 out and Ṗout = Ṁout v out , respectively.We obtain Ėout = 1.2 × 10 43 erg s −1 and Ṗout = 1.47 × 10 35 g cm s −2 with uncertainties of one order of magnitude.The outflow kinetic rate is 0.1% of the bolometric luminosity.This is ∼ 50(5) times smaller than the theoretical values ( Ėout = 0.05(0.005)Lbol ) expected for quenching massive galaxies (Di Matteo et al. 2005;Choi et al. 2012;Costa et al. 2018;Harrison et al. 2018).This can suggest that the obscured AGN has not yet reached its maximum activity.
We then calculate the mass loading factor, defined as η = Ṁout /SFR, where we adopt the SFR = 196 ± 10 M ⊙ yr −1 for ALESS073.1 obtained with a SED fitting (see Appendix A).We obtain a mass loading factor on the order of 6%, meaning that only a minor part of the gas present in the central region is expelled compared to the gas used to create stars, and the outflow strength is probably not sufficient to remove gas and halt the vigorous, ongoing SF in the system.We note that using other SFR estimations result in a low mass loading factor as well (SFR De Breuck et al. 2014).

Gas kinematics
In this Section, we investigate the gas kinematics traced by the Hα narrow component that maps warm ionized gas in the interstellar medium of the galaxy and we compare it with the results obtained by studying the [Cii] emission.The first two panels on the left of Figure 10 illustrate the velocity and velocity dispersion maps of the Hα narrow component.The velocity and velocity dispersion maps are computed as the displacement between the centroid position of the line in that spaxel with respect to the centroid position computed in the central region, and the standard deviation of the emission line deconvolved of the instrumental spectral resolution computed on a spaxel-by-spaxel level (see Section 3.2), respectively.The velocity map shows a velocity gradient that spans a range of velocities between -120 km s −1 and 120 km s −1 with respect to the systemic redshift of the galaxy.The velocity pattern is consistent with that observed in [Cii] by Lelli et al. (2021) and indicates the presence of a regularly rotating disk.We note that the velocity gradient is not symmetric on the red-shifted and blue-shifted sides, with the blue side having higher velocities closer to the center with respect to the red-shifted one, and the blue-shifted side being less extended.On the northwestern side, we also note that the velocity pattern is irregular, similar to what is found in the [Cii] kinematic analysis (Lelli et al. 2021) that might indicate the presence of a spiral arm.
Following Parlanti et al. (2023), we model the Hα kinematics using the publicly available python library KinMS (Davis et al. 2013) that creates mock data cubes based on flux, velocity, and velocity dispersion radial profiles.We set up KinMS to simulate our NIRSpec observations, setting an angular resolution of 0.202 ′′ × 0.167 ′′ (see Appendix B) and a spectral resolution of 49 km s −1 .We then generate the moment maps to com-  2021) find indeed that the de-convolved velocity is as high as ∼ 400 km s −1 in the nuclear region (< 1 kpc) of the galaxy which is comparable with what is observed in the outskirt of the galaxy.These kinematics cannot be reproduced by the velocity curve of an exponential disk, as we show and discuss in Appendix C. The flat profile of the [Cii] velocity indeed suggests the presence of a bulge component in addition to the classic exponential kinematic component (Lelli et al. 2021;Parlanti et al. 2023).
Thanks to the flexibility of KinMS, we adopt a nonparametric model to reproduce the observed velocity and velocity dispersion profile as a function of the radius.In particular, the Hα kinematics is modeled as a series of concentric circular rings having a width of 0.6 kpc (i.e., 0.09 ′′ ), which is comparable to the HWHM of the PSF of NIRSpec IFS at the Hα wavelength.In each ring, we assume that the emitting clouds have the same radial velocity and velocity dispersion.For the velocity fitting we assume Gaussian priors with a standard deviation of 100 km s −1 centered around the velocity of the previous ring.The velocity of the first ring was left free to vary between 0 and 1000 km s −1 .Similarly, for the velocity dispersion of each ring, we assume Gaussian priors with a standard deviation of 50 km s −1 centered on the velocity dispersion of the previous ring leaving the first ring free to vary between 0 and 500 km s −1 .Using Gaussian priors allows us to ensure the continuity of the velocity and velocity dispersion profiles and that the discontinuities in the profiles are driven by a real increase in the likelihood.We note that assuming flat priors for each ring does not change the maps in Figure 10, but introduces discontinuities in the intrinsic profiles, that are then washed out by the beam smearing process when making the moment maps.
We let the disk position angle free to vary with flat priors between 0 and 90 deg, while we fixed the inclination angle of the galaxy to 22 deg as found by Lelli et al. (2021) and compatible within the errors with the value of inclination found in Section 3.3.The emitting clouds are distributed over the rings following the surface brightness profiles obtained from the best-fit results of the flux map (Section 3.3).We note that the resolution and sensitivity do not allow us to determine stringent constraints on the velocity and velocity dispersion profiles if we adopt a disk model with a ring size smaller than the PSF HWHM due to the beam smearing.We also observe that the last ring is only probed by the outer region of the redshifted side of the galaxy as it is more extended than the blueshifted one.
The best-fit results are reported in Table 2 and the best-fit model and residual maps of the velocity and velocity dispersion are reported in Figure 10.In Figures C.2 and C.3 we also report the corner plot to highlight the best-fit parameters and the posterior distributions that show degeneracy between the parameters.We note that the velocity and velocity dispersion of each ring degenerate with the ones of the previous ring.In particular if one increases the other one decreases, and vice versa.This is expected and consistent with the beam smearing that acts to average out these differences.The inferred disk position angle is in agreement with the one found for the [Cii] kinematics by Lelli et al. (2021).
The left panel of Figure 11 shows the best-fit velocity curve of Hα, which is in agreement, within the errors, with the profile inferred from [Cii] observation.The velocities reach a value of ∼ 350 km s −1 at small scales and slightly decrease at large radii down to values on the order of ∼ 250 km s −1 .This result supports the scenario that ALESS073.1 has already formed a bulge at its center that dominates the dynamics of the gas.
The velocity dispersion profile is reported in the right panel of Figure 11 and spans a range between 30 and 170 km s −1 .Differently from the velocity curve, here we note a discrepancy between [Cii] and Hα gas kinematics.The velocity dispersion of the warm ionized gas is, on average, higher than the velocity dispersion of the cold gas mapped by the carbon line.In the in- ner 1.5 kpc the velocity dispersion observed in Hα is even three times larger than the one traced by [Cii].At larger radii, the discrepancy decreases and in the two rings between 1.5 and 3.0 kpc the velocity dispersion of [Cii] is comparable (within the uncertainties) with the one derived from Hα.
As the discrepancies between our results and those found by Lelli et al. (2021) might depend on the different approaches and tools used to fit the data, we compare directly the velocity dispersion maps in the same regions of the galaxy.As the ALMA and NIRSpec data cubes of ALESS073.1 have a slightly different size PSF (ALMA: 0.17 ′′ × 0.14 ′′ , JWST: 0.20 ′′ × 0.17 ′′ ) we used the Python library Photutils to create a matching kernel between ALMA and JWST PSFs.Hence we convolved the ALMA cube, which is the higher resolution one with the matching kernel to obtain 2 data cubes with the same PSF.After this process the two data cubes are affected by the same level of beam smearing, enabling us to compare the velocity dispersion maps directly.The [Cii] moment maps are created by fitting spaxel-byspaxel singular Gaussian component to the ALMA spectrum.We also try to include a second component to verify the presence of outflows in [Cii] but the double Gaussian model returns a higher χ 2 in all spaxel and the BIC test supports the single component fit.
Figure 12 illustrates the velocity dispersion of the narrow Hα component and [Cii] emission line with the same velocity range.We stress that the line broadening due to the line spread function of NIRSpec is corrected during the Hα line fitting (Sec-tion 3.2) and the Hα velocity dispersion map does not include the broad blue-shifted Hα component associated with outflows.Comparing the kinematic maps, we find that the velocity dispersion mapped by the hydrogen line is on average higher than that traced by [Cii].Both maps have high σ at the center (σ [CII] ∼ 130 km s −1 , σ Hα ∼ 200 km s −1 ) mainly due to beam smearing, but as shown in the σ Hα − σ [CII] map (right panel) the discrepancy between the two increases at large radii up to 0.3-0.4′′ .In particular, the difference in terms of velocity dispersion in the region is evident in the direction north-south from 0.1 ′′ to 0.4 ′′ from the center.In these regions ⟨σ Hα − σ [CII] ⟩ = 100 km s −1 that is more than two times larger than the uncertainties on the velocity dispersion estimates.

Discussion
Investigating galaxy dynamics at high redshifts is fundamental to understand how galaxies grow their stellar mass at their early stages of formation.Measurements of small ratios between the rotational velocity and the gas velocity dispersions are usually interpreted as evidence for turbulence in the disk due to past or ongoing strong feedback mechanisms and merging processes.Conversely, observations revealing a low velocity dispersion and V/σ ∼ 10 in the gas kinematics suggest a less turbulent gas accretion and evolution characterized by a limited number of extreme events.The trends of velocity dispersion and the V/σ evolution with redshift at z > 4 are still a matter of debate today.On one hand, a large number of kinematic studies at 1 < z < 4 suggest that high redshift galaxies are more turbulent than local ones (Cresci et al. 2009;Förster Schreiber et al. 2009;Epinat et al. 2010;Gnerucci et al. 2011;Ianjamasimanana et al. 2012;Green et al. 2014;Wisnioski et al. 2015;Mogotsi et al. 2016;Di Teodoro et al. 2016;Harrison et al. 2017;Swinbank et al. 2017;Turner et al. 2017;Förster Schreiber et al. 2018;Johnson et al. 2018;Übler et al. 2019;Girard et al. 2021) with V/σ reaching values close to unity at z ∼ 3.5.On the other hand, the results from kinematic studies at z > 4 lead to contrasting results showing the presence of both turbulent (Tsukui & Iguchi 2021;Herrera-Camus et al. 2022;Parlanti et al. 2023;de Graaff et al. 2023) and kinematically cold galaxies (Sharda et al. 2019;Neeleman et al. 2020;Rizzo et al. 2020;Jones et al. 2021;Fraternali et al. 2021;Lelli et al. 2021;Rizzo et al. 2021;Posses et al. 2023;Pope et al. 2023) with values of V/σ ranging from 20 to 0.1 across the redshift range 4 < z < 8.However, we note that results at high redshift are limited either by small sample sizes for galaxies observed with high angular resolution or by significant uncertainties due to the low angular resolution observations used for the larger samples of galaxies.
We also note that so far galaxies up to z ∼ 3.5 are principally studied by exploiting Hα and [Oiii] emission lines that arise from HII regions around massive and young stars and trace the warm ionized medium with a temperature of ∼ 10 4 K.Studies at z > 4 target mainly the [Cii] line that principally arises from photodissociation regions (PDR) tracing the cold neutral medium at a temperature of ∼ 100 K and only ∼ 30% of its emission is associated with the warm diffuse ionized gas (e.g., Stacey et al. 1991Stacey et al. , 2010;;Croxall et al. 2017).Therefore, several studies have concluded that the cosmic evolution of the velocity dispersion obtained from the optical lines cannot be directly compared with that from the far-infrared [Cii] line because they are mapping different gas phases (Rizzo et al. 2022).One possible solution is to exploit the far-infrared [Oiii]λ88µm line that traces the warm ionized medium and can be observed with ALMA at z > 6.Unfortunately, long exposure times are necessary to obtain an accurate measurement of the gas velocity dispersion, and most of the current observations have coarse angular resolution resulting in large uncertainties (Parlanti et al. 2023).Moreover, [Oiii]λ88µm would still leave a gap in observations in the redshift range 3.5 < z < 6.
With the advent of JWST NIRSpec and in particular, thanks to the IFS mode observations, we can finally compare the kinematics determined from the rest-frame optical lines with that of [Cii] in the same z > 4 galaxies and verify if there is a discrepancy between these tracers or not.ALESS073.1 is the first massive (M ⋆ ∼ 10 11 M ⊙ ) galaxy for which we have both ALMA and NIRSpec high-resolution observations and the results presented in Section 7 show that the velocity curves of the two tracers are consistent within the errors.The data however highlight a difference between Hα and [Cii] in terms of velocity dispersion.The velocity dispersion inferred from Hα is systemically larger by more than 50 km s −1 in the central 1.5 kpc than the one determined from the carbon line.This difference cannot be associated only with the difference in thermal velocity dispersion of the gas mapped by the two tracers as this is on the order of 20 km s −1 .The difference between the two gas kinematics is not uniform over the field of view but is higher at the center ⟨σ Hα − σ [CII] ⟩ = 100 km s −1 and reaches values comparable with the errors (∼ 30 km s −1 ) at larger radii.Significant differences between the Hα and [Cii] kinematics have been also recently reported both from observations (Arribas et al. 2023) and simulations (Kohandel et al. 2023) To understand whether the galaxy is rotation or dispersionsupported we compute the ratio between the rotational velocity and the velocity dispersion.When the ratio is larger than √ 3.36 the galaxy is considered supported by rotation, on the contrary, is supported by the turbulent random motion of the gas (Förster Schreiber et al. 2018).In Figure 13 we report the ratio V/σ computed for each ring for the Hα (red) and [Cii] (blue) as small symbols along with other high and low redshift galaxies.The big blue triangle and red circle represent the mean value of V/σ across all the rings for the [Cii] and Hα profiles, respectively.If we compare the value we obtain with the literature we see that the results with both tracers lie in the rotation-dominated region, even if the Hα and [Cii] points create two different clouds around the central value of V/σ = 10 for the [Cii], and V/σ = 3.5 for the Hα.The ratio between the rotational velocity and the velocity dispersion derived from the Hα line is V/σ ∼ 10 in the outer regions as also found with the [Cii] tracer for this galaxy (Lelli Fig. 13.Evolution of the ratio between the rotational velocity and the velocity dispersion with redshift.Red and blue small symbols represent the value of V/σ for each ring (see also Figure 11) for Hα and [Cii], respectively.While large symbols are the mean V/σ value across the galaxy with the two tracers.Gray symbols are other literature results, in particular, circles represent the results from kinematic studies that exploited ionized gas tracers (i.e., Hα, [Oiii]) (Green et al. 2014;Turner et al. 2017;Förster Schreiber et al. 2018;Wisnioski et al. 2019;Parlanti et al. 2023;de Graaff et al. 2023), while triangles represent the V/σ values for galaxies studied though the molecular or neutral emission lines (i.e., [Cii], CO) (Rizzo et al. 2020;Fraternali et al. 2021;Jones et al. 2021;Girard et al. 2021;Lelli et al. 2021;Tsukui & Iguchi 2021;Rizzo et al. 2021Rizzo et al. , 2023)).The dashed line is the demarcation between rotational supported (upper) and dispersion supported (lower).et al. 2021) and other high-redshift dusty star-forming galaxies (e.g., Rizzo et al. 2020Rizzo et al. , 2021;;Fraternali et al. 2021), and decreases down to values of ∼ 2 in the central region, consistent with a thick turbulent rotating disk.
The enhanced velocity dispersion of the warm ionized gas can be caused by the outflows driven by the central AGN.Theoretical models predict that outflows might remove gas from the galaxy as well as inject energy into the interstellar medium and kinematically heat the gas.The increase in the turbulence in the ionized phases of the ISM due to the feedback effect has been recently observed in Marasco et al. (2023) at z = 0, while a correlation between the increase of the turbulence in the galaxy and the presence of a central AGN that can power a nuclear outflow has been found in Übler et al. (2019).The Hα gas might be affected by the galactic winds and its kinematics reflects the impact of the outflows on the host galaxy.We thus compare the kinetic energy of the gas with that injected by the outflows into the galaxy.We estimate the kinetic energy of the warm gas traced by the narrow Hα component, E Hα = M Hα σ 2 Hα /2 = 1 × 10 54 erg, by employing the Hα flux and velocity dispersion in the region with σ Hα − σ [CII] > 70 km s −1 that is the region where the difference between the velocity dispersions in the third panel of Figure 12 is two times larger than the median error.We compare the inferred kinetic energies with the outflow energy, E out = ⟨ Ėout ⟩τ out = 2.3 × 10 57 erg, where we assume that the outflow kinetic rate is constant over time and the feedback mechanism started 10 Myr ago (τ out = 10 Myr).We find that the ratio between the energy provided by the outflow, and the energy necessary to increase the turbulence of the ISM is E out /E Hα = 0.05%.The result suggests that the energy of the outflows is powerful enough to provide the kinetic energy of the ionized gas observed in the line-width enhancement region.In the outer regions of the galaxy, instead, the warm and cold gas are coupled as the velocity dispersion of both tracers reaches values of 30 km s −1 and we do not find evidence of the feedback from the weak outflow or the accreting BH enhancing the turbulence at larger radii.The level of turbulence in the outskirts of the galaxy can easily be sustained by star-formation feedback, gravitational instabilities due to the accretion of gas on the disk or due to the transport of gas from outer to inner radii (Krumholz & Burkhart 2016;Krumholz et al. 2018;Ginzburg et al. 2022).

Conclusions
In this work, we have presented the JWST/NIRSpec Integral Field Spectrograph (IFS) observation of the AGN-host galaxy ALESS073.1 at z = 4.755.The observations of the highresolution gratings have allowed us to study for the first time the rest-frame optical emission lines of a dusty-obscured SMG hosting an AGN.In particular, we exploited the Hα and [Nii] emission lines to trace the host galaxy kinematics, determine the presence of a BLR, and investigate the properties of ionized outflows.Our main results are the following: -We observe a broad Hα component with a FWHM of ∼ 9000 km s −1 arising from the BLR around the accreting supermassive black hole.The presence of the BLR unambiguously classifies the target as a type 1 AGN, in contrast to what was believed before due to the high observed column density of N H ∼ 10 24 cm −2 .The broad line emission implies a BH mass of log(M BH /M ⊙ ) > 8.7 that is slowly accreting at a smaller rate than the Eddington rate (λ Edd <0.18).-On the M BH − M ⋆ plane, the target lies on the relation for massive quiescent local ellipticals, classic bulges, and luminous QSO at high redshift.But M BH /M ⋆ is more than one order of magnitude higher with respect M BH /M ⋆ observed in local AGN with similar stellar masses.-We find hints of a weak, marginally resolved, ionized outflow with a mass loading factor of ∼0.06, implying that the outflow is not able to eject away a large amount of gas to halt the SF ongoing in the galaxy.-By measuring the ratio between [Nii] and Hα we have found that the AGN hard radiation is the dominant source of ionization of the ISM, especially in the central region of the galaxy and for the outflow component.At larger radii we find, instead, softer radiation, compatible with emission from young, bright stars.-Despite the low mass loading factor, the ionized outflow seems to be sufficiently energetic to increase the turbulence in the system.In fact, the kinematic analysis of the Hα line shows that gas turbulence in the central region is 2-3 times higher than the rest of the galaxy.This increase in turbulence might be the initial effect of the outflow on the host galaxy.
However, as such high-velocity dispersion is observed only in the Hα kinematics and not in the map of [Cii], we conclude that the outflow is injecting turbulence in the warm and diffuse ionized gas, but it is not sufficiently powerful to disrupt the dense gas and quench star formation.
The complex scenario of galaxy-black hole coevolution is still far from being constrained and firmly established.In this work, we have highlighted how JWST with its high spatial resolution, spatially resolved spectroscopy capability, and an infrared wavelength range probed with unprecedented sensitivity, significantly enhance previous studies based on observations at various wavelengths.In particular, its ability to probe the rest frame optical emission lines at high redshift with high spectral and spatial resolution has allowed us to study the first phases of the interplay between the accreting black hole and the host galaxy and connect kinematic measurements of high-z tracers ([Cii]) with the more traditional rest-frame optical emission line tracers at lower redshift.Future JWST IFS observations, alongside other groundbased facilities (e.g., ALMA) will allow us to better observe and understand the phenomena involved in the complex environment of dusty massive galaxies in the early Universe and their interaction with the massive black holes they host.ALESS073.1 observations suggest that SMGs could be the evolutionary stage preceding the active QSO phase.The accreting BHs in SMG have not yet reached the Eddington limit and outflows are not powerful enough to remove gas from the galaxy, but they are injecting energy into the system and increasing the turbulence of less dense gas.(Circosta et al. in prep.).The black circles are all the available photometric points, while the white-filled circles represent the 3σ upper limits.The black solid line is the best-fit result.The yellow, red, and green solid lines are the best models for the attenuated stellar emission, dust emission, and AGN emission, respectively.
is spread over a larger area on the detector.This effect has a consequence both in the spatial and the spectral axis due to the reconstruction of the final cubes.In particular, on the moment maps the flux map will have a larger extension, the velocity map will have a shallower velocity gradient, and emission lines will be artificially broadened resulting in a larger velocity dispersion, especially in the central region of the galaxy.To correct these artifacts and recover the intrinsic model parameters we convolve our kinematic model with the telescope PSF, in that way, we can directly compare the model and observations.The webbpsf python tool (Perrin et al. 2014) generates the theoretical PSF by taking into account instrumental properties and the optics based on test results.Since JWST has wavelengthdependent PSF (see also D'Eugenio et al. 2023) we simulated the PSF at the wavelength at which we observe the Hα line (3.778 µm).The generated PSF presents spikes and irregularities and fitting it with a 2D Gaussian we obtain a FWHM of 0.11×0.10′′ .
Since webbpsf only simulates the NIRSpec PSF if used in imaging mode without taking into account the IFS optics, the dithering schemes, and the cube reconstruction processes, we also have fitted with a Gaussian function the BLR flux map.The broad line emission arises from a region close to the black hole (Wandel et al. 1999), hence it will not be spatially resolved and its shape should reflect the instrumental PSF.The results of the fitting are shown in Figure B.1.The results of the fitting provide an elliptical PSF with a FWHM of 0.202 ′′ × 0.167 ′′ where the major axis is rotated counterclockwise in respect to the x-y axis of an angle θ = 52 ± 4 deg.The resulting PSF is elongated in the cross-dispersion direction (i.e., slices direction).As the BLR emission should reflect the real telescope PSF we assumed the BLR best fitting results as our fiducial model for the PSF to use in our kinematic analysis.
Fig. 1.RGB image of ALESS073.1.The colors are a combination of HST/ACS and HST/WFC3 images.White contours highlight the position and size of the NIRSpec IFS field of view.

Fig. 3 .
Fig. 3. Spectrum of the central region from the G235H/F170LP cube.In the upper panel, we report in blue the spectrum around the [Oiii]λ5007Å emission, extracted from a circular aperture with a radius of 0.15 ′′ in the central region with the associated error (gray-shaded region).In black we report the best-fit model resulting from the sum of the dashed lines.The red dashed line represents the best-fit result of the Gaussian component tracing [Oiii]λ5007Å, the blue dashed line represents the best-fit continuum emission.The solid vertical lines on the top represent the expected position of the [Oiii] and Hβ lines.In the lower panel, we report in as a solid gray line the residuals of the fit and as a gray-shaded region the errors associated with the data.

Fig. 4 .
Fig. 4. Flux maps created from the results of the spaxel-by-spaxel multi-Gaussian fitting.Upper row: the flux of the Hα and [Nii]λ6584Å narrow components representing the host galaxy emission from left to right, respectively.Bottom left panel: the flux of the outflow component traced by the [Nii]λ6584Å broader component.Bottom right panel: the broad Hα component that traces the emission from the BLR.The x and y axes are the displacement in arcseconds from the galaxy center at RA = 03:32:29.3,Dec = -27:56:19.6.
ered marginal positive evidence in favor of the model with lower BIC value(Kass & Raftery 1995).In the other cases, we adopted the model with one component.We note that selecting the model with the lower BIC value allows us to select the best-fit model with the highest statistical significance without overfitting the data.The low S/N on a spaxel-by-spaxel analysis does not allow us to use a larger ∆BIC threshold to obtain more stringent constraints on the outflow statistical significance.Figure 4 illustrates the flux maps for the narrower component of Hα and [Nii]λ6584Å tracing the host-galaxy, the flux map of the [Nii] λ6584Å broad component associated with outflows as it is much stronger than the one traced by Hα (see Figure 2 and Table

Fig. 5 .
Fig. 5. Observed Hα flux map, best-fitting model, and residuals, from left to right, respectively.In the left panel, we show the Hα flux, overlaid with the [Cii] contours at 3, 6, and 9 σ.The two maps were aligned by centering them on the brightest spaxel.The gray ellipse represents the PSF size and shape.In the central panel, we show the best-fit model composed by the sum of a 2D Gaussian component that includes all the unresolved or marginally resolved components (BLR, outflow, bulge) and a 2D Sérsic component with the Sérsic index n = 1 that represent the galactic disk.In the right panel, we show the residuals, calculated as the observed flux minus the model, divided by the error.The color bar stretches between −3σ and +3σ.

Fig. 6 .
Fig. 6.BPT diagram of the target.The star marker in magenta shows the position of the central spaxels of ALESS073.1inthe BPT diagram for the narrow component.The position of the outflow component based only on the log([Nii]/Hα) detection is reported as a purple vertical line and the error associated with it is reported as the purples-shaded area.The dashed and the dotted lines are the predictions from Kewley et al. (2001) and Kauffmann et al. (2003), respectively, for the separation between star-forming (below) and AGN (above) dominated regions at z = 0.The dash-dotted line is the theoretical redshift evolution of the separation curve for galaxies at z = 3 by Kewley et al. (2013).The gray shaded area represents the position in the BPT of SDSS galaxies at z ∼ 0. Blue diamonds are the results from Cameron et al. (2023) for high-redshift (z > 5) galaxies.Overplotted with the BPT diagram we report the histogram representing the number of spaxels of the narrow component as a function of log([Nii]/Hα).
shows the relation between the stellar mass and the black hole mass for ALESS073.1,where we use the stellar mass of M ⋆ = (4.7 ± 1.6) × 10 10 M ⊙ computed byLelli et al. (2021) with a kinematic analysis, and M ⋆ = 9.5 +4.3 −2.9 × 10 10 M ⊙ estimated from the SED fitting (see Appendix A).We compare our results with those obtained from other AGN-host galaxies both at high redshift (z > 1)(Pensabene et al. 2020;Neeleman et al. 2020;Harikane et al. 2023;Kocevski et al. 2023;Larson et al. 2023;Maiolino et al. 2023;Übler et al. 2023) and in the nearby Universe (z < 0.055)(Reines & Volonteri 2015), and the local massive "red-and-dead" elliptical galaxies and classical bulges(Kormendy & Ho 2013).The estimated BH mass places ALESS073.1 above the relation byReines & Volonteri (2015) yielding a M BH /M ⋆ ratio > 10 − 30 times higher than those estimated in local AGNs by using the stellar mass estimated by SED and kinematical fitting, respectively.However, the inferred M BH /M ⋆ is consistent within the uncertainties with the relation determined for massive quiescent local galaxies and high redshift luminous quasars(Kormendy & Ho 2013;de Nicola et al. 2019;Pensabene et al. 2020).As with other high-redshift AGNs and QSOs, it lies above the relation M BH = 0.01 × M ⋆(Decarli et al. 2010;de Nicola et al. 2019;Pensabene et al. 2020;Neeleman et al. 2021).These high-redshift observations suggest that the BH growth dominated early on, with the galaxy catching up later.This requires that feedback and self-regulation are somehow different at early times with respect to what is observed in local AGNs.

Fig. 7 .
Fig. 7. Resolved BPT diagram for ALESS073.1 host galaxy (narrow) and outflow (broad) component.Upper panels: spatially resolved map of the observed values of log([Nii]/Hα)for the narrow and broad components from left to right, respectively.Lower panels: the spatially resolved BPT where every spaxel in the galaxy is color-coded according to its value of log([Nii]/Hα) for the narrow and broad components, from left to right, respectively.as in the histogram in the upper right panel.

Fig. 8 .
Fig. 8. Relation between black hole mass and host galaxy stellar mass.The galaxy targeted in this work (ALESS073.1)is marked with a star, the star filled in red is by assuming the stellar mass derived from the kinematic fitting by Lelli et al. (2021), while the white filled star considers the stellar mass derived from SED fitting (see Appendix A).The dark gray diamonds and crosses represent z ∼ 0 broad line AGNs presented in Reines & Volonteri (2015) and the massive black holes hosted in ellipticals and spirals at z ∼ 0 by Kormendy & Ho (2013), respectively.In light green diamond, we report the QSO at z ∼ 5.5 studied in Übler et al. (2023).In green, we show the broad line AGNs discovered with JWST at 4 < z < 7 presented in Harikane et al. (2023), Maiolino et al. (2023) and Kocevski et al. (2023) as light blue circles, triangles and pentagons, respectively.The yellow triangle represents the AGN at z ∼ 8.7 identified in Larson et al. (2023).Blue squares are the results obtained by Neeleman et al. (2021) for QSOs at z ∼ 6. Blue triangles are the results for QSOs at z > 2 from Pensabene et al. (2020).Blue crosses are the results for two luminous QSOs at z ∼ 6.8 from Marshall et al. (2023).We note that Pensabene et al. (2020), Neeleman et al. (2021) and Marshall et al. (2023) report the dynamical mass rather than the stellar mass.The black and green solid lines are the best-fit for Reines & Volonteri (2015) and Kormendy & Ho (2013), respectively, and the shaded areas are their 1σ uncertainties.The black dashed line is the relation M BH = 0.01 × M ⋆ .

Fig. 9 .
Fig.9.Spatially resolved kinematic maps for the outflow component.From left to right, the velocity map of the outflow component, the v10 that corresponds to the velocity at which the outflow component is at the 10th percentile, and the velocity dispersion map of the outflow, respectively.Zero velocity corresponds to the velocity of the narrow line.

Fig. 10 .
Fig. 10.Nonparametric best-fit results for the modeling of the Hα velocity and velocity dispersion maps.From left to right on the top the observed velocity map, the best-fit model map, and the residuals.From left to right on the bottom the observed velocity dispersion map, the best fitting model dispersion map and the residuals.The colorbars of the residual range between −3σ and +3σ, and the black lines indicate ±1σ.

Fig. 11 .
Fig. 11.Velocity (left) and velocity dispersion (right) profile derived from Hα and [Cii] from Lelli et al. (2021).Blue points are the results with the [Cii] tracer by Lelli et al. (2021).Red solid lines are the results of the nonparametric fitting of the Hα maps from this work with associated uncertainties.In gray is the region affected by the PSF.

Fig. 12 .
Fig. 12. Observed velocity dispersion maps and difference between Hα and [Cii].In the left panel the Hα velocity dispersion map.In the central panel, the [Cii] velocity dispersion map rebinned to have the same spaxel size as Hαand matched to have the same PSF.In the right panel the spaxel-by-spaxel difference between the velocity dispersion of Hα and [Cii].

1 .
Fig. A.1.SED fitting of ALESS073.1 in rest wavelengths(Circosta et  al. in prep.).The black circles are all the available photometric points, while the white-filled circles represent the 3σ upper limits.The black solid line is the best-fit result.The yellow, red, and green solid lines are the best models for the attenuated stellar emission, dust emission, and AGN emission, respectively.

Fig. B. 1 .
Fig. B.1.Results of the Gaussian fitting of the BLR to recover the observed point spread function.On the top, the BLR observed flux maps, the best-fit model, and residual from left to right.On the bottom, the corner plot reports the best-fit results and uncertainties.The amplitude is a normalization constant, x 0 and y 0 are the displacement in arcsec of the centroid from the galaxy center coordinates, σ x and σ y are the minor and major axis standard deviations in arcsec, θ is the rotation angle in radians measured counterclockwise from the positive y-axis.Above each histogram, we show the best value for each parameter and the 2σ interval.

Fig. C. 1 .
Fig. C.1.Posterior distribution for the free parameters in the fitting of the flux map.The dashed lines represent the 16th, 50th and 84th percentile.The galaxy is modeled as a 2D Sérsic profile, while the AGN is modeled as a two-dimensional Gaussian.The coordinates of the centers are expressed in arcseconds where the position (0,0) corresponds to the coordinates of RA = 03:32:29.3,Dec = -27:56:19.6.The radius of the galaxy and the standard deviation of the Gaussian component σ x and σ y are measured in arcseconds and θ gal is the angle measured in radians between the y-axis and the semi-major axis of the galaxy.

Table 1 .
Results from the fitting of the spectrum extracted from a circular aperture of radius 0.15 ′′ centered on the spatial peak of the emission.
Measurement [km/s] FWHM Hα host galaxy 519 ± 58 FWHM Hα outflow 1491 ± 206 FWHM Hα BLR 9008 ± 407 FWHM [Oiii]λ5007Å 553 ± 105 ∆v Hα host galaxy,outflow -441 ± 129 ∆v Hα host galaxy,BLR 427 ± 201 Notes.On top: FWHM of the detected lines and velocity shift between the systemic velocity identified by the narrow component of the Hα relative to the host and the centroid of the outflow (broader component) and BLR ∆v Hα host galaxy,outflow = v Hα host galaxy − v Hα outflow .On the bottom: fluxes of the detected lines.

Table 2 .
Best fitting kinematic parameters obtained by the nonparametric fitting.