11institutetext: INAF - Osservatorio Astrofisco di Arcetri, largo E. Fermi 5,
50127 Firenze, Italy
11email<EMAIL_ADDRESS>22institutetext: Università degli Studi di
Firenze, Dipartimento di Fisica e Astronomia, via G. Sansone 1, 50019 Sesto
F.no, Firenze, Italy 33institutetext: Centro de Astrobiología (CAB,
CSIC–INTA), Departamento de Astrofísica, Cra. de Ajalvir Km. 4, 28850 –
Torrejón de Ardoz, Madrid, Spain 44institutetext: Dipartimento di Fisica e
Astronomia, Alma Mater Studiorum, Universitá di Bologna, via Gobetti 93/2,
40129 Bologna, Italy 55institutetext: INAF - Osservatorio di Astrofisica e
Scienza dello Spazio di Bologna, via Gobetti 93/3, 40129 Bologna, Italy
66institutetext: Scuola Normale Superiore, Piazza dei Cavalieri 7, I–56126,
Pisa (PI), Italy 77institutetext: INAF - Istituto di Radioastronomia, via
Gobetti 101, 40129 Bologna, Italy 88institutetext: Cavendish Laboratory,
University of Cambridge, 19 J.J. Thomson Ave., Cambridge, UK 99institutetext:
Kavli Institute for Cosmology, University of Cambridge, Madingley Road,
Cambridge, UK 1010institutetext: Instituto de Astrofísica, Facultad de Física,
Pontificia Universidad Católica de Chile, Casilla 306, Santiago 22, Chile
1111institutetext: Space Telescope Science Institute, 3700 San Martin Drive,
Baltimore, MD 21218, USA 1212institutetext: INAF – Osservatorio Astrofisico di
Torino, Via Osservatorio 20, 10025 Pino Torinese, Italy

# Bubbles and outflows: The novel JWST/NIRSpec view of the z=1.59 obscured
quasar XID2028

G. Cresci 11 G. Tozzi 2211 M. Perna 33 M. Brusa 4455 C. Marconcini 2211 A.
Marconi 2211 S. Carniani 66 M. Brienza 5544 M. Giroletti 77 F. Belfiore 11 M.
Ginolfi 2211 F. Mannucci 11 L. Ulivi 2211 J. Scholtz 8899 G. Venturi 101011 S.
Arribas 33 H. Übler 8899 F. D’Eugenio 8899 M. Mingozzi 1111 B. Balmaverde 1212
A. Capetti 1212 E. Parlanti 66 T. Zana 66

(Received ; Accepted)

Quasar feedback in the form of powerful outflows is invoked as a key mechanism
to quench star formation in galaxies, although direct observational evidence
is still scarce and debated. Here we present Early Release Science JWST
NIRSpec IFU observations of the $z$=1.59 prototypical obscured Active Galactic
Nucleus (AGN) XID2028: This target represents a unique test case for studying
quasar feedback at the peak epoch of AGN-galaxy co-evolution because extensive
multi-wavelength coverage is available and a massive and extended outflow is
detected in the ionised and molecular components. With the unprecedented
sensitivity and spatial resolution of the JWST, the NIRSpec dataset reveals a
wealth of structures in the ionised gas kinematics and morphology that were
previously hidden in the seeing-limited ground-based data. In particular, we
find evidence of an interaction between the interstellar medium of the galaxy
and the quasar-driven outflow and radio jet that produces an expanding bubble
from which the fast and extended wind detected in previous observations
emerges. The new observations confirm the complex interplay between the AGN
jet, wind and the interstellar medium of the host galaxy, highlighting the
role of low-luminosity radio jets in AGN feedback. They also clearly show the
new window that NIRSpec opens for detailed studies of feedback at high
redshift.

###### Key Words.:

Galaxies: active – Galaxies: individual (XID2028) - ISM: jets and outflows –
Techniques: imaging spectroscopy

##  1 Introduction

Massive and fast outflows are almost ubiquitous in luminous active galaxies.
They are now observed in different gas phases and on different physical
scales, from parsec-scale ultrafast ($\sim$10% of the speed of light) outflows
detected in X-rays (e.g. Pounds et al. 2003, Nardini et al. 2015, Matzeu et
al. 2022) to kiloparsec-scale outflows observed in atomic, molecular, and
ionised gas with velocities up to $\sim 1000\ \rm km\ s^{-1}$ (e.g. Rupke &
Veilleux 2011, Cicone et al. 2014, Carniani et al.2015, Harrison et al. 2016,
Förster-Schreiber et al. 2019, Kakkad et al. 2020).

Several models have suggested that these massive and fast outflows are able to
suppress star formation (SF) activity in the host galaxy by removing and
heating the interstellar medium (ISM; e.g. King 2005, Costa et al. 2015).
However, from the observational point of view, it has not yet been clearly
assessed whether and how outflows affect the global level of SF in the host
galaxies (see e.g. Balmaverde et al, 2016, Woo et al. 2017).

Active galactic nucleus (AGN) feedback is expected to reach its maximum
efficiency at $z\sim 2$, where both SF and black hole accretion histories peak
(e.g. Madau & Dickinson 2014). Essentially all massive early-type galaxies are
already formed at $z\sim 1$ (Cimatti, Daddi & Renzini 2006), with indications
that quenching must have taken place several billion years before (Renzini
2006). The impact of AGN-driven outflows critically depends on the ISM
properties and gas content, which are different at $z\sim 2$ from what is
observed in local analogues (e.g. Kewley et al. 2013; Steidel et al. 2014;
Coil et al. 2015, Tacconi et al. 2020): It is thus fundamental to study AGN
outflows at these redshifts.

Figure 1:  H$\beta$ (left panel) and H$\alpha$ (right panel) spectral region
fit in $5\times 5$ NIRSpec spaxels centred on the QSO position. In the upper
panels, the data are shown in black, and the best-fitting broad line region
components (e.g. H$\beta$, H$\alpha$, and FeII) are shown in blue, the narrow
line region/host galaxy components (i.e. H$\beta$, [OIII], H$\alpha$, [NII],
and [SII]) are shown in green, and the total flux is presented in red. The
lower panels show the fit residuals.

XID2028 ($z$=1.5930), originally discovered in the XMM-COSMOS survey (RA
10:02:11.27 DEC +01:37:06.6; Brusa et al. 2010), is one of the best-studied
objects in this crucial feedback phase. From an optical point of view, XID2028
is classified as a type 1.8-1.9 quasar on the basis of the significant
extinction in the broad line region (BLR) component as measured from the
H$\alpha$/H$\beta$ ratio (see Brusa et al. 2015a), and the lack of broad
emission lines in the rest-frame UV. Based on the large multi-wavelength
coverage, it has reasonably well-constrained bolometric luminosity $\rm
log(L_{AGN})=46.3$, stellar mass $\textrm{M}_{*}\sim 4.5\times 10^{11}\
\textrm{M}_{\sun}$ , and SFR $\sim 250\ \textrm{M}_{\sun}/\textrm{yr}$, which
were derived from fits of the UV to far-IR (FIR) spectral energy distribution
(Lusso et al. 2012, Brusa et al. 2015a, Perna et al. 2015).

The presence of a massive ($\rm\dot{M}\sim 300\
\textrm{M}_{\odot}/\textrm{yr}$) and extended ($\sim 13$ kpc) ionised outflow,
traced by [OIII]$\lambda$5007 emission, has independently been advocated by
X-shooter slit spectroscopy (Brusa et al. 2015a, Perna et al. 2015) and
SINFONI111Spectrograph for INtegral Field Observations in the Near Infrared
IFU observations (Cresci et al. 2015a; C15). The outflow appeared to be
located at the centre of a cavity traced by rest-frame U-band HST/ACS imaging
and possibly by a narrow H$\alpha$ component attributed to star formation in
the host by C15 (but this scenario has been contested by Scholtz et al. 2021).
This suggested that the wind removes the gas from the host galaxy (so-called
negative feedback; see e.g. also Cano-Díaz et al. 2012, Carniani et al. 2016,
Cresci & Maiolino 2018), and simultaneously triggers SF by induced pressure at
the outflow edges or directly in the wind (so-called positive feedback; see
e.g. Silk 2013, Cresci et al. 2015b, Maiolino et al. 2017, Perna et al. 2020).

Perna et al. (2015) also reported the detection of sodium
NaID$\lambda\lambda$5890,5896, magnesium MgII$\lambda\lambda$2796,2803, and
MgI$\lambda$2853 absorption lines in the X-Shooter optical and near-IR (NIR)
spectra of the source. The absorption lines can be ascribed to an outflowing
neutral/low-ionisation gas component along the line of sight, with roughly the
same velocity shift as the [OIII] emission line.

XID2028 was then observed by Plateau de Bure Interferometer in CO(3-2) (Brusa
et al. 2015b), and by ALMA in CO(2-1), CO(5-4) and in the 1.3 mm continuum by
Brusa et al. (2018; B18). They derived a total molecular gas mass of $\sim
10^{10}\ \textrm{M}_{\odot}$, which corresponds to a very low gas fraction
($<$5%) and depletion times (i.e. the timescale during which all the gas
reservoir is converted into stars) of only $40-75$ Myr, supporting a scenario
in which outflows affect the gas reservoirs. While the bulk of the CO emission
appeared to be located in a rotating disk in the inner few kiloparsec, B18
also detected extended high-velocity gas, which they interpreted as a
signature of a galaxy-scale molecular outflow that is spatially coincident
with the blueshifted ionised gas outflow. A possible redshifted counterpart of
the outflow was also detected on the opposite side of the galaxy (see Fig. 8).
Summing up both the ionised and molecular gas components, B18 estimated a
total mass rate of the outflowing gas of $\rm\dot{M}_{\textrm{out,tot}}\sim
500-800\ \textrm{M}_{\odot}/\textrm{yr}$.

The available numerous multi-wavelength observations and its unique properties
therefore make XID2028 one of the best objects for characterising outflows and
feedback at z$\sim$1-2. It was thus selected to be observed with the Near
Infrared Spectrograph (NIRSpec) and Mid-Infrared Instrument (MIRI) integral
field units (IFU) on board the James Webb Space Telescope (JWST) during Early
Release Science (ERS) observations. In this paper, we present the NIRSpec
dataset, which reveals a wealth of structures that were hidden in the previous
seeing-limited ground-based data.

The paper is organised as follows: in Sect. 2 we present the observations,
data reduction, and analysis. The results are presented in Sect. 3, and we
draw our conclusions in Sect. 4. Throughout the paper, we adopt the
cosmological parameters $\textrm{H}_{0}=70\ \textrm{km}\ \textrm{s}^{-1}\
\textrm{Mpc}^{-1}$, $\Omega_{\textrm{m}}=0.3,$ and $\Omega_{\lambda}=0.7$
(Spergel et al. 2003).

##  2 NIRSpec observations, data reduction, and analysis

###  2.1 Observations

XID2028 was observed with NIRSpec (Jakobsen et al. 2022, Böker et al. 2022) on
board the JWST on 20 November 2022 as part of the Q-3D ERS Proposal 1335 (PI
D. Wylezalek). The program aims at studying the properties of three well-
studied and luminous quasars (QSOs) with the unprecedented sensitivity,
spatial resolution, and spectral coverage provided by the JWST (Wylezalek et
al. 2022).

The NIRSpec IFU observations provide spatially resolved imaging spectroscopy
in the wavelength range $0.97-1.82\ \mu m$ at $\textrm{R}\sim 2700$ with the
G140H/F100LP grating/filter pair over a $3"\times 3"$ field of view with
$0.1"\times 0.1"$ spatial elements. The FWHM of the NIRSpec PSF is $\sim 45$
mas at 1.3 $\mu m$, which is the observed wavelength of [OIII]$\lambda$5007
emission for this source. The observations were taken with an NRSIRS2 readout
pattern with 16 groups, using a nine-point medium-cycling dither pattern. This
resulted in a total exposure time of 2.95 hours.

Figure 2: Representative examples of our multi-component fitting of H$\beta$
and [OIII] emission lines in single spaxels from three different regions of
the galaxy, as marked with coloured crosses in Fig. 3. From the top, the
panels refer to a nuclear region spaxel, to a filament spaxel (see text), and
to a spaxel in the most blueshifted wind region, respectively. The NIRSpec
data are shown in black, the best-fitting systemic narrow components are shown
in cyan, and the Gaussian components used to reproduce the outflow are shown
with different tones of green. We classified the Gaussian components in each
spaxel with a velocity shift —v—¡300 km/s (with respect to the systemic
velocity of the QSO, marked with a dotted line) and $\sigma<300$ km/s as
narrow (see Sect. 3.2 and 3.6.1). Systemic narrow line emission dominates in
the nucleus and becomes fainter along the filaments, while only outflow
emission is detected in the wind region (see Sect. 3.1). Residuals of the
emission-line fitting are shown in blue in a sub-panel for all regions.
Figure 3:  Channel maps for the [OIII]$\lambda$5007 line emission. The six
panels show the [OIII] line flux in different velocity bins. In the upper left
panel, the highest blueshifted velocities are shown, highlighting the fastest
part of the wind. In the following panels, the filaments connecting the
outflow to the quasar location become evident. In contrast, the lower right
panel shows the possible redshifted outflow emission to the north-east. The
maps show the spaxel with S/N$>$3 on the total [OIII] flux. At the distance of
the target, the scale is $\sim 8.5$ kpc/”. The coloured crosses show the
location of the extraction spaxels for the spectra shown in Fig. 2.

###  2.2 Data reduction

The data reduction was performed with the JWST calibration pipeline version
1.8.2, using the context file jwst_1014.pmap. The individual raw images were
first processed for detector-level corrections using the Detector1Pipeline
module of the pipeline (stage1). Then, the resulting count-rate images were
calibrated using the Calwebb_spec2 module (stage2). At this stage, WCS-
correction, flat-fielding, and flux calibrations were applied to convert the
data from units of count rate into flux density. At the time of writing, many
NIRSpec calibration files contain placeholder values: for this reason, no
proper flat-fielding has been applied to the data. This has implications for
the shape of the background continuum emission, but does not affect the
analysis of the bright line emission from the QSO. The individual stage2
images were then resampled and coadded onto a final data cube through the
Calwebb_spec3 processing (stage3). We also applied a number of additional
steps (and custom modifications of the current pipeline code) to improve the
data reduction quality. Moreover, different configurations were also used to
obtain additional data products and test the pipeline robustness (e.g. of flux
and spatial resolution recovery). In particular:

  * •

We applied the snowball flagging for the jump step during the first stage of
the pipeline. This removes data artefacts caused by large cosmic-ray impacts
on the detector.

  * •

We further processed each individual count-rate frame to correct for varying
zero-level in the dithered frames at the end of stage1. The median value
(computed considering the entire image) was subtracted from each frame to
obtain a base level consistent with zero counts per second.

  * •

The obtained count-rate images were also processed to remove the $1/f$ noise
pattern, a low-frequency noise for which the noise power is inversely
proportional to the frequency. This appears as correlated vertical noise,
which we modelled in each column (i.e. along the spatial axis) with a low-
order polynomial function after flagging and removing all bright pixels (e.g.
associated with the observed target) through a sigma-clipping algorithm. The
modelled $1/f$ noise was then subtracted before proceeding with stage2 of the
pipeline.

  * •

The flux calibration was performed using the commissioning observations of the
standard star TYC 4433-1800-1 (PID 1128, o009). We reduced the star
observations with the same pipeline version and context file as the science
frames, and we used it to derive the response curve of the instrument required
to convert count rates into flux densities.

  * •

The outlier_detection step in stage3 was used to identify and remove the
numerous remaining spikes from the datacubes that are due to cosmic rays and
other artefacts and that were not corrected for in the previous reduction
steps. Unfortunately, this step cannot be used in the current version of the
pipeline, as it identifies too many false positives and seriously compromises
the data quality by removing too many spaxels with reliable signal. We
therefore used a post-processing correction based on a $\sigma$-clipping cut
to exclude all spikes in the reduced datacubes (at the spaxel level).

  * •

We applied the cube_build step to produce a final datacube with a spaxel size
of $0.05"$, obtained with the drizzle weighting. In this step, we patched the
cube_build script by fixing a bug affecting the drizzle algorithm, as
implemented in version 1.8.5222cube_build code changes in
https://github.com/spacetelescope/jwst/pull/7306.

  * •

As the final astrometric solution calculated by the pipeline resulted in an
offset of $\sim 0.1"$ of the QSO position, we corrected the WCS using the
archival HST/ACS imaging in the F814W filter as reference (Koekemoer et al.
2007, C15) to align the bright QSO nucleus in NIRSpec with its expected HST
position.

###  2.3 Data analysis

For the spectral analysis of the data, we adopted the fitting code presented
in Marasco et al. (2020, 2022) and Tozzi et al. (2021) to analyse IFU data of
local and high-z AGN and star-forming galaxies. The first step of the data
analysis models the unresolved BLR. We fit a BLR model to a nuclear spectrum
extracted in a $5\times 5$ spaxel ($0.25"\times 0.25"$) region around the AGN
location. Following C15, we fitted the asymmetric broad lines (H$\alpha$,
H$\beta$, H$\gamma$ and HeII$\lambda$4686) with a broken power-law
distribution convolved with a Gaussian. The broad line centroids, power-law
indices, and widths were tied to be the same for all lines, with the exception
of the H$\alpha$ Gaussian width, which required a slightly different width to
reproduce the observed spectrum. The FeII emission was reproduced using a
combination of different CLOUDY (Ferland et al. 1998) models. Two additional
BLR Gaussian components were added at $\sim-8000$ and $\sim-16000$ km
$\textrm{s}^{-1}$ from H$\alpha$ (see e.g. Carniani et al. 2016). The
forbidden lines and the narrower components of the BLR lines were fitted with
two Gaussian components each, forcing each of them to have the same velocity
and width for all species. The resulting fit is shown in Fig. 1 for the
spectral regions around H$\beta$ and H$\alpha$.

We used the derived BLR model to subtract its unresolved emission from the
data spaxel by spaxel. We used the pPXF code (Cappellari 2017) to fit the
entire datacube, allowing the BLR template to vary only in flux, as expected
for unresolved emission that is distributed as the instrumental point spread
function (PSF). We additionally modelled the residual continuum with a first-
order polynomial and the extended line emission of the main species as
combinations of multiple Gaussian components, ranging from one to three. A
Kolmogorov-Smirnov test on the residuals was used to select the optimal number
of components in each spaxel to reproduce the observed profile while
minimising the number of free parameters (see Marasco et al. 2020). The
rescaled best-fitting BLR template and continuum models were then subtracted
spaxel by spaxel, obtaining a BLR-subtracted cube containing only the extended
line emission due to the NLR, the host galaxy, and the outflow components.

We then refit the BLR-subtracted cube to carefully model the residual line
emission. We fixed the ratio of the lines in each of the
[OIII]$\lambda$4958,5007 and [NII]$\lambda$6548,6584 doublets to the
theoretical value of 3, and we forced each Gaussian component of all emission
lines to have the same line profile shape (centroid velocity and velocity
dispersion), assuming that they are emitted from regions with similar
kinematics. This reduced the degeneracy in the fitting procedure, especially
for the H$\alpha$ and [NII] multiplet in spaxels with complex kinematics. As
before, we performed the fit with up to three Gaussian components for each
line, and we selected the optimal number of components to reproduce the line
profile spaxel by spaxel. Three examples of multi-component line fitting in
different regions of the galaxy are shown in Fig. 2. Their locations are
marked with coloured crosses in Fig. 3.

Figure 4:  Kinematics of the [OIII] line emission. The upper panels show
moment 0, moment 1, and moment 2 of the whole line profile velocity. The lower
panels instead show $\textrm{v}_{10}$, the velocity at the 10th percentile of
the overall emission-line profile in each spaxel, $\textrm{v}_{90}$, the 90th
percentile, and $\textrm{W}_{40}$, i.e.,
$\textrm{W}_{40}=\textrm{v}_{50}-\textrm{v}_{10}$, the line width containing
40% of the emission line flux. Solid contours represent 30% and 50% of the
peak emission in the bluest channel map ($\textrm{v}<-1000$ km/s) in Fig. 3,
and the dashed contours the 5% and 10% emission of the third channel
($-600<\textrm{v}<-250$ km/s). The maps show the spaxel with S/N$>$3 on the
total [OIII] flux, and the cross marks the position of the QSO.

##  3 Results

###  3.1 Bubbles and outflows: Extended ionised gas kinematics

Fig. 3 shows channel maps of the modelled [OIII] line flux in different
velocity bins. The same pattern is also present in H$\alpha$, H$\beta,$ and
the other emission lines, but we present a kinematic analysis of [OIII] here
because this bright forbidden line is not contaminated by emission from the
BLR and is thus an ideal tracer of outflows. Prominent collimated blue
emission is evident in the most blueshifted velocity bin ($\textrm{v}<-1000$
km/s). This highly blueshifted emission is located in the same region as the
[OIII] outflow detected in earlier SINFONI data by C15, as well as the
blueshifted molecular outflow discovered by B18. Notably, as seen also in the
SINFONI data, the blue [OIII] outflow emission peaks at a projected distance
of $\sim 0.7"$ ($\sim 6$ kpc) west of the QSO.

At lower blueshifted velocity ($-600<\textrm{v}<-250$ km/s), a cavity of
suppressed line emission becomes evident between the QSO and the outflow,
surrounded by bright [OIII] filaments that seem to connect the QSO with the
region with the fastest outflow emission. This filamentary structure around
the cavity roughly corresponds to the elongated regions detected in HST/ACS
imaging and narrow H$\alpha$ in C15.

The kinematics maps in Fig. 4 also highlight the features discussed above. The
figure shows the moment 0, moment 1, and moment 2 maps of the total
[OIII]$\lambda$5007 emission line velocity in the upper panels, with
superimposed contours of the moment 0 maps. The bottom panels show
$\textrm{v}_{10}$, the velocity at the 10th percentile of the overall
emission-line profile in each spaxel, $\textrm{v}_{90}$, the 90th percentile,
and $\textrm{W}_{40}$, defined as the line width containing 40% of the
emission line flux $\textrm{W}_{40}=\textrm{v}_{50}-\textrm{v}_{10}$. In this
parametrisation, v10 traces the highest-velocity blueshifted gas in the data,
with velocities as high as $\textrm{v}_{10}\sim-1000$ km/s towards our line of
sight. These velocities are not compatible with either rotational motions in
the host, given their high values and their direction perpendicular to the
velocity gradient in the host (see B18), or with star formation-driven
outflows, which are generally characterised by lower peak velocities (see e.g.
Arribas et al. 2014, Fluetsch et al. 2019, Förster Schreiber et al. 2019). We
also exclude the possibility that they may be caused by a companion galaxy
because the line width at the highest blueshifted velocities is
$\sigma_{\textrm{blue}}\sim 450$ km/s, which would imply a very high-mass
companion galaxy that is undetected in continuum333A continuum feature,
located nearby ($\sim 0.3"$) but not coincident with the blueshifted emission,
is in fact detected only in one of the two detectors covering the spectra at
that location, at wavelengths between 1.465 $\mu m$ and 1.89 $\mu m$, while it
suddenly disappears at wavelengths $\lambda<1.43\ \mu m$ in the other
detector. Because the data were not flat-field corrected for the differential
pixel-to-pixel variations because a reliable flat-field calibration is not yet
available for NIRSpec, and because multiple similar continuum structures were
present in the datacube with different spatial positions in wavelength regions
covered by different detectors, we attribute this feature to pixel-to-pixel
continuum level variations and consider the interpretation as a companion
galaxy very unlikely., nor is it dected in deep ACS imaging of the source (see
C15).

Figure 5: Moment 1 and moment 2 maps for the component of the
[OIII]$\lambda$5007 line emission tracing gas in the host galaxy (see text).
The contours are the same as in Fig. 4. An S/N$>$40 threshold has been used to
mask lower S/N spaxels and isolate the rotation in the central spaxels. A
rotational pattern is suggested by the data that is compatible with the
CO(5-4) rotation detected by B18.

Similarly, $\textrm{v}_{90}$ traces the highest redshifted velocities, up to
v${}_{90}\sim 800$ km/s. These are located towards the north-east, opposite to
the QSO location with respect to the blueshifted velocities, suggesting that
the gas is outflowing in opposite directions from the central source, as
predicted in bi-conical outflow models. This emission spatially corresponds to
the redshifted molecular outflow discussed in B18 (see Fig. 4 and Fig. 8,
upper panel), and thus might represent the receding ionised outflow that was
undetected in the earlier SINFONI data. The location of the redshifted
outflows also shows the highest line widths, with W${}_{40}\sim 700$ km/s.

The moment 0 [OIII] map also shows the filamentary structure that was
identified in the channel maps, which spatially and in velocity connects the
QSO to the collimated blueshifted outflow. As shown in the moment 1 map, the
filaments show a symmetric intermediate velocity of $\sim-300$ km/s between
the QSO and the outflow (see Fig. 2), and their velocity gradient shows
increasing blueshifted velocities towards the position of the collimated
outflow. This suggests that the two are physically connected in some way.

Highly blueshifted emission, but with fainter [OIII] emission, is also
detected towards the northwest, with v10 velocities comparable to those
observed in the outflow. The detection of a possible continuum source in the
NIRSpec datacube at that location may suggest a companion galaxy. Reliable
flat-field calibration for NIRSpec is required to test this hypothesis, but is
not yet available, however.

###  3.2 Possible rotating disk in the host galaxy

We also attempted to measure the gas kinematics in the disk of the host
galaxy. To do this, we selected the Gaussian component with the highest peak
in each spaxel, excluding those with a velocity shift $|\textrm{v}|>300$ km/s
or with $\sigma>300$ km/s. We show the resulting moment 1 and moment 2 maps in
Fig. 5, where we impose a cut of S/N$>$40 in the displayed spaxels to better
isolate the rotation-like pattern in the central region. Interestingly, the
velocity gradient obtained in this way is consistent with the one revealed by
B18 in CO(5-4), which was interpreted as a compact rotating disk in the host
galaxy. In particular, both the position angle PA $\sim 125^{\circ}$ and the
magnitude of the velocity gradient ($-200\lesssim\textrm{v}\lesssim 200$ km/s)
are consistent with the molecular gas measurements. Assuming that this
velocity gradient traces the major axis of the host galaxy, this orientates
the red- and blueshifted outflow components roughly along the minor axis.
Assuming an inclination of 30∘, B18 derived a dynamical mass
$\textrm{M}_{dyn}\sim 6-8\times 10^{10}\ \textrm{M}_{\odot}$ within the region
in which CO is detected (r$\sim 0.25"$), using both the Cresci et al. (2009)
and the Di Teodoro & Fraternali (2015) kinematical fitting codes.

Figure 6: Electron density map from [SII]$\lambda\lambda$6716,30 line ratio
(left panel) and extinction map from the Balmer decrement H$\alpha$/H$\beta$
(right panel). The contours used in Fig. 4 are overplotted. Both maps display
only the spaxels with at least S/N$>$3 in all the emission lines involved, and
the QSO location is marked with a cross.  Figure 7: Resolved [SII]-BPT (upper
panels) and [NII]-BPT (lower panels) diagrams for each spaxel with S/N$>$3 in
each line. The data points in the BPT diagrams (left panels) are colour coded
as a function of their [SII]/H$\alpha$ and [NII]/H$\alpha$ ratios for the
[SII]-BPT and [NII]-BPT, respectively. The same colours are used in the
corresponding maps (right panels), where the contours shown in Fig. 4 are
overplotted. The QSO location is marked with a black cross.

###  3.3 ISM properties

We explored the physical conditions of the ISM using the observed line ratios
in the NIRSpec datacube. The electron density ne was derived using the
[SII]$\lambda\lambda$6716,30 line ratio. The [SII]$\lambda\lambda$6716,30
lines were fitted using the same components used for the brighter forbidden
lines, although the deblending into different components is very uncertain due
to the low S/N of the broader components. We therefore used the total
[SII]$\lambda\lambda$6716,30 line flux, the parametrisation of Sanders et al.
(2016) and assumed Te=$10^{4}$ K. In the left panel of Fig. 6, we show the
resulting $\rm n_{e}$ map, which is noisy because the [SII] line emission is
faint and the modelling of the two distinct components of the doublet may have
degeneracy issues. However, the density seems uniform across the sampled
regions where the [SII] lines are detected with S/N$>$3, with values between
$200-500$ cm-3. This is compatible with similar measurements at $z\sim 2$ in
star-forming and active galaxies (see e.g. Sanders et al. 2016, Masters et al.
2014, Steidel et al. 2014, Kaasinen et al. 2017, Kashino et al. 2017).
Although a higher density has been found in the outflow component than in the
disk component in the hosts of local and high-z AGNs (see e.g. Mingozzi et al.
2019, Förster Schreiber et al. 2019), we derive similar values in the central
region dominated by the QSO and on the outflow region within the errors. In
the outflow region, where the entire emission is blueshifted and samples the
outflow, we measured a median of n${}_{e}\sim 360\pm 180$ cm-3, which is
consistent with the value derived by an integrated spectrum centred on the
outflow region with a radius of 3 spaxels.

We used the Balmer decrement, H$\alpha$/H$\beta$, to derive the dust
attenuation map, assuming a Calzetti et al. (2000) attenuation law and a
temperature of $10^{4}$ K. The map we obtained is shown in the right panel of
Fig. 6. The optical extinction seems to decrease towards the nucleus despite
the red colours of the QSO. The reason probably is that the rest frame optical
lines cannot probe deep in the highly embedded regions around the QSO (see
e.g. Cresci et al. 2010, Cresci et al. 2017). MIRI observations will be useful
to probe the dust that is located deeper in this region. Some unsubtracted
residuals of the BLR in H$\beta$ might also contribute. The extinction also
decreases along the filaments, where shocks might contribute to destroying the
dust. A ratio H$\alpha$/H$\beta$$\sim 4$ is measured in the highly blueshifted
outflow region, corresponding to a median extinction
$\textrm{A}_{\textrm{V}}=1.1\pm 0.5$, where the error is given by the standard
deviation of $\textrm{A}_{\textrm{V}}$ over the outflow region.

We also used the so-called BPT diagrams (Baldwin et al. 1981, Kewley et al.
2006, Lamareille et al. 2010, Law et al. 2021) to investigate the dominant
ionisation source in the different regions of the galaxy. As shown in Fig. 7,
all spaxels are AGN dominated according to [SII]-BPT
([SII]$\lambda\lambda$6716,30/H$\alpha$ vs [OIII]$\lambda$5007/H$\beta$) and
[NII]-BPT ([NII]$\lambda$6584/H$\alpha$ vs [OIII]$\lambda$5007/H$\beta$)
relative to the full line emission. We find that the low-ionisation lines
[SII] and [NII] are enhanced along the outflow and in a cavity surrounded by
the filaments (see Sect. 3.4), possibly by shocks. The BPT diagrams shown in
Fig. 7 are relative to the full line profile after removing the BLR component.

By inspecting fit residuals, we searched for a residual narrow H$\alpha$
component tracing star formation in the host galaxy. In previous SINFONI data,
C15 detected such an SF component in two distinct spots around the outflow
location, where NIRSpec has now unveiled gas filaments. The component was
later not detected in the different analysis performed by Scholtz et al.
(2021). As the residuals of the H$\alpha$ fit are critically dependent on the
assumptions made to perform the fit and the BLR subtraction, we fixed the
parameters of the fitting of the BLR-subtracted H$\alpha$ to the [OIII] best-
fitting parameters to minimise the degeneracy. Most of the (small) residuals
in the H$\alpha$ flux after subtracting the best fit are detected in the wings
of the line, and we do not identify them as additional kinematic components,
but rather with an imperfect parametrisation of the BLR wings. Nevertheless,
the elongated filaments observed in the rest-frame U band HST/ACS imaging
still suggest that the presence of an outflow and an expanding bubble of hot
gas, as discussed in the next section, might have a significant effect in
triggering star formation and compressing the ISM at their edges. The related
H$\alpha$ emission might be difficult to disentangle from the bright QSO
contribution, however.

Figure 8: Channel map of [OIII]$\lambda$5007 emission with v$<-300$ km/s with
contours of the CO(5-4) molecular blue- (cyan) and redshifted (magenta)
outflow components (upper panel, from B18), and 3 GHz radio emission (white)
from VLA observations (lower panel, from Vardoulaki et al. 2019). The location
of the QSO is marked with a black cross, and the filaments and the cavity are
shown with arrows. The VLA radio emission and the blueshifted ionised and
molecular outflows appear to be co-spatial, as well as the possible redshifted
counterparts.

###  3.4 Jet-ISM interaction in XID2028

Vardoulaki et al. (2019) presented radio images of XID2028 among the multi-
components radio sources identified in the VLA-COSMOS Large Project at 3 GHz
(0.75” resolution, 2.3 $\mu\textrm{Jy}\ \textrm{beam}^{-1}$ RMS). The galaxy
is identified as number 10964 in their catalogue, but is surprisingly
classified as a star-forming galaxy based on the lack of a radio excess with
respect to the IR-radio correlation of star-forming galaxies (Delvecchio et
al. 2017), as well as on the lack of clear jets at 3 GHz and of any compact
nuclear radio emission with VLBA observations ($10\ \mu\textrm{Jy}\
\textrm{beam}^{-1}$, at $16.2\times 7.3\ \textrm{mas}^{2}$, Herrera Ruiz et
al. 2017). In this respect, we note that subsequent deeper VLBA+GBT
observations in different parts of the COSMOS field have produced 40% more new
AGN detections (Herrera Ruiz et al. 2018), suggesting a possible sensitivity
limit of current observations.

Figure 9: Comparison between the $\mathrm{MOKA^{3D}}$ weighted model (upper
panels) and the NIRSpec [OIII] kinematics (bottom panels) for moment 0, moment
1, and moment 2 in the higher S/N region of the galaxy, showing how the
assumed bubble+jet geometry is able to reproduce the observed kinematics. The
cross marks the location of the QSO.

However, in the 3 GHz VLA image, the emission peak is not coincident with the
QSO and host galaxy centre, as expected if the radio emission were produced by
star formation alone (see Fig. 8, bottom panel). Instead, we observe two blobs
that extend north-east and south-west, in opposite directions from the centre,
at locations roughly corresponding to the red- and blueshifted outflows
detected in the ionised and molecular data (see Fig. 8, upper panel). The two
blobs can thus be interpreted as radio lobes that extend beyond the host
galaxy, produced by a bi-polar jet that is co-spatial with the ionised and
molecular wind (see e.g. Cresci et al. 2015b, Jarvis et al. 2019, Venturi et
al. 2021). We also note that based on the integrated source flux density at
1.4 GHz (162 $\mu$Jy, Schinnerer et al. 2007) and at 3 GHz (92 $\mu$Jy,
Smolčić et al. 2017, Vardoulaki et al. 2019), the radio spectral index is
$\alpha=0.74$, consistent with synchrotron emission from AGN radio lobes. The
corresponding radio luminosities, 6.8$\times 10^{23}$ W/Hz and 3.9$\times
10^{23}$ W/Hz at 1.4 and 3 GHz, respectively, would indeed also be compatible
with a low-power radio galaxy.

The coexistence of a bi-polar radio jet with ionised and molecular gas
outflows has been observed in several sources (see e.g. Morganti et al. 2015,
Cresci et al. 2015b, Venturi et al. 2021, Ulivi et al. in prep., Giroletti et
al. 2017 and reference therein). In this scenario, the filamentary structure
detected with NIRSpec between the QSO and the outflow can be interpreted as a
hot expanding bubble filled with low surface brightness emission that is
brightened at the edges, which is inflated and dragged into the host galaxy
ISM by the jet. According to the Ks LUCI+ARGOS imaging presented in B18, this
is compatible with the host galaxy effective radius $\rm r_{e}=10$ kpc, as the
bubble extends to distances $\rm D_{bubble}=0.7"\sim 6$ kpc. Similar
structures are also observed in low radio luminosity systems at intermediate
redshift (e.g. the Teacup galaxy; Harrison et al. 2015, Venturi et al. in
prep.).

Ionised gas filaments surrounding radio lobes have also been observed in
higher-luminosity radio galaxies in the nearby Universe, where complex
networks of line-emitting streams are observed around the hot X-ray cavities
(e.g. 3C317; Balmaverde et al. 2018). However, in these bright radio galaxies,
the filaments show complex kinematics due to the combination of emission from
receding and approaching streams. In contrast, the filaments in XID2028 appear
to have consistent kinematics, compatible with tracing the brightened edge of
an expanding bubble. According to simulations (e.g. Gaibler et al. 2012), the
expanding bubble and the jet excavate the central region of the galaxy disk,
creating a cavity and possibly triggering star formation at the edges of the
expanding shell by a blast wave causing strong compression and cooling in the
ISM (e.g. Silk et al. 2013, Ishibashi & Fabian 2012, C15, Cresci et al. 2015b,
Cresci & Maiolino 2018, Duggal et al. 2021, Bessiere & Ramos Almeida 2022).

At a larger radius ($\sim 0.7"$, corresponding to 6 kpc), the jet appears to
have pierced the shell, and it propagates faster through the lower-density
environment together with the wind. This creates the collimated and highly
blueshifted outflows already detected with SINFONI by C15, with approaching
projected velocities as high as 1100 km/s.

A similar scenario is reproduced in 3D hydrodynamic simulations of jets
interacting with a clumpy ISM (see e.g. Gaibler et al. 2011, Gaibler et al.
2012, Wagner et al. 2012, Dugan et al. 2017, Mukherjee et al. 2018, Huško et
al. 2022). In these simulations, the jet initially inflates energy-driven
bubbles in the host galaxy ISM that expand perpendicularly to the galaxy disk.
The jet eventually pierces the bubble and propagates into the halo. The
asymmetry observed between the approaching and receding side of the outflow
can also be ascribed to differences in the ISM density and clumpiness between
the two sides of the bipolar jet (Gaibler et al. 2011).

Figure 10: Best-fitting $\mathrm{MOKA^{3D}}$ model used to reproduce the
XID2028 bubble+jet kinematics in Fig. 9. Upper panel: Projected unweighted
moment 1 map. Middle panel: 3D unweighted model, i.e. the intrinsic cloud
distribution before weighting them on the observed datacube. A schematic 3D
view of the assumed toy-model geometry and kinematics is superposed. The QSO
position is shown with a star, and the arrows show the outflow direction
together with the expansion velocity that inflates the bubble (see text).
Lower panel: 3D model. Each cloud is weighted based on the observed flux for
its location and velocity in the data.

###  3.5 Modelling the observed gas kinematics with $\mathrm{MOKA^{3D}}$

To test whether the shell and piercing jet geometry predicted by simulations
can be applied to XID2028, we compared the NIRSpec XID2028 data to the 3D AGN
outflow model generated using the modelling framework called modelling
outflows kinematics in AGN 3D ($\mathrm{MOKA^{3D}}$) that is presented in
Marconcini et al. (in prep.). The model assumes a simple velocity field
described by an analytical expression, but populates it with a 3D distribution
of emitting clouds weighted using the observed line flux in all spaxels and
spectral channels of the data cube. This procedure improves upon the
assumption of a smooth distribution of gas-emitting regions used in the
literature because it is able to reproduce the clumpy appearance observed in
the data. Details of the model implementation are discussed in Marconcini et
al. (in prep.). Here we just show that the observed gas geometry can be
reproduced using an expanding bubble dragged by the jet and the wind, plus a
collimated outflow where the radio jet escapes the host galaxy ISM. A
comparison between the $\mathrm{MOKA^{3D}}$ model with observations is
presented in Fig. 9, highlighting that this configuration can simultaneously
reproduce the filaments, the cavity and the collimated and fast outflow beyond
the bubble in flux and velocity. The best-fitting unweighted 3D model, that
is, the starting model with the intrinsic cloud distribution before weighting
the clouds on the observed datacube, is shown in Fig. 10 (middle panel), along
with its 2D velocity distribution (upper panel). The weighted model, in which
each cloud has been weighted based on the observed flux for its location and
velocity in the data, is shown in the bottom panel.

The estimate of the intrinsic deprojected outflow velocity in the best-fitting
model is $\textrm{v}_{out,\mathrm{MOKA^{3D}}}=1027$ km/s, with an aperture
angle $\alpha=10^{\circ}$, while the bubble has a radial velocity of
$\textrm{v}_{rad,\mathrm{MOKA^{3D}}}=509$ km/s, an expansion, a tangential
velocity $\textrm{v}_{exp,\mathrm{MOKA^{3D}}}=517$ km/s, and an intrinsic
dispersion velocity $\textrm{v}_{disp,\mathrm{MOKA^{3D}}}=117$ km/s. The
bubble and outflow system is reproduced with an inclination of $\sim
27^{\circ}$ with respect to the observer’s line of sight.

###  3.6 Outflow energetics and mass outflow rate

####  3.6.1 Ionised outflow

Through the new NIRSpec data with their higher resolution and sensitivity, we
were also able to revise the estimate of mass outflow rate based on SINFONI
data by C15. We can directly measure the electron density in the outflow from
the [SII]$\lambda\lambda$6716,30 line ratio and used the H$\alpha$ flux
associated with the outflow to derive the line-emitting ionised gas mass.
Following Cresci et al. (2017), for instance,

$\textrm{M}_{out}=3.2\cdot
10^{5}\left(\frac{\textrm{L}_{out}(\textrm{H}\alpha)}{10^{40}\
\textrm{erg}/\textrm{s}}\right)\ \left(\frac{100\
\textrm{cm}^{-3}}{\textrm{n}_{e,out}}\right)\ \textrm{M}_{\sun},$ (1)

where $\textrm{L}_{out}(\textrm{H}\alpha)$ is the H$\alpha$ luminosity of the
outflow component, and $\textrm{n}_{e}$ is the electron density measured in
the outflow. To isolate the outflow emission, we excluded the Gaussian
components with a velocity shift —v—¡300 km/s (with respect to the systemic
velocity of the quasar) and $\sigma<300$ km/s (see Fig. 2) from each spaxel,
exhibiting disk-like kinematics as shown in Fig. 3.2. In this way, we derived
a total outflow H$\alpha$ flux $\textrm{F}_{out}(\textrm{H}\alpha)\sim 1.28\pm
0.04\times 10^{-16}$ erg s-1 $\textrm{cm}^{-2}$. We corrected this flux for
extinction using the H$\alpha$/H$\beta$ line ratio for the broad component, as
derived over the outflow region ($\textrm{A}_{\textrm{V}}=1.1\pm 0.45$, see
Sect. 3.3), obtaining $\textrm{F}_{out,corr}(\textrm{H}\alpha)=3\pm 1\times
10^{-16}$ erg s-1 $\textrm{cm}^{-2}$. Using the electron density for the
outflow component measured on the outflow region, $\textrm{n}_{e,out}=360\pm
180\ \textrm{cm}^{-3}$, we derived a total outflow mass
$\textrm{M}_{out}=4.5\pm 3\times 10^{7}\ \textrm{M}_{\sun}$. Assuming an
outflow velocity $\textrm{v}_{out}=-1100$ km/s, which is the observed maximum
$\textrm{v}_{10}$, consistent with the $\mathrm{MOKA^{3D}}$ intrinsic value,
and a radius $\textrm{R}_{out}=1.6"$ ($\sim$ 13.6 kpc), deprojected using the
inclination of $\sim 27^{\circ}$ from $\mathrm{MOKA^{3D}}$ to 3.5”, we obtain
an outflow rate (e.g. B18)

$\rm\dot{M}_{out,ion}=\textrm{v}_{out}\
\frac{\textrm{M}_{out}}{\textrm{R}_{out}}=6\pm 3\
\textrm{M}_{\sun}/\textrm{yr}.$ (2)

The computed value is consistent with the estimate from the
$\mathrm{MOKA^{3D}}$ model, from which we derive a median outflow rate in the
outflow volume of $\rm\dot{M}_{out,ion}(\textrm{$\mathrm{MOKA^{3D}}$})\sim 8\
\textrm{M}_{\sun}/\textrm{yr}$.

The main difference with respect to the analysis presented in C15
($\rm\dot{M}_{out,ion}\sim 100$) follows from our measured value of ne, which
they assumed to be n${}_{e}=100$ cm-3, thus pushing the outflow rate towards
higher values, by the deprojection of the outflow radius based on the
$\mathrm{MOKA^{3D}}$ modelling, and by the better definition of the outflow
component allowed by the NIRSpec data. The corresponding kinetic power of the
ionised outflow is given by

$\rm\dot{E}_{kin,ion}=\frac{1}{2}\rm\dot{M}_{out,ion}\textrm{v}^{2}=2.1\pm
1.4\times 10^{42}\ \textrm{erg}\ \textrm{s}^{-1},$ (3)

and the momentum rate $\rm\dot{P}_{kin,ion}=\dot{M}_{out,ion}\
\textrm{v}_{out}=4\pm 2\times 10^{34}$ dyne.

An estimate of the ionised outflow rate can be also derived from the observed
[OIII] flux in the broad component (e.g. Cano-Diaz et al. 2012). We measure a
dust corrected total [OIII] emission from the outflow
$\textrm{F}_{out}(\textrm{[OIII]})\sim 10\pm 5\times 10^{-16}\ \textrm{erg}\
\textrm{s}^{-1}\ \textrm{cm}^{-2}$, which translates into a total outflowing
mass $\rm M_{out}([OIII])=1.2\pm 0.8\times 10^{7}\ \textrm{M}_{\sun}$ assuming
solar metallicity and the same values for the density as were used for
H$\alpha$. This translates into $\rm\dot{M}_{out,ion}([OIII])=3\pm 2\
\textrm{M}_{\sun}/\textrm{yr}$, consistent with the H$\alpha$ value, but
slightly lower. However, given the more uncertain set of assumptions needed to
compute the ionised gas mass from [OIII] (e.g. gas metallicity and ionisation
structure; see Carniani et al. 2015), we chose to use the H$\alpha$ based
measurement as the fiducial mass.

####  3.6.2 Neutral outflow

The neutral sodium absorption lines NaID$\lambda\lambda$5890,5896 are clearly
detected in the NIRSpec nuclear spectrum. We fitted the observed intensity
profiles of the sodium lines in an integrated spectrum on a region centred on
the QSO with a radius of 10 spaxels (0.5”), where the absorption is detected,
considering a model parametrised in the optical depth space (e.g. Rupke et al.
2002). As for the emission lines, a multiple-component fit was required to
reproduce the observed profiles in absorption. The best-fit model reported in
Fig. 11 shows three distinct blueshifted kinematic components in absorption,
at $\rm v_{0,1}=-215$ km/s (FWHM=93 km/s, $\rm
log(N_{H,1}/\textrm{cm}^{-2})=20.68$), $\rm v_{0,2}=-260$ km/s (FWHM=234 km/s,
$\rm log(N_{H,2}/\textrm{cm}^{-2})=20.89$), and $\rm v_{0,3}=-640$ km/s
(FWHM=94 km/s, $\rm log(N_{H,3}/\textrm{cm}^{-2})=20.68$); no emission line
components are required for either the NaID or the HeI$\lambda$5877 line to
reproduce the observed profile in the fit. Following the prescriptions
indicated by Perna et al. (2015), we derived a new estimate for the neutral
outflow energetics, assuming that the neutral gas flow extends up to 1 kpc,

$\rm\dot{M}_{out,neut}=7\cdot\sum_{i=1}^{3}\left(\frac{\textrm{N}_{\textrm{H},i}}{10^{20}\
\textrm{cm}^{-2}}\right)\left(\frac{\rm R}{5\
\textrm{kpc}}\right)\left(\frac{\textrm{v}_{0,i}}{300\
\textrm{km}/\textrm{s}}\right)\sim 30\ \textrm{M}_{\odot}/\textrm{yr},$ (4)

where we sum the contributions of the three different velocity components.
This measurement should be considered as an order-of-magnitude estimate
because we cannot constrain the spatial extension of the neutral outflow on
the line of sight; nevertheless, this mass rate is broadly consistent with the
lower limit reported in Perna et al. (2015), $\rm\dot{M}_{out,neut}\gtrsim 80\
\textrm{M}_{\sun}/\textrm{yr}$, based on MgII absorption in their X-Shooter
spectra after accounting for the different assumed outflow radius (in Perna et
al. 2015 the same extension as the ionised outflow, 11 kpc).

Figure 11: Portion of NIRSpec integrated spectra of XID2028 including the NaID
doublet. The best-fitting components are superposed. The feature has been
reproduced with three blueshifted velocity components (see text). The vertical
lines indicate the systemic wavelengths of the HeI$\lambda$5877 and
NaID$\lambda$5890,5896 doublet lines.

####  3.6.3 Total outflow rate

The molecular outflow phase has been studied in detail by B18. Rescaling their
outflow rate taking the increased deprojected radius based on the
$\mathrm{MOKA^{3D}}$ inclination into account, we derived a molecular outflow
rate $\rm\dot{M}_{out,mol}\sim 20-120\ \textrm{M}_{\sun}/\textrm{yr}$
depending on the assumed CO-H2 conversion factor $\alpha_{\rm CO}$ (see B18
for details). Summing up all the three gas phases, we derive a total mass
outflow rate

$\rm\dot{M}_{out,tot}\sim 60-160\ \textrm{M}_{\sun}/\textrm{yr}.$ (5)

Based on this revised value for the total outflow rate and using the molecular
gas mass derived by B18 ($\textrm{M}_{gas}=1.1\pm 0.5\times 10^{10}\
\textrm{M}_{\sun}$) as well as the SFR$\sim 250\
\textrm{M}_{\sun}/\textrm{yr}$ derived from the total IR luminosity (Perna et
al. 2015), we can estimate the total depletion time for this source. Using the
numbers above, we derive
$\rm\tau_{depl}=\textrm{M}_{gas}/(\textrm{SFR}+\dot{M}_{out,tot})\sim 30$ Myr.
This timescale is remarkably shorter than typical gas depletion times in
normal star-forming galaxies with comparable stellar mass and redshift ($\sim
0.5-1.5$ Gyr; see e.g. Sargent et al. 2014, Tacconi et al. 2018). This further
supports the significant effect of AGN feedback in this source.

##  4 Conclusions

We have presented JWST NIRSpec IFU observations of XID2028, a prototypical
outflowing obscured QSO at z=1.59. A massive and extended outflow has
previously been detected in this source in different gas phases in previous
SINFONI (ionised), X-Shooter (neutral atomic), and ALMA (neutral molecular)
observations. The unprecedented sensitivity and resolution of NIRSpec allowed
us to reveal new details of the morphology, kinematics, and wind/jet
interaction with the ISM in this iconic source. In particular:

  * •

we detected a rotational pattern in the narrowest component of the emission
lines, which we interpreted as tracing a rotating disk in the host galaxy (see
Fig. 5). The velocity gradient is consistent with the inner kinematics
observed in CO(5-4) line emission observed with ALMA (B18), supporting the
estimate of a dynamical mass $\textrm{M}_{\textrm{dyn}}\sim 10^{11}$ M☉ within
the observed CO and ionised gas rotating structure.

  * •

A highly blueshifted and collimated outflow, with velocities up to $\sim-1100$
km/s, is detected $\sim 6$ kpc apart from the QSO, extending up to $\sim 13$
kpc in projection (see the upper left panel of Fig. 3). The ionised wind is
co-spatial with the molecular outflow detected in B18 and consistent with the
detection in [OIII] by C15. A possible redshifted outflow component is also
detected in the ionised phase. It is spatially consistent with the molecular
redshifted outflow counterpart in the opposite direction with respect to the
QSO and the blueshifted outflow, with velocities up to $\sim 700$ km/s (see
the lower right panel of Fig. 3 and the upper panel of Fig. 8).

  * •

The outflow appears to be connected to the QSO by a system of two line-
emitting filaments, surrounding a region with suppressed emission (see the
moment 0 map in Fig. 3.1), resembling shell-like features that have been
observed in local low-luminosity radio galaxies. The filaments have low
projected velocities close to the QSO, and show progressively more blueshifted
velocity in regions closer to the fast collimated outflow (see the moment 1
map in Fig. 3.1).

  * •

Archival VLA 3 GHz observations show two extended radio lobes that are broadly
coincident with the bi-polar outflows, suggesting low-luminosity radio jets in
the galaxy (see Fig. 8). We thus interpret the filamentary structure between
the QSO and the outflow as a hot expanding bubble filled with low surface
brightness emission brightened at the edges that is inflated and dragged into
the host galaxy ISM by the jet. At a larger distance from the QSO ($\sim 6$
kpc projected), the jet pierces the shell, and it propagates faster through
the lower-density environment.

  * •

We compared the NIRSpec data with the expectation from our
$\mathrm{MOKA^{3D}}$ (Marconcini et al. in prep.) AGN outflow 3D model
assuming this expanding bubble + collimated outflow geometry. We showed that
this scenario can reproduce the observed kinematics (see Fig. 9).

  * •

The proposed scenario is compatible with a contribution of both negative
feedback along the path of the outflow, where the gas is heated and removed
from the host galaxy, and positive feedback at the bubble edges, where the
compression of the ISM during the expansion can trigger star formation within
the feedback-driven shell, as traced by rest-frame U-band emission. However,
we do not detect a residual H$\alpha$ component ascribable to star formation
in the host.

  * •

We recomputed the ionised outflow energetics and mass outflow rate by directly
measuring the extinction and electron density in the outflow (presented in
Fig. 6) and based on the more detailed kinematics and structural information
allowed by the new NIRSpec data. We derived an ionised gas outflow rate
$\rm\dot{M}_{out,ion}=6\pm 3\ \textrm{M}_{\sun}/\textrm{yr}$ and a kinetic
power $\rm\dot{E}_{kin,ion}\sim 2\times 10^{42}\ erg\ s^{-1}$.

  * •

We estimated a neutral outflow rate of $\rm\dot{M}_{out,neut}\sim 30\
\textrm{M}_{\odot}/\textrm{yr}$ based on the detection of blueshifted velocity
components of the NaID doublet on the nuclear region (see Fig. 11). Although
the estimate is uncertain due to the unknown extent of the absorbing material,
the results confirm the estimate by Perna et al. (2015) based on MgII
absorption.

  * •

Including the contribution of the molecular outflow from B18 in the ionised
and neutral phases derived by NIRSpec, we derived a total mass outflow rate of
$\rm\dot{M}_{out,tot}\sim 110\ \textrm{M}_{\sun}/\textrm{yr}$. This
corresponds to a depletion time of just
$\rm\tau_{depl}=\textrm{M}_{gas}/(\textrm{SFR}+\dot{M}_{out,tot})\sim 30$ Myr,
which supports the key role of AGN feedback in this source.

These observations add to the growing evidence that the radio emission and
ionised gas outflows in active galaxies are connected. The also show the
relevant role of radio jets in AGN feedback on the host galaxy ISM even at
modest luminosities (see e.g. Jarvis et al. 2019, 2021, Venturi et al. 2021,
Girdhar et al. 2022).

Altogether, the NIRSpec observations analysed in this work demonstrate that
NIRSpec is capable of providing unprecedented information on the role of
quasar feedback at high redshift. Future NIRSpec and MIRI datasets will be
thus fundamental not only to enable our first glance at AGN feedback at z$>$4,
but also to overcome the current limitations of ground-based IFU in the study
of QSOs at z$\sim 1.5-2.5$. They will provide the required major leap in
resolution and sensitivity to obtain a sharper answer regarding the role of
quasar feedback on galaxy evolution at the epoch of peak of AGN activity in
the Universe, the golden epoch for AGN feedback.

###### Acknowledgements.

GC, GT, MP, MB, AM, FM, FB and GV acknowledge the support of the INAF Large
Grant 2022 ”The metal circle: a new sharp view of the baryon cycle up to
Cosmic Dawn with the latest generation IFU facilities”. GC, MB, AM and MB
acknowledge support from PRIN MIUR project “Black Hole winds and the Baryon
Life Cycle of Galaxies: the stone-guest at the galaxy evolution supper”,
contract #2017PH3WAT. FM, AM, and GC acknowledge support from the INAF Large
Grant 2022 ”Dual and binary SMBH in the multi-messenger era”. SA and MP
acknowledge support from the research project PID2021-127718NB-I00 of the
Spanish Ministry of Science and Innovation/State Agency of Research
(MICIN/AEI). GV acknowledges support from ANID program FONDECYT Postdoctorado
3200802. HÜ gratefully acknowledges support by the Isaac Newton Trust and by
the Kavli Foundation through a Newton-Kavli Junior Fellowship. Funded by the
European Union (ERC, WINGS, 101040227). Views and opinions expressed are
however those of the author(s) only and do not necessarily reflect those of
the European Union or the European Research Council Executive Agency. Neither
the European Union nor the granting authority can be held responsible for
them.

## References

  * (1) Arribas, S., Colina, L., Bellocchi, E., et al. 2014, A&A, 568, A14.
  * (2) Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5
  * (3) Balmaverde, B., Marconi, A., Brusa, M., et al. 2016, A&A, 585, A148.
  * (4) Balmaverde, B., Capetti, A., Marconi, A., et al. 2018, A&A, 619, A83.
  * (5) Bessiere, P. S. & Ramos Almeida, C. 2022, MNRAS, 512, L54.
  * (6) Böker T., Arribas S., Lützgendorf N., Alves de Oliveira C., Beck T. L., Birkmann S., Bunker A. J., et al., 2022, A&A, 661, A82.
  * (7) Brusa, M., Civano, F., Comastri, A., et al. 2010, ApJ, 716, 348
  * (8) Brusa M., Bongiorno A., Cresci G., Perna M., Marconi A., Mainieri V., Maiolino R., et al., 2015, MNRAS, 446, 2394.
  * (9) Brusa, M., Feruglio, C., Cresci, G., et al. 2015, A&A, 578, A11.
  * (10) Brusa, M., Cresci, G., Daddi, E., et al. 2018, A&A, 612, A29.
  * (11) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682.
  * (12) Cano-Díaz, M., Maiolino, R., Marconi, A., et al. 2012, A&A, 537, L8
  * (13) Carniani, S., Marconi, A., Maiolino, R., et al. 2015, A&A, 580, A102.
  * (14) Carniani, S., Marconi, A., Maiolino, R., et al. 2016, A&A, 591, A28.
  * (15) Cappellari, M. 2017, MNRAS, 466, 798.
  * (16) Cicone, C., Maiolino, R., Sturm, E., et al. 2014, A&A, 562, A21.
  * (17) Cimatti, A., Daddi, E., & Renzini, A. 2006, A&A, 453, L29.
  * (18) Coil, A. L., Aird, J., Reddy, N., et al. 2015, ApJ, 801, 35.
  * (19) Costa, T., Sijacki, D., & Haehnelt, M. G. 2015, MNRAS, 448, L30.
  * (20) Cresci, G., Hicks, E. K. S., Genzel, R., et al. 2009, ApJ, 697, 115
  * (21) Cresci, G., Vanzi, L., Sauvage, M., et al. 2010, A&A, 520, A82.
  * (22) Cresci, G., Mainieri, V., Brusa, M., et al. 2015, ApJ, 799, 82.
  * (23) Cresci, G., Marconi, A., Zibetti, S., et al. 2015, A&A, 582, A63.
  * (24) Cresci, G., Vanzi, L., Telles, E., et al. 2017, A&A, 604, A101.
  * (25) Cresci, G. & Maiolino, R. 2018, Nature Astronomy, 2, 179.
  * (26) Delvecchio, I., Smolčić, V., Zamorani, G., et al. 2017, A&A, 602, A3.
  * (27) Dugan, Z., Gaibler, V., & Silk, J. 2017, ApJ, 844, 37.
  * (28) Duggal, C., O’Dea, C., Baum, S., et al. 2021, Astronomische Nachrichten, 342, 1087.
  * (29) Di Teodoro, E. M. & Fraternali, F. 2015, MNRAS, 451, 3021.
  * (30) Ferland, G. J., Korista, K. T., Verner, D. A., et al. 1998, PASP, 110, 761.
  * (31) Fluetsch, A., Maiolino, R., Carniani, S., et al. 2019, MNRAS, 483, 4586.
  * (32) Gaibler, V., Khochfar, S., & Krause, M. 2011, MNRAS, 411, 155.
  * (33) Gaibler, V., Khochfar, S., Krause, M., et al. 2012, MNRAS, 425, 438.
  * (34) Girdhar, A., Harrison, C. M., Mainieri, V., et al. 2022, MNRAS, 512, 1608
  * (35) Förster Schreiber, N. M., Übler, H., Davies, R. L., et al. 2019, ApJ, 875, 21
  * (36) Giroletti, M., Panessa, F., Longinotti, A. L., et al. 2017, A&A, 600, A87.
  * (37) Harrison, C. M., Thomson, A. P., Alexander, D. M., et al. 2015, ApJ, 800, 45.
  * (38) Harrison, C. M., Alexander, D. M., Mullaney, J. R., et al. 2016, MNRAS, 456, 1195.
  * (39) Herrera Ruiz, N., Middelberg, E., Deller, A., et al. 2017, A&A, 607, A132.
  * (40) Herrera Ruiz, N., Middelberg, E., Deller, A., et al. 2018, A&A, 616, A128.
  * (41) Huško, F. & Lacey, C. G. 2022, arXiv:2208.09393
  * (42) Ishibashi, W., & Fabian, A. C. 2012, MNRAS, 427, 2998
  * (43) Jakobsen P., Ferruit P., Alves de Oliveira C., Arribas S., Bagnasco G., Barho R., Beck T. L., et al., 2022, A&A, 661, A80.
  * (44) Jarvis, M. E., Harrison, C. M., Thomson, A. P., et al. 2019, MNRAS, 485, 2710.
  * (45) Jarvis, M. E., Harrison, C. M., Mainieri, V., et al. 2021, MNRAS, 503, 1780.
  * (46) Kaasinen, M., Bian, F., Groves, B., et al. 2017, MNRAS, 465, 3220.
  * (47) Kakkad, D., Mainieri, V., Vietri, G., et al. 2020, A&A, 642, A147.
  * (48) Kashino, D., Silverman, J. D., Sanders, D., et al. 2017, ApJ, 835, 88.
  * (49) Kewley, L. J., Groves, B., Kauffmann, G., et al. 2006, MNRAS, 372, 961.
  * (50) Kewley, L. J., Dopita, M. A., Leitherer, C., et al. 2013, ApJ, 774, 100.
  * (51) King, A. 2005, ApJ, 635, L121
  * (52) Koekemoer, A. M., Aussel, H., Calzetti, D., et al. 2007, ApJS, 172, 196
  * (53) Lamareille, F. 2010, A&A, 509, A53.
  * (54) Law, D. R., Ji, X., Belfiore, F., et al. 2021, ApJ, 915, 35.
  * (55) Lusso, E., Comastri, A., Simmons, B. D., et al. 2012, MNRAS, 425, 623
  * (56) Madau, P. & Dickinson, M. 2014, ARA&A, 52, 415.
  * (57) Maiolino, R., Russell, H. R., Fabian, A. C., et al. 2017, Nature, 544, 202.
  * (58) Marasco, A., Cresci, G., Nardini, E., et al. 2020, A&A, 644, A15.
  * (59) Marasco, A., Belfiore, F., Cresci, G., et al. 2022, arXiv:2209.02726
  * (60) Masters, D., McCarthy, P., Siana, B., et al. 2014, ApJ, 785, 153.
  * (61) Matzeu, G. A., Brusa, M., Lanzuisi, G., et al. 2022, arXiv:2212.02960.
  * (62) Mingozzi, M., Cresci, G., Venturi, G., et al. 2019, A&A, 622, A146.
  * (63) Morganti, R., Oosterloo, T., Oonk, J. B. R., et al. 2015, A&A, 580, A1.
  * (64) Mukherjee, D., Bicknell, G. V., Wagner, A. Y., et al. 2018, MNRAS, 479, 5544.
  * (65) Nardini, E., Reeves, J. N., Gofford, J., et al. 2015, Science, 347, 860.
  * (66) Perna, M., Brusa, M., Cresci, G., et al. 2015, A&A, 574, A82.
  * (67) Perna, M., Arribas, S., Catalán-Torrecilla, C., et al. 2020, A&A, 643, A139.
  * (68) Pounds, K. A., King, A. R., Page, K. L., et al. 2003, MNRAS, 346, 1025.
  * (69) Renzini, A. 2006, ARA&A, 44, 141.
  * (70) Rupke, D. S., Veilleux, S., & Sanders, D. B. 2002, ApJ, 570, 588.
  * (71) Rupke, D. S. N., & Veilleux, S. 2011, ApJ, 729, L27.
  * (72) Sanders, R. L., Shapley, A. E., Kriek, M., et al. 2016, ApJ, 816, 23.
  * (73) Sargent, M. T., Daddi, E., Béthermin, M., et al. 2014, ApJ, 793, 19.
  * (74) Schinnerer, E., Smolčić, V., Carilli, C. L., et al. 2007, ApJS, 172, 46.
  * (75) Scholtz, J., Harrison, C. M., Rosario, D. J., et al. 2021, MNRAS, 505, 5469.
  * (76) Silk, J. 2013, ApJ, 772, 112
  * (77) Smolčić, V., Novak, M., Bondi, M., et al. 2017, A&A, 602, A1.
  * (78) Spergel, D. N., Verde, L., Peiris, H. V., et al. 2003, ApJS, 148, 175.
  * (79) Steidel, C. C., Rudie, G. C., Strom, A. L., et al. 2014, ApJ, 795, 165.
  * (80) Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179.
  * (81) Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157.
  * (82) Tozzi, G., Cresci, G., Marasco, A., et al. 2021, A&A, 648, A99.
  * (83) Vardoulaki, E., Jiménez Andrade, E. F., Karim, A., et al. 2019, A&A, 627, A142.
  * (84) Venturi, G., Cresci, G., Marconi, A., et al. 2021, A&A, 648, A17.
  * (85) Wagner, A. Y., Bicknell, G. V., & Umemura, M. 2012, ApJ, 757, 136.
  * (86) Wylezalek, D., Vayner, A., Rupke, D. S. N., et al. 2022, ApJ, 940, L7.
  * (87) Woo, J.-H., Son, D., & Bae, H.-J. 2017, ApJ, 839, 120.