© The Authors 2023

1. Introduction

Supermassive black holes (SMBHs) in the hearts of galaxies are believed to interact with their host galaxy. The enormous amount of energy released by an active galactic nucleus (AGN) is injected into the interstellar medium (ISM) where it effectively delays the gas cooling. This deprives the galaxy of its cold gas reservoir (Zubovas & King 2012) and eventually leads to star formation (SF) quenching, an effect that is referred to as ‘negative feedback’. Such AGN-induced negative feedback is an important ingredient to reproduce galaxy properties in cosmological hydrodynamical simulations (e.g. Crain et al. 2015; Dubois et al. 2016; Weinberger et al. 2017) and semi-analytic models (e.g. Somerville et al. 2008; Fontanot et al. 2020). In contrast, AGNs have also been suggested to locally promote the SF in its host galaxy by enhancing the gas pressure in the ISM (‘positive feedback’, Silk 2013; Cresci et al. 2015).

The physical mechanisms through which AGN feedback acts are poorly constrained by observations. Open questions include when and how much energy is released during the AGN phase, how the energy is transported from the accretion disc to galaxy scales, and how it couples to the surrounding ISM. One important signature of AGNs are powerful outflows, which are a product from the energy released by the accreting SMBH and its interaction with the ISM. Such AGN-driven outflows have a multi-phase nature. In the X-ray, they can be identified in the form of fast radiation-driven winds, which originate from scales of the BH accretion disc (Tombesi et al. 2010, 2013; Gofford et al. 2013). However, atomic, molecular, and ionised gas outflows are typically observed at distances from parsecs up to several kilo-parsecs from the galaxy centre. Ionised outflows are often identified from the systematically asymmetric emission line shape of [O III]λλ4960,5007. The blue shoulder is usually interpreted as a bi-conical outflow of AGN-ionised gas where the receding side is obscured by dust from the host galaxy (e.g. Heckman et al. 1981; Bischetti et al. 2017; Bae & Woo 2016; Wylezalek et al. 2020). Both observations (Silk & Rees 1998; King 2003; Holt et al. 2006, 2008; Cicone et al. 2014) and hydrodynamical simulations (e.g. Di Matteo et al. 2005; Gaspari et al. 2020) suggest that multi-phase outflows represent a key channel through which the AGN feedback acts on the host galaxy. Moreover, such feedback can significantly stimulate turbulent condensation, leading to gas precipitation and chaotic cold accretion (CCA; e.g. Gaspari et al. 2019).

Due to the compactness of the AGN-dominated region from which the outflows are launched, it is still a matter of debate how the evolution of the central BH impacts the host galaxy on different scales. An increasing number of studies provide a detailed and spatially resolved analysis of ionised gas outflows in local AGNs (e.g. Greene et al. 2011; Husemann et al. 2016a, 2019; Revalski et al. 2018; Venturi et al. 2018; Riffel et al. 2020). One big challenge is to determine the outflow location, orientation, and intrinsic size. The works of Greene et al. (2011), Liu et al. (2013), Harrison et al. (2014) and McElroy et al. (2015) have suggested that galaxy-scale outflows might be prevalent in the majority of luminous AGNs. In contrast, a recent study by Singha et al. (2022) has shown that among the luminous AGNs in the representative sample of the Close AGN Reference Survey (CARS, Husemann et al. 2022), a large fraction (64%) of them have ionised gas outflows with an extent of less than 100 pc. This disagreement can partially be explained by beam smearing in seeing-limited observations (Husemann et al. 2016b; Villar-Martín et al. 2016; Davies et al. 2019). Another problem is that the inferred outflow properties such as electron density, masses, and kinetic energies are poorly constrained as they depend on the outflow geometry, which is often unknown in luminous type 1 AGNs (e.g. Rakshit & Woo 2018).

Constraining the fundamental spatial, energy, and timescales is crucial to understand how AGN-driven multi-phase outflows are launched and how they couple to the host galaxy ISM. For AGNs in the Sloan Digital Sky Survey (SDSS), Mullaney et al. (2013) and Woo et al. (2016) have shown that outflows occur more frequently in AGNs with a high Eddington ratio. In this regime, radiation pressure from the accretion disc couples to the dense ISM and can drive the outflow (Fabian 2012), suggesting that ionised gas outflows are launched during the AGN phase. Among the general AGN population, narrow line Seyfert 1 (NLS1) galaxies show low BH masses at high accretion rates (e.g. Boroson 2002; Grupe et al. 2010, review by Komossa & Xu see review by 2007). NLS1s often host ultra-fast outflows (UFOs; e.g. Gupta et al. 2013; Parker et al. 2017; Reeves et al. 2018; Kosec et al. 2018; Reeves & Braito 2019; Xu et al. 2022). Furthermore, NLS1s are located at the extreme end of the Eigenvector 1 (EV1) correlation, which involves the widths and strengths of [O III]λ5007, Hβ, and Fe II emission lines (Boroson & Green 1992). Wang et al. (1996) extended the EV1 plane by the soft X-ray photon index where NLS1s typically exhibit a prominent, soft X-ray excess (Gliozzi & Williams 2020). The EV1 correlation is thought to be primarily driven by the Eddington ratio (Sulentic et al. 2000), which suggests that there might be a connection between the BH accretion rate and the launching mechanism of a multi-phase outflow. For a comprehensive understanding, we need to resolve the multi-phase outflow both spatially and in its kinematic components.

In this work we present a detailed analysis of the ionised gas outflow in the centre of Mrk 1044. Mrk 1044 is a nearby, luminous (L_bol = 3.4 × 10⁴⁴ erg s⁻¹, Husemann et al. 2022) NLS1 with a stellar mass of $log(M_{\star }/M_{\odot })=9.{92}_{-0.12}^{+0.17}$ (Smirnova-Pinchukova et al. 2022). Its central engine is powered by a BH with a reverberation mapped mass of M_• = 2.8 × 10⁶ M_⊙ (Du et al. 2015). Mrk 1044 shows a soft X-ray excess, which Mallick et al. (2018) explained by relativistic reflections from a high density accretion disc. From the several narrow absorption lines that are present in Mrk 1044’s XMM-Newton spectra, Krongold et al. (2021) identified four kinematically distinct UFOs. Their velocity and density structure suggest that the two light absorbers may originate from the same multi-phase outflow.

The results presented in this work are based on the analysis of Mrk 1044’s host galaxy presented in Winkel et al. (2022, hereafter Paper I) where we found that the host galaxy SF is concentrated in a circumnuclear ellipse (CNE) at ∼300 pc from the centre. Furthermore, we detected tentative signatures of ionised gas channelling towards the centre, which might be a signature of ongoing BH fuelling. Despite the high BH accretion rate, the host galaxy does not exhibit any signs of disturbance even at a short distance from the nucleus (∼160 pc). Based on the velocity field of the stellar component, we have constrained Mrk 1044’s systemic velocity to cz = 4913.4 ± 0.2 km s⁻¹, corresponding to a redshift of z = 0.0164. In the very centre of Mrk 1044, we detected a spatially resolved ionised gas component that is blue-shifted, which indicates an ionised gas outflow. In this work we use optical IFU data, UV spectroscopic data, and radio imaging to constrain the properties of Mrk 1044’s multi-phase outflow. Combining the observations, we aim to attain a detailed understanding of the outflow location, geometry, and its physical driver. Throughout this paper we assume a flat Λ cold dark matter cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_M = 0.3, and Ω_Λ = 0.7. In this framework, 1″ corresponds to 333 pc at the galaxy redshift, where Mrk 1044’s associated luminosity distance is 70.0 Mpc.

2. Observations and data reduction

2.1. Optical IFU observations

We use integral field spectroscopic data obtained with the adaptive optics assisted narrow field mode (NFM-AO) of the Multi Unit Spectroscopic Explorer (MUSE, Bacon et al. 2010, 2014) at the Very Large Telescope (VLT). The data reduction is described in Paper I. The reduced data cube consists of 369 × 378 spaxels, corresponding to a FOV of $9\stackrel{″}{.}23\times 9\stackrel{″}{.}45$. The spectral resolution is almost constant with a full width at half maximum (FWHM) of 2.54 ± 0.10 Å across the entire wavelength range 4750 Å–9350 Å corresponding to 160.4 km s⁻¹ and 81.5 km s⁻¹ at the blue and the red end of the spectrum, respectively. In Paper I we describe the modelling of the adaptive-optics-shaped point spread function (PSF) using a hybrid approach between the empirical PSF and the analytical PSFAO19 model from Fétick et al. (2019). The spatial resolution measured from the PSF extracted at the broad Hβ emission is 89 mas which corresponds to 30 pc in the galaxy system. We have deblended the AGN from the host emission in an iterative process that is described in Paper I. In this work, we aim to extract the properties of the outflow from the original data cube that contains the blended emission from the AGN and the host.

2.2. UV spectroscopy

We employ an archival UV spectrum of Mrk 1044 that was acquired with the Space Telescope Imaging Spectrograph (STIS) on the Hubble Space Telescope (HST). This data has first been presented in Fields et al. (2005) where the authors identified two outflowing absorbers with super-solar metallicity. Since we have a robust measure of Mrk 1044’s systemic velocity from the stellar rotational field, we re-analyse the outflow kinematics. We retrieve the archival data from the Hubble Legacy Archive¹ and stack the two data sets o8k4010-50 and o8k4010-60 taken with the G140M grating. The target was observed on UTC 2003 June 28 with the 52X0.2 aperture centred on Mrk 1044’s nucleus. The two spectra were acquired with an exposure time of 1294 s and 1440 s respectively and have a resolution of 0.053 Å px⁻¹. The stacked spectrum extents from 1194.57 Å to 1249.10 Å with a central wavelength of 1222 Å.

2.3. Radio imaging

Mrk 1044 has been observed with Karl G. Jansky Very Large Array (VLA) in the C-band (6 GHz) on October 30, 2016 and X-band (10 GHz) on January 15, 2017. The JVLA observations were taken in the A configuration which has a maximum baseline of 36.4 km and a minimum baseline of 0.68 km. The A configuration observations are sensitive to emission on scales up to 9″ in C-band and up to 5$\stackrel{″}{.}$3 at X-band. The three C-band scans have an integration time of 245 s each resulting in a total integration time 735 s. The two X-band scans have an integration time of 325 s each. Together with two existing additional scans (330 s), the X-band has a total integration time of 1310 s. 3C48 and 3C138 were used as flux density calibrators. We reduced the JVLA data using the Common Astronomy Software Applications (CASA, McMullin et al. 2007) pipeline version 6.2.1.7 and the dedicated CASA tools for VLA observations.

3. Analysis and results

3.1. Optical AGN spectrum

In order to understand the ISM properties of Mrk 1044’s outflow on nuclear scales, we aim to extract the ionised gas kinematic components present in the innermost few parsecs. In Paper I we have described the iterative deblending of AGN and host emission. For the following analysis we use the AGN spectrum from the deblended AGN data cube which contains the point-like emission from Mrk 1044’s nucleus exclusively. Due to the finite width of the PSF, it probes the innermost 30 pc around the nucleus. The spectrum in the observed frame is shown in Fig. 1 for both the Hβ+[O III] and the Hα+[N II]+[S II] regions. It exhibits prominent emission lines Hβ/Hα and [O III]λλ4959,5007 and especially from Fe II, whereas the [N II]λλ6548,6583 doublet is barely visible.

Fig. 1. Modelling Mrk 1044’s AGN spectrum in the Hβ–[O III] region (left) and Hα–[N II] region (right). To reproduce the observed spectrum (black) we use a multi-Gaussian model for Hβ (blue), [O III] (green), Fe II (purple), Hα (blue, right panel) and [N II] (red). The best-fit spectrum is shown as a red line and well reproduces the prominent emission lines, including the blue shoulder of the [O III] narrow line.

We simultaneously model the AGN spectrum in the Hβ+[O III] and the Hα+[N II] region. Our model consists of several Gaussian components and a linear approximation of the local continuum in each of the windows. To reproduce the Hβ and Hα line shape, we require a broad+medium component and a narrow core component. The broad and medium component both have a kinematically tied component in the prominent broad Fe IIλλ4924, 5018 doublet. For the forbidden doublet emission lines [N II]λλ6548,6583 and [O III]λλ4959,5007 we couple the flux ratios to the theoretical prediction of 2.96 (Storey & Zeippen 2000; Dimitrijević et al. 2007). To reproduce the narrow [O III] line shape we require an additional blue-shifted ‘wing’ component which has a kinematically tied counterpart in Hβ and Hα. To find the best-fit parameters of the model we minimise the χ² of the residuals with a Levenberg-Marquardt algorithm. We estimate the uncertainties of the parameters by modelling 1000 artificial spectra generated from the 1σ flux density errors in the AGN spectrum. We note that the identification of the components in the emission line spectrum can also be reproduced with software packages that involve a Fe II-template to account for the strong Fe II contribution in the Hβ–[O III] window. However, the systematic uncertainty of the optional Fe II template provides a poorer description of the AGN spectrum than our multi-Gaussian model, especially in the [O III] region where we perform most of our diagnostics. We note further that the results presented in the following do not change within the uncertainties if the emission from ions with different ionisation potentials are not kinematically coupled. In particular, the results do not change qualitatively if we release the kinematic coupling between the Hβ and [O III] components and the Hα-wing and [N II] components respectively.

Figure 1 shows that the AGN spectrum in both spectral regions is well reproduced by our four component model. Considering the high signal-to-noise ratio (S/N), also the Hα line shape is well-reproduced close to its centre. The broad wings, however, are contaminated by the broad Fe IIλ6517, He I emission and Fe complex (cf. Fig. 7, Véron-Cetty et al. 2004) which are not included in our model. The kinematic parameters of the best-fit model of the AGN spectrum are summarised in Table 1. Both Hα+[N II] and Hβ+[O III] narrow core components have radial velocities significantly smaller than zero by v_r(Hβ + [O III]) = −144 ± 5 km s⁻¹ and v_r(Hα + [N II]) = −140 ± 7 km s⁻¹, indicating that the region from which the narrow line emission emerges is moving towards the observer. In addition to the blue-shifted narrow core component, we detect a wing component in the Hβ and [O III] emission with a velocity shift of −560 ± 20 km s⁻¹. This feature is often recovered in AGN spectra and is typically interpreted as an ionised gas outflow moving towards the observer, where the receding side is obscured by dust around the nucleus (e.g. Heckman et al. 1981; Boroson 2005; Zakamska et al. 2003; Woo et al. 2016).

3.2. Spectroastrometry

After having identified two outflowing systems in both Mrk 1044’s Hβ- and [O III] emission, we now aim to constrain their location and spatial extent. For this purpose, we have developed the Spectroastrometric analysis tool for Ionized gas Emission in Nearby AGN (Siena3D, Winkel & Singha 2022). Siena3D locates and traces kinematic components in 3D spectroscopic observations with high precision, using a spectroastrometric approach (for discussion see Bailey 1998; Gnerucci et al. 2010, 2011). The software package is available on GitHub² under a OSI³ approved licence. The concept behind Siena3D is described in Singha et al. (2022) where the technique has been applied and tested for the entire CARS sample of nearby unobscured AGNs. Here, we give a brief description of the concept.

(1) As a first step, we use the kinematic components identified in the AGN spectrum to generate normalised base spectra (top left panel of Fig. 2). All emission lines that belong to the kinematic components are combined in the respective base spectrum. Their flux ratios and kinematics are tied to what we measured in the highest signal-to-noise spaxel. This step involves the assumption that both line ratios and kinematics are constant throughout the structure that is traced (e.g. a compact fast moving ionised gas cloud or the BLR).

Fig. 2. Spectroastrometric analysis of Mrk 1044’s central region using the original MUSE NFM-AO data cube. The left panels show an arbitrary example spectrum that is picked from the small window around the nucleus shown in the panels on the right. To fit the spectrum, we keep the kinematics and line ratios of the kinematics fixed to that of the components found in the AGN spectrum. The panels on the right show the surface brightness within the central 150 mas for the kinematic components where the flux maps are normalised to their peak flux. Here we only show the measured light distribution from the BLR Hβ (left), [O III]-wing (middle), and narrow [O III]-core (right) component. From top to bottom the maps show the 2D light distributions, the best-fit model of the PSF and the residual maps. The red cross indicates the PSFAO19 centroid for the best-fit PSFAO19 fit to the BLR Hβ. The green and white crosses indicate the centroids for the PSFAO19 models to the surface brightness maps of the [O III]-wing and core respectively, in which the amplitude and position of the Moffat model were varied only. The offset of both [O III] core and wing component from the AGN position are smaller than 0.1 px.

(2) In the next step, we use the base spectra to fit the spaxels in a 14 × 14 pixel (0$\stackrel{″}{.}$35 × 0$\stackrel{″}{.}$35) aperture around the AGN position using a non-negative least squares optimiser. For one arbitrary spaxel, the best-fit model together with the contribution from the components are shown in the bottom left panel of Fig. 2. Our results do not depend on the exact size of the aperture, since the high S/N core of the PSF (FWHM = 1.6 px to 3.6 px) is captured entirely.

(3) From the best-fit model, we map the flux that originates from the individual kinematic components. We fit the PSF to the light distributions to locate their centroids relative to the AGN position. For this task we adopt the hybrid PSF model that involves the empirical PSF and the analytic PSFAO19 model for the adaptive-optics-shaped PSF measured in Paper I, using the analytic model from Fétick et al. (2019).

(4) Finally, we evaluate whether the light distribution deviates from a point-like emission to test whether the component has an intrinsic extended that is spatially resolved. The methodology is described in Appendix A where we constrain the projected extent of a resolved component (see Sect. 3.3).

The precision by which different components can be kinematically and spatially disentangled depends on multiple parameters including the size and sampling of the PSF, the accuracy of the PSF model, the spectral resolution, the velocity offset between the kinematic components and the S/N of the emission lines and their dispersion. In order to estimate the systematic uncertainty of our method specific to our data set, we have generated a set of mock data cubes by simulating the emission line profiles at different velocities, locations and sizes. We manually displaced the [O III]-wing component while keeping the location of the broad and core components fixed. At the S/N of [O III] core and wing component, the offset of both can be detected down to 0.1 px, corresponding to 2.5 mas or 0.85 pc in the galaxy system. Furthermore, we have simulated spatially resolved structure by convolving the PSF with a Gaussian kernel. For Mrk 1044 we can constrain the intrinsic extent of the outflow down to a size of 0.2 px (1.7 pc). At smaller scales, the precision is limited by the relatively low S/N of the wing component that is blended with the core component in each of the emission lines [O III]λλ4959,5007.

For Mrk 1044, the result of the spectroastrometric measurement in the Hβ–[O III]-region is shown in Fig. 2, with the corresponding centroid positions listed in Table 2. The broad component originates from the μpc-scale BLR that is unresolved and follows the PSF light profile by construction. Therefore, the residual maps of the broad component reflect the systematic uncertainty of the PSF. But also the residuals of the [O III]-core and wing component show no extended structure. Neither [O III] outflowing component has an intrinsic extent that is larger than 0.2 px (1.7 pc). Furthermore, both core and wing component are located at the exact AGN position within the systematic uncertainties. We therefore conclude that neither of the [O III]-outflows have a projected distance from the nucleus that is larger than 0.85 pc.

3.3. The spatially resolved outflow

During the analysis of the star-forming CNE in Paper I, we have reported a high surface brightness spot in the centre that is present after subtracting the AGN emission. In the following we refer to this feature as H0 as it is most pronounced in Hα and Hβ surface brightness maps. The left panel of Fig. 3 shows that in the host data cube the feature is constrained to the innermost $0\stackrel{″}{.}$5. In Paper I the analysis with a single Gaussian emission line model yielded blue-shift of −140 km s⁻¹ and a peak velocity dispersion that is significantly higher than that of the CNE (σ = 234 ± 9 km s⁻¹, see Fig. 7 in Paper I). The high velocity-offset with respect to the galaxy systemic velocity suggests that the ionised gas is outflowing. We detect the gas in each of the narrow emission lines Hβ, [O III], Hα and [N II] whereas we cannot trace a significant component [S II] given the low spectral signal-to-noise of the spatially resolved analysis.

Fig. 3. Origin of the forbidden [S II] emission in Mrk 1044’s centre. The left panel shows the Hα surface brightness map where the H0 outflow is located the central 0$\stackrel{″}{.}$5. The right panel shows the two component model for the inner 0$\stackrel{″}{.}$5 integrated spectrum where we fixed the host and outflow kinematics to what we found from the 3″ integrated spectrum. For the outflow, we do not find a significant contribution of [S II] and [N II] suggesting that the gas density lies above their critical densities.

Although the outflow appears circularly symmetric, it is not an artefact from the PSF subtraction. To better understand this we recall the PSF subtraction process. Hα and Hβ are the wavelengths at which we have extracted the empirical PSF. By definition, the empirical PSF contains all the point-like emission. After the PSF cube construction, our PSF subtraction method is an iterative routine to clean the host emission from the AGN contamination. Since we interpolate the PSF between Hα and Hβ, an inaccurate PSF model would cause much stronger residual features at wavelengths at which the empirical PSF was not directly measured, that is [O III]. Furthermore, the central feature is present after each iteration of the PSF subtraction process. In particular, it is present after the first iteration where the host-contaminated AGN-spectrum is subtracted. The circularly symmetric surface brightness profile can thus be explained by a compact source together with the effect of beam smearing. To understand the origin of this spatially resolved outflow, we constrain its extent as described Appendix A. Our analysis yields a projected intrinsic extent of 0.55 ± 0.10 px, corresponding to 4.6 ± 0.6 pc in the galaxy frame.

3.3.1. Multi-component analysis

The spectrum around the nucleus is a superposition of the emission from the host and that of the spatially resolved outflow H0. The two signals stem from intrinsically distinct structures, potentially located at different regions in the galaxy and having different ionisation conditions. To disentangle and analyse their signal, we follow the methodology described in Singha et al. (2022), which contains the following steps:

First, we fitted Hβ+O IIIλλ4959,5007 with a two-component (host and outflow) Gaussian model each in a within a 3″ aperture. This gives us a robust estimate of the integrated host and outflow kinematic parameters.

We then fitted each of the Hβ+[O III]λλ4959,5007 and Hα+[N II]λλ6548,83+[S II]λλ6716,31 with the two-component model, this time keeping the kinematic parameters fixed to the values retrieved in step 1. Furthermore, we couple the line ratios among the [O III] and [N II] doublet lines to the theoretical prediction of 2.96 (Storey & Zeippen 2000; Dimitrijević et al. 2007). Due to the relatively low S/N of the blended emission lines, this task could only be performed on the $0\stackrel{″}{.}5$-integrated spectrum.

The result of the process in shown in Fig. 3. Each of the strong emission lines is detected at a > 3σ confidence level with their kinematic parameters listed in Table 3. Except for the high signal-to-noise centre of the Hα emission line, the residual spectrum is uniformly distributed around zero, with a scatter of 0.22σ (FWHM), where σ is the formal uncertainty from the AGN-blended cube. We conclude that the multi-Gaussian model provides a good description of the nuclear host spectrum. We detect a wing component that is present in Hβ, Hα and [O III] exclusively. In contrast, the corresponding wing components in [N II] and [S II] are faint or missing such that we could only estimate upper flux limits. We constrain and discuss the underlying ionisation mechanism from the emission line ratios in Sect. 4.1.

3.4. Radio emission

The EVLA broad-band images in the C- and X-band are shown in Fig. 4. The peak frequencies are located at 5 GHz (C-band) and 10 GHz (X-band). At these frequencies, we detect the AGN with a peak flux of 758 ± 19 μJy beam⁻¹ (C-band) and 496 ± 10 μJy beam⁻¹ (X-band). The radio emission appears to be concentrated in the beam aperture in both bands. To quantify the size of the source, we use the CASA routine imfit to measure the AGN emission within a 10 px radius from the centre which is sufficient to capture the source’s entire emission. We find that the source is spatially unresolved by VLA and the centroid of its radio emission matches the location of the AGN within the uncertainties. Its integrated flux density amounts 714 ± 32 μJy in the C-band and 484 ± 18 μJy in the X-band respectively. Assuming that the spectrum between 5 GHz and 10 GHz is described by a power law S_ν ∝ ν^α, we estimate the spectral index from the flux density. Its uncertainty depends both on the statistical error from thermal noise of the image and the systematic uncertainties $\sqrt{{\sigma }_{\mathrm{rms}}^2+{\sigma }_{\mathrm{sys}}^2}$. Combining both, we use Eq. (1) from Sect. 3.2 from Ramírez-Olivencia et al. (2022) to estimate the uncertainty throughout our calculation. From Mrk 1044’s integrated flux densities we obtain a spectral index of α_int = −0.56 ± 0.23, whereas using the peak flux densities yields α_peak = −0.61 ± 0.16. Both values are consistent with each other, indicating once more that the source is spatially unresolved by VLA. As discussed by Ramírez-Olivencia et al. (2022), the peak flux is not affected by the systematic uncertainty from the different beam sizes between the C- and the X-band. We therefore consider the more robust result to be α_peak and adopt its value for the following analysis.

Fig. 4. Mrk 1044’s nuclear region seen in the optical (left) and radio emission (right). Panel (a) shows the surface brightness map of the Hα narrow line emission extracted as described in Paper I. The panels on the right show the broad-band EVLA images of the radio emission in the C-band (b) and X-band (c) with the corresponding FWHM-contour of the beam shown in the lower left corner respectively. The contours of the star forming CNE encompassing the nucleus are shown as black contours (see Paper I). In both C- and X-band the nuclear radio emission is unresolved by VLA.

3.5. Ly-α absorption

Fields et al. (2005) have identified two outflowing systems in Mrk 1044’s UV spectrum from Ly-α, N Vλλ1239,1243 and C IVλλ1548,51 absorption. Using the line strengths and a photo-ionisation model they estimated several times solar metallicities for both absorbing clouds and concluded that the absorbers were ejected from the BLR local to the galaxy. In contrast to Fields et al. (2005), we have used a kinematic model of a thin rotating disc to describe the stellar rotation in Paper I. This represents a robust measure for the galaxy systemic velocity cz = 4913.4 ± 0.2 km s⁻¹ which allows for a self-consistent comparison between Mrk 1044’s AGN, outflow and host galaxy properties after re-analysing the kinematics of the two UV absorbing systems.

We model Mrk 1044’s UV spectrum in the Ly-α-window using a set of Gaussians and a linear approximation of the continuum. To reproduce the broad emission line shape, we require three components whereas a single Gaussian is sufficient to model the absorption line shapes. To find the best-fit parameters, we use the same technique that is described in Sect. 3.1. The result of this process is shown in Fig. 5. For the absorbers that we employ for the following discussion, the kinematic parameters for the individual components are listed in Table 3. The broad component appears blue-shifted in the spectrum since Ly-α is a resonant line where the photons gain momentum by scattering on the fast-moving BLR clouds. Another reason might be that the BLR geometry might not be entirely symmetric. If the receding side of the BLR is obscured by dust, the observed Ly-α emission line is skewed towards the blue.

Fig. 5. Modelling Mrk 1044’s UV spectrum in the Ly-α window. We use a multi-Gaussian model for the Ly-α emission (blue), and the absorbing systems (orange). The best-fit model (red line) reproduces the spectrum (black line) well in the selected wavelength range. Since the Ly-α systems L1 and L2 are accompanied by high metallicity absorption lines at the same velocity (Fields et al. 2005), the absorbing clouds must be local to the galaxy.

There are two absorbers with velocities −1937 ± 6 km s⁻¹ and −2736 ± 4 km s⁻¹ which can be explained by foreground extinction from the Ly-α forest (Fields et al. 2005). The different systematic velocity explains the nature of the Ly-α absorbing system L0 that is located at the galaxy rest-frame with v_r(L0)=0.9 ± 3.8 km s⁻¹. Most importantly, we recover the two outflowing systems L1 and L2 for which we measure systematically higher velocities for the absorbing features compared to Fields et al. (2005). Both systems are local to the galaxy as they show corresponding absorption from N V and C IV. We measure the outflow velocity of the absorbing systems L1 and L2 as v_r(L1)= − 278 ± 43 km s⁻¹ and v_r(L2)= − 1118 ± 2 km s⁻¹, respectively.

4. Discussion

There are multiple outflowing systems present in Mrk 1044’s outflow, both in optical emission and UV absorption. In Table 3 we summarise their kinematic properties and their distance from the nucleus. In the following we discuss the ionised gas excitation mechanism close to the galaxy nucleus and limitations of the diagnostics. Furthermore, we discuss whether we can link the different outflow velocities between [O III]-emitters and Ly-α-absorbers using a geometric model. Finally, we estimate the mass outflow rates and energetics and discuss the future impact of the multi-phase outflow on the host galaxy star formation.

4.1. Nuclear excitation mechanism

In Paper I we have shown that Mrk 1044 has a high star formation rate that is concentrated in a circumnuclear ellipse (CNE). The 3″-integrated emission line ratios have shown that all spaxels in Mrk 1044’s CNE are located in the star forming regime. But also in the very centre at the location of the outflow, the line ratios seem to favour SF as the dominant excitation mechanism. However, after disentangling the host emission from the resolved outflow H0, a closer look onto the forbidden line ratios suggests that the situation is more complicated.

4.1.1. Nuclear star formation hidden by the outflow

From the multi-component model of the central region in Sect. 3.3 we compile the BPT diagnostic diagram that is shown in Fig. 6. It contains the data from the spatially resolved 8 × 8 co-added spaxels that belong to the star-forming CNE (see Paper I) together with the 0$\stackrel{″}{.}$5-integrated values for the host and outflow component. By definition, the host component is located in the galaxy rest-frame. It occupies the same area as the resolved CNE-spaxels. Although the host component may contain a significant amount of beam-smeared emission from the CNE, the host [S II] emission peaks in the centre, which implies that it cannot be beam-smeared emission alone. Thus, even within the central 0$\stackrel{″}{.}$5 there must be gas present at relatively low densities n_e <  4 × 10³ cm⁻³ the critical density of [S II] (De Robertis & Osterbrock 1986).

Fig. 6. BPT diagnostics of the extended outflow H0. The grey data points correspond to the 8 × 8 co-added spaxels that belong to the star-forming CNE (see Paper I). The stars represent the two components that are present in the $0\stackrel{″}{.}5$-integrated spectrum that is shown in Fig. 3. The host component (red star) is consistent with SF-classification of the line ratios recovered in the CNE. For the outflowing component (blue star), both [O III] and [N II] emission might be suppressed due to collisional de-excitation in the high density gas. The BPT classification of the outflow as SF-ionised might therefore be misleading.

Combining the above findings, it indicates that the host galaxy star-formation in the CNE is not limited to the distance of ∼300 pc from the nucleus. It extends down to < 0$\stackrel{″}{.}$5 (< 160 pc) from Mrk 1044’s nucleus, but is heavily blended by the beam-smeared emission of the ionised gas outflow. Following the same procedure as in Paper I, we use the $0\stackrel{″}{.}$5-integrated line flux from the host galaxy to estimate the $\mathrm{SFR}=5_{-3}^{+5}\times {10}^{-2}{M}_{\odot }$ yr⁻¹ that is ‘hidden’ by the beam-smeared emission from the luminous ionised gas outflow.

4.1.2. High outflow gas densities

Compared to the host components, the line fluxes associated with the outflowing ‘wing’-component H0 have a larger systematic uncertainty in the region around the narrow Hα+[N II]. This is caused by the blended emission from host and outflow components between [N II] and Hα and the relatively low signal-to-noise. While we detect the [O III]-wing, the kinematically corresponding [N II] emission is weak such that we could only estimate an upper limit from the level of the spectral noise (see Sect. 3.3). The strong line ratios appear to unambiguously locate Mrk 1044’s ionised gas outflow in the SF regime of the BPT diagnostic. We argue, however, that the BPT-classification of the outflow is misleading.

The total ionised gas mass ejected by the outflow from < 4.6 pc is 10³ M_⊙ yr⁻¹. Comparing with the relation for SF-driven outflows, the mass outflow rate exceeds the prediction by four orders of magnitude (Fluetsch et al. 2019; Stuber et al. 2021, see also Sect. 4.3). We conclude that Mrk 1044’s outflow cannot be star-formation driven but must be powered by the radiation field of the luminous AGN. In this scenario, however, the AGN photo-ionisation is expected to dominate the emission line ratios.

A physical explanation that resolves the contradiction might be related to the compactness of the outflow. Considering that it has not yet expanded to host-galaxy scales beyond a few pc, the entire gas must be contained in a significantly smaller volume compared to the kiloparsec-scale outflows of the extended narrow line region (ENLR). At less than a few pc distance from the nucleus we therefore expect significantly higher outflow gas densities. If it exceeds the critical density of [N II] and [O III], the bulk of the ions collisionally de-excite and thus suppress the forbidden emission lines. Assuming an electron temperature of T_e = 10⁴ K we estimate the electron density to be between the critical densities of [N II] and [O III] respectively, that is 8.7 × 10⁴ cm⁻³ <  n_e <  7 × 10⁵ cm⁻³ (De Robertis & Osterbrock 1986). In this scenario the BPT classification of the outflow as SF-ionised can be explained by the high gas density, even if the actual excitation mechanism is related to the AGN radiation field. Therefore, the BPT diagnostic is not a reliable tool to constrain the underlying ionisation mechanism in the presence of high density ionised gas dominating an unresolved structure. Combining the above arguments we suggest that Mrk 1044’s ionised gas outflow is likely AGN-driven.

4.2. Exploring a geometric model for the unresolved outflow component

In Sects. 3.1 and 3.5 we have described that in each [O III]-emission and Ly-α absorption, two spatially unresolved outflowing systems are present. A summary of the derived kinematic parameters is listed in Table 3. We notice that velocity ratios among the [O III]-emitters (v_O2/v_O1 = 3.9) and the Ly-α absorbers (v_L2/v_L1 = 4.0) are similar. Furthermore, the velocity ratios between the faster and the slower components respectively match (v_L2/v_O2 = 2.0, v_L1/v_O1 = 1.9). Motivated by this constant factor, we aim to identify the systems with each other and discuss whether the integrated [O III] emission line shape can be explained through a geometric alignment.

In order to constrain the geometry, we use a forward modelling approach to predict the shape of the [O III] emission lines based on an expanding shell geometry. We parameterise the model by the inclination i of the ionised gas outflow velocity with respect to the line-of-sight, the ‘narrowness’ α/θ₀ that describes the radial dependence of the shells’ intrinsic luminosity profile in units of the half-opening angle θ₀, and the flux ratio between the two shells r. We use a forward modelling approach to constrain the parameters by fitting the [O III] emission line shape. A detailed description of the method is described in the Appendix B.

The geometric model predicts highly inclined shell caps (i = 79 ± 8°) which are centrally concentrated α/θ₀ = 0.2 ± 0.1. Within the uncertainties, the result on the inclination is independent of the exact value of the BLR opening angle and how centrally concentrated the outflow is (0.01 <  α/θ₀ <  1).

4.2.1. Expanding shell model versus dust torus model

The AGN unification model predicts that an active nucleus possesses a compact strong radiation source around the SMBH that is surrounded by a high velocity emission-line region, the BLR. The relatively small scatter of the half-opening angles of the AGN dust-tori of ∼60° (Ricci et al. 2017) suggests that the optical classification is a reliable indicator for the orientation of the obscuring dust torus. For Mrk 1044’s multi-phase outflow, the expanding shell model fails to reproduce the [O III] emission line shape for i ⪅ 60°. At larger inclination, however, the dust obscuration from the torus encompassing the nucleus is expected to block the view into the BLR, which contradicts the type-1 nature of Mrk 1044.

One way to resolve the tension is that Mrk 1044’s BLR obscuration may differ from the nearby low-luminosity Seyfert galaxies for which the AGN unification model was originally developed. From the Swift-Burst Alert Telescope (BAT) selected local AGN sample, Ricci et al. (2017) reported that the AGN obscured fraction is mainly determined by the Eddington ratio λ_Edd. They concluded that the main physical driver of the torus diversity is the accretion rate, which regulates the torus covering factor by means of radiation pressure. A clumpy torus (Schartmann et al. 2008; Stalevski et al. 2019) provides a more accurate description of the relationship between X-ray and mid-IR luminosities (Nenkova et al. 2008; Stalevski et al. 2012; Hönig & Kishimoto 2017) and also been resolved for individual systems with high-resolution mid-IR imaging (e.g. Isbell et al. 2021). It is unclear whether the same orientation-based unification model involving a smooth torus is directly applicable to a highly accreting system like Mrk 1044. Nonetheless, Marin & Antonucci (2016) argued that the probability to see a type-1 source at high inclination is almost zero. Based on the polarisation properties of AGNs in the 3CRR catalogue the authors concluded that the optical classification is a reliable indicator for the BLR orientation.

A less drastic solution for the mismatch between the inclination and Mrk 1044’s unobscured nucleus is to discard our toy model’s rigorous assumptions on the geometry multi-phase outflow. Since we have not spatially resolved any of the components, the Ly-α absorbing clouds may not be linked to the [O III] emitting clouds. The different kinematics between the Ly-α absorbers and [O III] emitters could thus be a reflection of the inhomogeneity of the ISM on sub-pc scales regarding the ionisation conditions and the density distribution. Compared to the [O III] ionised gas, the Ly-α absorbers have higher velocity but smaller dispersion. This can be understood if the outflow is launched from accretion disc scales and entrains the ambient gas, thus gradually decreasing its bulk velocity while increasing the dispersion/turbulence due to the interactions with the surrounding ISM (Gaspari & Sądowski 2017). This can also potentially trigger a CCA condensation rain (see Sect. 4.5). In an alternative interpretation, the outflow is intrinsically elongated but aligned with the line-of-sight. Due to the optically thin [O III] clouds, a velocity gradient within the medium would increase the integrated observed line width and mimic a high dispersion of the spatially unresolved structure. Both explanations favour that the [O III] emission occurs at larger distance from the nucleus than the Ly-α absorption. A qualitative illustration of this alignment is shown in Fig. 7 where the fast Ly-α absorbing clouds are located closer to the nucleus. We note that the inclination and geometry of the [O III] clouds remain unconstrained. For the following discussion, we therefore do not adopt the specific value that predicted by the expanding shell model.

Fig. 7. Qualitative illustration of a possible alignment of the unresolved outflowing systems. Combining the dust torus model with the expanding shell model, the Ly-α absorption cannot be spatially linked with the [O III] emitting cloud. Hence, the inclination of the system as well as the geometry of the [O III] emitting structure is unconstrained. Compared to the extended [O III] emitting region, the Ly-α clouds are likely to be located closer to the AGN engine due to their higher velocity.

4.3. Outflow energetics and mass outflow rates

Despite the compactness of the ionised gas outflows at present day, it will propagate through the ISM of the host galaxy. Once the shock front reaches the molecular gas from which the CNE produces stars, the outflow’s energy and momentum will induce turbulence into the ISM (e.g. Wittor & Gaspari 2020). To estimate the future impact on the host galaxy star formation, we compute the current mass, energy and momentum ejection rates released through the ionised gas outflow. We stress that the following derivations require an assumption about the morphology of the [O III] emitting clouds and their velocity as derived in Sect. 4.2. As discussed in the previous section, the Ly-α absorbers may not be spatially associated with the ionised gas outflow. We therefore use the velocities of the [O III]-emitters for the following derivations.

As a first step, we derive the total Hα-flux of the outflowing components O1 and O2 by integrating over the individual components’ flux maps extracted in Sect. 3.2. From the Balmer-decrement Hα/Hβ and the Milky Way-like attenuation curve Cardelli et al. (1989), we derive the attenuation analogously to Paper I. As described in Husemann et al. (2016a), the ionised gas mass can then be estimated from the extinction-corrected intrinsic luminosity L_Hα as

$$\begin{array}{c}\hfill M_{\mathrm{ion}}=3.4\times {10}^6\left(\frac{100{\mathrm{cm}}^{-3}}{n_{\mathrm{e}}}\right)\left(\frac{L_{\mathrm{H}\alpha }}{{10}^{41}{\mathrm{erg}\mathrm{s}}^{-1}}\right)M_{\odot }.\end{array}$$(1)

Due to the missing [S II] and [N II] emission (see Sect. 4.1), we assume a lower limit on the electron density n_e  =  8.7  ×  10⁴ cm⁻³, the critical density of [N II]. In this way we ensure that the outflow masses and energies are not underestimated.

The estimated mass outflow rate depends on the assumed outflow morphology (Cicone et al. 2014; Veilleux et al. 2017). Since the UV absorbers have relatively low velocity dispersion, we consider thin shells to be the most probable description of the [O III] outflow morphology. However, in order to achieve a consistent comparison between our results and the scaling relations from Fiore et al. (2017), we assume a conical outflow geometry which lowers the derived values by a factor of approximately three. We compute the ionised gas mass outflow rate as

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{OF},\mathrm{cone}}=3\times \left(\frac{v_{\mathrm{OF}}}{100\mathrm{km}{\mathrm{s}}^{-1}}\right)\left(\frac{M_{\mathrm{ion}}}{{10}^7{M}_{\odot }}\right)\left(\frac{1\mathrm{kpc}}{R_{\mathrm{OF}}}\right)M_{\odot }{\mathrm{yr}}^{-1}\end{array}$$(2)

where v_OF is the outflow velocity and R_OF is the size of the homogeneously filled cone.

In addition to the mass outflow rate, we estimate the momentum injection rate ${\dot{p}}_{\mathrm{OF}}=v_{\mathrm{OF}}{\dot{M}}_{\mathrm{OF}}$ and the kinetic energy injection rate as ${\dot{E}}_{\mathrm{kin}}=0.5v_{\mathrm{OF}}^2{\dot{M}}_{\mathrm{OF}}$. Finally, we compare the SFR on different spatial scales and the mass injection rate from the AGN-driven outflow by computing the mass loading factor as η = Ṁ_OF/SFR.

The results for each of the outflowing components are listed in Table 4. The energetics of Mrk 1044’s outflow properties are affected by several systematic uncertainties. Especially the electron density n_e and the size of the system might underestimate each of the derived quantities, as we could only estimate lower and higher boundaries respectively. Nonetheless, the derived outflow properties hold some important implications especially for the future impact on the host galaxy discussed in Sect. 4.5. The two unresolved [O III] are luminous and dominate the outflow mass and energy budget by four orders of magnitude. The total kinetic energy of the outflows is 7.3 × 10⁴¹ erg s⁻¹ which corresponds to 0.2% of Mrk 1044’s bolometric luminosity (3.4 × 10⁴⁴ erg s⁻¹, Husemann et al. 2022). This favours the possibility of an energy-driven ionised gas outflow, if the photons couple with a comparable efficiency to the gas in the host galaxy ISM.

4.4. Origin of Mrk 1044’s radio emission

Due to the radio-quite nature Mrk 1044 and the lack of extended radio jets, its nuclear radio emission could be produced by a multitude of processes. This includes the accretion disc coronal activity, an AGN driven wind, free-free emission from photo-ionised gas, or nuclear star formation. The radio spectra of H II regions is shaped by the free-free emission which, in the optically thin limit, produces a nearly flat radio spectrum with α  ∼   − 0.1. Together with the steep synchrotron spectrum from supernova remnants, the integrated radio emission of H II regions is expected to have a characteristic spectral index of ∼ − 0.7 which is consistent with the observations of star-forming galaxies (Condon 1992; Panessa et al. 2019; Pérez-Torres et al. 2021). To test whether the unresolved nuclear SFR is sufficient to produce the Mrk 1044’s nuclear radio emission, we estimate the SFR from the observed radio continuum luminosity and the calibration from Murphy et al. (2011):

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{\mathrm{SFR}}{M_{\odot }{\mathrm{yr}}^{-1}}=& {[2.18{\left(\frac{T}{{10}^4\mathrm{K}}\right)}^{0.45}{\left(\frac{\nu }{\mathrm{GHz}}\right)}^{-0.1}+15.1{\left(\frac{\nu }{\mathrm{GHz}}\right)}^{{\alpha }^{\mathrm{NT}}}]}^{-1}\hfill \\ \hfill & \times \frac{L_{\nu }}{{10}^{27}{\mathrm{erg}\mathrm{s}}^{-1}{\mathrm{Hz}}^{-1}}·\hfill \end{array}\end{array}$$(3)

Here, α^NT corresponds to the non-thermal spectral index for which we use the value α_peak = −0.61 ± 0.16 measured between 5 GHz and 10 GHz (see Sect. 3.4). Using Mrk 1044’s C-band luminosity and an electron temperature of 10⁴ K we use this relation to estimate a star formation rate of SFR = (0.81 ± 0.45) × 10⁻³ M_⊙ yr⁻¹. The associated uncertainty involves both the error of our measured parameters (see Sect. 3.4) together with the 50% uncertainty from relation (3). Since the SFR estimated from the radio luminosity is exceeded by the 0$\stackrel{″}{.}$5 unresolved SFR derived from Mrk 1044’s Hα-luminosity by one order of magnitude (see Sect. 4.1.1), we cannot exclude that the radio emission is partially produced by SF-related processes. Thus, we conclude that neither the spectral index nor the amount of flux allows us to distinguish whether Mrk 1044’s radio emission originates exclusively from AGN-related processes or the nuclear star formation.

4.5. Outflow expansion from nuclear to galaxy scales

4.5.1. Kinematic signature of the outflow

To test whether the individual components of the multi-phase outflow have the potential to reach galactic scales, we use the ‘k-plot’ diagnostic from Gaspari et al. (2018). The k-plot compares the line-of-sight (LOS) velocity dispersion with the LOS velocity shift. It is useful to assess the role of chaotic and turbulent motions in addition to laminar and bulk motions, and related key physical processes such as chaotic cold accretion (CCA) and AGN outflows. We use the pencil-beam approach, which better traces individual clouds or small-scale gas elements (e.g. Rose et al. 2019; Maccagni et al. 2021; North et al. 2021; Temi et al. 2022).

Figure 8 shows how the CNE gas (black points) and nuclear components (red/green) of Mrk 1044 compare with the multi-phase k-plot retrieved via high-resolution hydrodynamical simulations (yellow 1-3σ contours; see Gaspari et al. 2018) and related observational sample of diverse galaxies (grey). The ionised gas in Mrk 1044 star-forming CNE with the highest velocity shows a symmetric spread of the dispersion within the 2σ contours. The tail towards low velocities at constant dispersion is more heavily associated with the rotating disc. Despite the dominance of rotation, several black points reside well within the 1-σ CCA region, suggesting that such gas elements can also be prone to turbulent condensation. This is a common scenario in spiral galaxies where the turbulent Taylor number is on average above unity, hence favouring a rotating CCA (Gaspari et al. 2015).

Fig. 8. Kinematical plot (k-plot) of the line broadening versus the line shift. It highlights the relative importance of chaotic/turbulent motions versus laminar/bulk motions and the connected physical processes. The yellow contours show the 1-3σ confidence intervals for the chaotic cold accretion (CCA) feeding mode found by high-resolution hydrodynamical simulations (Gaspari et al. 2018). The grey data points indicate measurements of warm and cold condensed gas in diverse galaxies (Gaspari et al. 2018). The ionised gas of Mrk 1044’s CNE is shown as black data points. The tail towards lower velocity is mainly driven by an inclined rotating disc, with several points overlapping with the main CCA area. In contrast, the four nuclear components (red and green points) have substantially higher bulk velocity indicating that they are not kinematically associated with the disc, but rather outflowing, in particular O2 and L2 (squares). O1 and L1 appear to interact more significantly with the surrounding medium, thus being more prone to turbulent condensation.

Compared to the region the resolved data points of the rotating CNE occupy, the unresolved nuclear components separate. The faster components – O2 and L2 – are located in the 3σ area and are thus expected to be less affected by condensation and CCA. This suggests that their momentum is currently little dissipated by the interaction with the surrounding medium, which enables the outflows to escape at least from the high-density nuclear region. Due to their lower velocities, systems O1 and L1 are located within the 1.5σ yellow confidence region. This can be interpreted that such outflowing components have started to deposit their kinetic energy into the ISM, increasing the turbulence of the ambient medium and potentially stimulating an initial CCA rain.

It is interesting to note that the [O III] outflows have systematically higher velocity dispersion than the rest of the system. This can be interpreted in a few ways. On the one hand, the highly ionised outflow component could simply have a stronger intrinsic turbulence along its spine, which is indicative of the entrainment of surrounding material that increases the turbulence in the ionised gas. On the other hand, it could be interpreted as a projection effect if the [O III] gas velocity is directed at a small inclination with respect to the line-of sight. In this case the [O III] emission could be affected by beam smearing as opposed to the kinematics from Ly-α foreground absorption. The observed [O III] line width might therefore reflect a systemic velocity gradient within the outflow and might not be interpreted as a dispersion.

4.5.2. Future impact on the host galaxy

Mrk 1044’s outflow carries a lot of mass in the ionised gas phase. To estimate its future impact on the host galaxy, we compare the summed energetics with the scaling relations of resolved galaxy-scale outflows. The left panel of Fig. 9 shows that compared to other AGNs at similar luminosity, Mrk 1044’s ionised gas outflow lies within the scatter of the outflowing mass–AGN bolometric luminosity correlation from Fiore et al. (2017). It is more massive by ∼1 dex than average, whereas the kinetic power is consistent with the ionised gas outflows in lower luminosity AGNs. To compare the ionised gas outflow momentum rate with the AGN radiation momentum rate L_bol/c, we compute the wind momentum load as ṧ_OF/ṧ_AGN = 3.3. It is consistent with the wide spread of the extended ionised gas outflows reported in Fiore et al. (2017) and close to the expectation of ṧ_OF/ṧ_AGN = 1 for a momentum conserving outflow predicted by the King (2003) model. Figure 9 highlights that multi-phase outflows often carry a large amount of their mass and kinetic power in the molecular gas phase, especially at low AGN luminosities (see also Fluetsch et al. 2019; Veilleux et al. 2020). We therefore expect that the AGN-driven mass injection into the host galaxy ISM is even higher if Mrk 1044’s outflow is accompanied by a molecular gas outflow.

Fig. 9. Comparing Mrk 1044’s integrated ionised gas outflow with AGN multi-phase outflow scaling relations. The coloured data points are taken from Fiore et al. (2017) and show the molecular gas outflows traced by CO, OH (blue), ionised gas outflows (green) and X-ray outflows (red). Panel (a) shows the outflowing ionised gas mass. The dashed coloured lines correspond to the correlations for the respective gas phase reported by Fiore et al. (2017). Panel (b) shows the AGN wind kinetic power where the dashed, dotted and dash-dotted line correspond to the ${\dot{E}}_{\mathrm{kin}}=0.01,0.1,1L_{\mathrm{bol}}$. Panel (c) shows the mass loading factor η = ṀOF/SFR. Here, the brown star represents Mrk 1044’s SF in the immediate vicinity (< 160 pc) from the centre, the orange point to that of the SF in the CNE and the yellow point to that of Mrk 1044’s galaxy-wide SFR. Panel (d) shows the kinetic outflow power as a function of SFR. The solid red line indicates the predicted SFR from the AGN bolometric luminosity using the relation from Netzer (2009) together with the average ${\dot{E}}_{\mathrm{kin}}/L_{\mathrm{bol}}=0.025$ for molecular outflows (Fiore et al. 2017). Compared to the ionised gas outflows in the sample, Mrk 1044’s outflow carries a lot of mass in the ionised gas phase. Its mass injection into the ISM exceeds the current SFR in the immediate vicinity of the nucleus, but also compared to the galaxy-wide values.

An important difference when comparing Mrk 1044 with ionised gas outflows in nearby galaxies is that they extend over different spatial scales within the host galaxy. While the relation from Fiore et al. (2017) was calibrated for ionised gas winds that extend over > 1 kpc, Mrk 1044’s ionised gas outflow is compact (< 4.6 pc) and located close to the galaxy nucleus (< 2.1 pc). Its proximity to Mrk 1044’s central engine is reflected in the mass loading factor η = Ṁ_OF/SFR which relates the SFR with the mass outflow rate. The intrinsic extent over which the galaxy property that is extracted, in this case the SFR, does affect the result dramatically, as it is shown in the third panel of Fig. 9. While the outflow is located at parsec-scale distance from the nucleus, the SFR can be measured in different apertures. If we compare the mass outflow rate with the SF in its immediate vicinity, that is the SF ‘hidden’ by the beam-smeared outflow (< 160 pc, see Sect. 4.1), ${\eta }_{0\stackrel{″}{.}5}=276$ suggests that the outflow currently dominates the ISM properties of the host galaxy on a few ∼100 pc scales. Such a high mass-loading factor is in line with the picture of a young outflow that has recently been launched by the powerful AGN. Furthermore, the ionised gas outflowing mass exceeds Mrk 1044’s present day SFR, even if we compare it with the galaxy-wide SFR. However, this comparison should be regarded with caution as both quantities trace different regions of the host galaxy and the evolution of the outflow’s properties during its propagation through the ISM are largely unknown.

Since Mrk 1044’s ionised gas outflow inclination and the morphology are unconstrained, it is not clear how the outflow will impact the host galaxy. To estimate the maximum future impact on the host galaxy, we assume that the outflow is launched at large inclination relative to the galaxy rotation axis and carries its present-day mass and energy out to a few 100 pc from the nucleus. In this case it directly impacts the CNE and the galaxy disc where the host galaxy star formation is concentrated. The total mass outflow rate also exceeds the integrated star formation rate in the CNE (0.19 ± 0.05 M_⊙ yr⁻¹) and the galaxy-wide SFR (0.70 ± 0.17 M_⊙ yr⁻¹) by more than one order of magnitude. Mrk 1044’s multi-phase outflow may therefore have the potential to deprive the host galaxy from their molecular gas reservoir and eventually quench the host galaxy star formation. Based on the mass loading factor, Fiore et al. (2017) argued that only powerful AGNs are able to drive the co-evolution between SMBHs and galaxies. In this interpretation Mrk 1044 would represent an uncommon example where the host galaxy evolution is regulated by the central SMBH despite the relatively low luminosity of its AGN if compared to the overall AGN population.

In the other extreme scenario the outflow is orientated perpendicular to the disc, the cold gas by which the SF in the CNE is fuelled will not be affected by the energy injection. Hence, we expect that in this case the outflow’s impact on the host galaxy is minor to none. The two extreme cases of outflow inclinations leave a wide range of possible scenarios in which the host galaxy SF may be quenched by the injection of turbulence into the ISM. Resolving the inclination of the accretion disc relative to the galaxy plane is therefore crucial to estimate the future impact of Mrk 1044’s outflow on the host galaxy.

4.6. Mrk 1044 – A young AGN

In Paper I we have touched on the idea that Mrk 1044 as a NLS1 could be archetypal for a population of ‘young’ AGN that host a rapidly growing SMBH. There we have shown that the circumnuclear SF enriches the ISM and shows potential signatures of ongoing BH feeding. SF is present even in the immediate vicinity of Mrk 1044’s nucleus (see Sect. 4.1), and might play a role in channelling the gas towards the centre as suggested by hydrodynamic simulations (Davies et al. 2007; Volonteri et al. 2015; Zhuang & Ho 2020). Since circumnuclear SF continues to enrich the ISM and has not yet been affected by the energy injection from the expanding multi-phase outflow, this process might continue until the outflow propagates through the ISM.

To estimate a frame for the look back time at which of the ionised gas outflow was launched, we assume a constant LOS velocity by which the outflow escapes from the central engine. Since the outflow inclination cannot be robustly constrained by the geometric toy model discussed in Sect. 4.2, we assume that the outflow is launched perpendicular to the accretion disc. Further, we assume that the rotation axis of the accretion disc is aligned with that of the dust torus. In this case, the accretion disc inclination i of a type 1 AGN cannot be larger than the half-opening angle of the dust torus, since the torus would obscure the view to the BLR. We therefore assume Mrk 1044’s outflow inclination to be smaller than 58°, the median of the dust torus half-opening angle distribution (Ricci et al. 2017). From the intrinsic velocity v = v_LOS/cosi and the projected distance from the nucleus d = d_proj/sini we estimate the timescale as

$$\begin{array}{c}\hfill t⪅\frac{d_{\mathrm{proj}}}{v_{\mathrm{LOS}}}tani.\end{array}$$(4)

With the maximum projected offset from the nucleus of 0.85 pc (see Sect. 3.2) we infer that the two ionised gas outflow components were launched no longer than 9900 yr (O1) and 3500 yr (O2) ago. Using the upper limit for the location of the spatially resolved outflow H0, we estimate an upper limit for its age of ∼25 000 yr. The derived ‘age’ estimates should be regarded as an upper limit since interaction of the fast-moving ionised gas outflow with the host galaxy ISM would slow down the outflow during its expansion. Furthermore, if the outflow velocity is directed closer to the observer (i.e. at smaller inclination i), the resulting timescale would be shorter. The picture of Mrk 1044 as a young AGN is consistent with the undisturbed kinematics that we could trace down to $0\stackrel{″}{.}5$ (∼160 pc) from the centre (see Paper I). Furthermore, the ongoing star formation in the CNE around Mrk 1044’s nucleus (∼300 pc) has not been quenched by the energy injection from the multi-phase outflow. Even within the innermost 160 pc where the ionised gas outflow dominates, the integrated line ratios that belong to the host are consistent with ionisation through SF (see Sect. 4.1).

Since the distance of an AGN-ionised gas cloud to the AGN engine directly translates into a light travel time, the ENLR size can be interpreted as a proxy for the lifetime of a single AGN episode (Lintott et al. 2009; Keel et al. 2017; Husemann et al. 2022). At fixed luminosity, AGNs with more massive BHs would be statistically observed at much later times in their episodic phase than lighter BHs. This interpretation implies an evolutionary sequence for AGNs where systems hosting less extended ENLRs have a shorter AGN lifetime or are ‘younger’. The EV1 relation could be explained as a reflection of the time-dependence at which the AGN phase is observed. NLS1s, which are located in the tail of the EV1 relation, would represent the extreme end of the rest-frame difference between [O III]-core and systemtic redshift. The interpretation of NLS1s as ‘young’ AGN would explain the relative weakness of their narrow [O III] emission lines by a combination of (i) the absence of an AGN-ionised ENLR in the host galaxy (ii) the compactness of the ionised gas outflow and (iii) the high gas densities in the compact outflow (see Sect. 4.1). In this scenario, the prevalence of strong broad Fe II emission in the spectra of NLS1s which can be interpreted as the finite timescale on which AGN quenching is able to shut off the host galaxy star formation on nuclear scales.

Another aspect of the interpretation is that it explains Mrk 1044’s offset of the narrow [O III]-core from the galaxy systemic velocity. As we have demonstrated, the [O III]-core component is outflowing with −140 km s⁻¹ (Sect. 3.1) which is fundamentally different from many AGNs where the narrow [O III] traces the ENLR associated the galaxy rest-frame (e.g. Harrison et al. 2016; Husemann et al. 2022). If a young AGN drives the expanding outflow, this mismatch can be understood since the [O III]-core component traces a young structure that has a smaller intrinsic extent. Hence, the integrated narrow [O III] may trace different structures, depending on the age of the system. This implies that the outflow velocities of ‘young’ AGN cannot be estimated from the integrated emission line spectrum alone. Spatially resolving young AGN host galaxies and their compact outflows is crucial to correctly interpret their spatial and kinematic structure.

A coherent picture for the recent activity of Mrk 1044’s BH growth emerges. Combining the arguments of (i) a compact outflow, (ii) unperturbed host kinematics and (iii) SF-dominated radiation field in the close to the nucleus, we conclude that Mrk 1044’s AGN phase has been triggered recently. Following the idea that NLS1s host rapidly growing SMBHs at an early stage of their evolution, we predict that highly accreting AGN host galaxies with low-mass SMBHs show similar signatures close to the nucleus.

5. Summary and conclusions

In this work we have combined IFU data from MUSE NFM-AO with the UV spectrum from HST to constrain the properties of Mrk 1044’s outflow. We have detected two outflowing systems in both ionised gas emission and Ly-α absorption that are unresolved at the spatial resolution of MUSE NFM-AO. We have used a spectroastrometric approach to locate the ionised gas outflows traced by [O III] emission and estimated their energetics. Our key results are the following.

• Mrk 1044’s narrow core component O1 in which we detect Hβ, Hα, and [O III] emission is outflowing with v_O1  =   − 144  ±  5 km s⁻¹. The [O III]-wing component O2 represents an additional component with v_O2  =   − 560  ±  20 km s⁻¹.

• Both ionised gas outflows O1 and O2 are spatially unresolved and have a projected offset from the AGN location of smaller than 1 pc. These two outflows carry the majority of the ionised gas mass and energy and they have been launched no longer than ∼10⁴ yr ago.

• Despite the compactness of Mrk 1044’s ionised gas outflow, its mass outflow rate and kinetic energy injection rate are comparable to typical AGN-driven ionised gas winds at a similar AGN luminosity. The high present-day outflow velocity and high mass loading suggest that the outflow has the potential to expand to galaxy scales and to impact the ISM properties.

• Using a geometric toy model, we could not identify Ly-α absorbers with the unresolved [O III] via their kinematics. This emphasises the complex nature of the multi-phase ISM on sub-parsec scales.

• The gas densities in Mrk 1044’s spatially resolved outflow H0 are above the critical density of the forbidden [S II] and [N II] line doublets, which implies that the classical BPT diagnostics cannot be used to constrain the underlying ionisation mechanism.

Our results demonstrate the importance of studying the AGN–host galaxy connection on different energy and density scales with spatially resolved observations. Although the ionised gas outflow can be identified from the AGN spectrum alone, Mrk 1044’s compact multi-phase outflow stands out from well-studied examples of kiloparsec-scale extended outflows. The different outflow components, kinematics, and energetics highlight the complexity of the AGN feedback cycle, especially on the smallest galactic scales from where the outflow is launched. To better understand the role of the warm ionised gas in the immediate vicinity of accretion-mode AGNs, more IFU observations with high spatial resolution are required. Performing a similar analysis on a larger sample will help to constrain the impact of young, compact outflows on the host galaxy as well as the underlying physical mechanisms behind AGN feeding and feedback.

Acknowledgments

We thank the anonymous referee for suggestions and comments that helped improve the presentation of this work. This work was funded by the DFG under grant HU 1777/3-1. N.W. and B.H. acknowledge travel support from the DAAD under grant 57509925 and the hospitality of MS at the University of Manitoba where large parts of this paper were written. B.H. was partially supported by the DFG grant GE625/17-1. V.N.B. gratefully acknowledges assistance from the National Science Foundation (NSF) Research at Undergraduate Institutions (RUI) grant AST-1909297. It is important to note that findings and conclusions do not necessarily represent views of the NSF. C.O. and M.S. are supported by the Natural Sciences and Engineering Research Council (NSERC) of Canada. M.G. acknowledges partial support by NASA Chandra GO9-20114X and HST GO-15890.020/023-A, and the BlackHoleWeather program. M.P.T. acknowledges financial support from the State Agency for Research of the Spanish MCIU through the “Center of Excellence Severo Ochoa” award to the Instituto de Astrofísica de Andalucía (SEV-2017-0709) and through the grant PID2020-117404GB-C21 funded by MCI/AEI/10.13039/501100011033, by ERDF “A way of making Europe”.

References

Appendix A: Size of the resolved outflow H0

The system H0 is one of the identified ionised gas outflows’ components. It is present in the AGN-subtracted host data cube and therefore spatially resolved. In Fig. A.1 we zoom in onto the central $0\stackrel{″}{.}1$ of the narrow Hβ surface brightness map where the outflow flux, velocity and velocity dispersion peak. To constrain the outflow’s location and projected size, we fit the surface brightness profile with a model light distribution of an intrinsically extended source. Our model consists of (i) the empirical PSF extracted at Hβ which we convolved with a Gaussian kernel r_OF, and (ii) a constant background flux b. To find the best-fitting parameters, we minimise the S/N-weighted χ² sum of the Hβ residual flux distribution.

Fig. A.1. Measuring the size of Mrk 1044’s spatially extended Hβ outflow. Panel (a) shows the flux map of the Hβ outflowing component which we extracted as described in Paper I. Panel (b) shows the best-fitting model that consists of the empirical PSF convolved with a Gaussian kernel and a constant background flux. The residuals in panel (c) show that broadened PSF describes the Hβ profile well. Panel (d) shows the χ2 distribution in the parameter plane where a global minimum is present around 0.55 px (4.6 pc).

The best-fit model is shown in the top right panel of Fig. A.1. It yields an intrinsic size of r_OF = 0.55 ± 0.10 px (4.6 ± 0.6 pc) from which we conclude that the source is spatially resolved at a > 3σ confidence level. This extended component of the outflow is also present in Hα emission. Here, we follow the same method and retrieve an intrinsic projected extent of 0.38 ± 0.02 px (3.2 ± 0.1 pc). Both values are consistent within the uncertainties, which suggest that both light profiles stem from the same structure. We measure the centroid of the light distribution which matches the AGN position with an uncertainty of 0.22 px (1.8 pc). Our results are robust against choosing different spectral regions of the broad line windows for the empirical PSF extraction that is described in Paper I. Furthermore, our results do not change within the uncertainty if we use the AGN-subtracted host emission line maps presented in Paper I instead of the Hβ and Hα surface brightness maps that we retrieve from the spectroastrometric analysis. The former were generated from fitting the AGN-subtracted host emission with a single-Gaussian component, which is independent but fundamentally different method compared to the spectroastrometric measurement.

Appendix B: Geometric expanding shell model

Here, we describe how we use a geometric model of two expanding shells to reproduce the [O III] emission line shape. With this model we can constrain the geometric alignment of the unresolved ionised gas outflow relative to the Ly-α absorbers. The implication of the results are discussed in Sect. 4.2.

B.1. Constraining the outflow geometry

In contrast to the unresolved but intrinsically extended [O III] outflow, the UV spectrum only probes the line-of-sight absorption towards the high density BLR clouds. We can link the emitters and absorbers by assuming that a small fraction of the [O III]-emitting gas is also responsible for the Ly-α absorption. In this scenario, the systematic velocity factor of ∼2 between [O III]-emitters and Ly-α-absorbers may be the result of a projection effect.

There exist different descriptions for outflow morphologies including uniformly filled cones (e.g. Cano-Díaz et al. 2012; Cresci et al. 2015; Fiore et al. 2017), hollow bi-cones (e.g. Müller-Sánchez et al. 2011) or thin shell structure (e.g. Rupke et al. 2005). In our case, the relatively narrow Ly-α absorption features have low velocity dispersion, which suggests that the absorbing structure is better described by a shell-like shock front scenario. We therefore assume that the ionised gas is distributed on two concentric shell caps with the AGN located in the origin. The shell caps expand radially with a velocity v. We further assume that the receding side is obscured by dust. Both shell caps have a half-opening angle θ₀ that is observationally constrained. Ricci et al. (2017) applied a torus model to a sample of Swift BAT AGN for which they constrained the median value of the torus half-opening angles to 58 ± 3°. We adopt this value for our model, although varying θ₀ by ±10° does not alter our results qualitatively. We assume that luminosity profile on the shell cap is centrally concentrated. We describe it by an exponential luminosity profile that declines radially with a half-light radius α in units of the half-opening angle θ₀. Thus, α/θ₀ = 1 means that the half-light radius coincides with the half-opening angle. Furthermore, we include the luminosity ratio between the shells as a free parameter. An illustration of the geometric model is shown in Fig. B.2. The two shells share the geometric parameters i and α, but not for the expansion velocity v that is fixed to the measured velocity of the Ly-α absorbers.

Fig. B.1. Illustration of the expanding shell geometry. The Ly-α occurs in the expanding shell, along the LOS towards the BLR with the velocity v. The [O IIII] emission originates from an intrinsically extended but spatially unresolved shell cap with a half-opening angle θ0, half-light radius α and inclination i. The observer measures the integrated flux and thus measures a smaller velocity.

B.2. Predicting the emission line shape

We use a forward-modelling approach to constrain the parameters of the expanding shell model from the [O III]λλ4959,5007 emission lines. We predict the observed emission line shape with a numerical approach that involves the following steps.

(1) As a first step, we sample the shell cap with 10⁵ equidistant cells using the Fibonacci sphere algorithm⁴.

(2) We incline the shell cap using a rotational coordinate transform.

(3) For each of the cells we compute the line-of-sight velocity component and performed a binning in velocity space. This provides us with the predicted emission line shape.

(4) To fit the rest-frame spectrum of Mrk 1044 we transform the distribution into the wavelength space using the expansion velocity of the shell together with the rest-frame wavelengths of the respective ion. With this method we generate a template of emission lines shapes, computed for a 10² × 10² parameter grid with i = [0, 90) and α = (0, 1].

Since there are two emitters and absorbers respectively, our model requires a fast and a slow component. To illustrate the effect of the geometric parameters in Fig. B.2, we show the dependence of the fast component’s peak emission on the inclination of the shell caps. Higher inclinations move the peak of the line closer to the systemic velocity v_sys = 0 km s⁻¹, whereas the effect of the central concentration α/Ω depends on the inclination of the system. The solution for the line shape for the slow component behaves equivalently.

To find the linear combination of the fast and slow component that fits the [O III] best, we use a non-negative least square optimisation algorithm. Since we only see the approaching side of the outflow, we only model the blue shoulder of the [O III]λλ4959,5007 in the with corresponding velocity range [v_L2, 0 km s⁻¹]. To get an estimate of the parameters’ uncertainties, we use a Monte Carlo approach. From the original AGN spectrum we generate 10⁴ mock spectra by randomly varying the flux density within the uncertainty. From the resulting posterior distribution, we adopt half of the 16th to 84th percentile range for the parameters’ uncertainties.

B.3. The best-fitting geometric model

Using the concentrically expanding shell model, we retrieve an inclination of i = 79.1 ± 8.4°. The best-fit concentration α/θ₀ = 0.09 ± 0.04 indicates a relatively narrow (‘jetted’) luminosity profile. The best-fit geometric parameters have an intuitive explanation. The width of the [O III] primarily determines the concentration of the luminosity profile. Since we do not observe a double-peaked emission line, both emission lines have to be close in wavelength. This can only be achieved at high inclinations, as shown in Fig. B.3.

Fig. B.2. Dependence of the emission line flux maximum corresponding velocity on inclination. Independent of how concentrated the light profile is (α/θ0), the geometric model requires high inclinations of the [O III]-emitting shell to significantly reduce the peak velocity by a factor of approximately two.

Fig. B.3. Best-fitting model for the blue shoulder of Mrk 1044’s [O III]λλ4959,5007 emission line. The left panel shows the AGN spectrum (black line) around the [O III] doublet emission. The best-fit model (orange line) reproduces the line shape well. Both the slow (red line) and the fast shell (blue line) contribute to the asymmetric line shape. The panels on the right show percentile regions of the posterior distribution for the geometric parameters (left) and the flux contribution (right) respectively.

Our results do not change within the uncertainties if the [O III] lines are fitted independently. Moreover, our results do not depend on how we define the ‘blue shoulder’, i.e. slightly adjusting the wavelength-borders around the emission line does not affect the results significantly. Introducing another geometric parameter leads to an over-fitting of the [O III] emission lines given their relatively low spectral S/N. Furthermore, we have tested different parameterisations to describe the emission line profile. In particular, different geometric descriptions like that of a hollow cone (e.g. Das et al. 2006; Müller-Sánchez et al. 2011; Venturi et al. 2018) require similar inclinations to explain the [O III] line shape. We conclude that the absence of a single-peaked [O III] emission line can only be reproduced by a concentrically expanding shell model if the system is highly inclined with respect to the LOS.

All Tables

All Figures

Fig. 1. Modelling Mrk 1044’s AGN spectrum in the Hβ–[O III] region (left) and Hα–[N II] region (right). To reproduce the observed spectrum (black) we use a multi-Gaussian model for Hβ (blue), [O III] (green), Fe II (purple), Hα (blue, right panel) and [N II] (red). The best-fit spectrum is shown as a red line and well reproduces the prominent emission lines, including the blue shoulder of the [O III] narrow line.

In the text

Fig. 2. Spectroastrometric analysis of Mrk 1044’s central region using the original MUSE NFM-AO data cube. The left panels show an arbitrary example spectrum that is picked from the small window around the nucleus shown in the panels on the right. To fit the spectrum, we keep the kinematics and line ratios of the kinematics fixed to that of the components found in the AGN spectrum. The panels on the right show the surface brightness within the central 150 mas for the kinematic components where the flux maps are normalised to their peak flux. Here we only show the measured light distribution from the BLR Hβ (left), [O III]-wing (middle), and narrow [O III]-core (right) component. From top to bottom the maps show the 2D light distributions, the best-fit model of the PSF and the residual maps. The red cross indicates the PSFAO19 centroid for the best-fit PSFAO19 fit to the BLR Hβ. The green and white crosses indicate the centroids for the PSFAO19 models to the surface brightness maps of the [O III]-wing and core respectively, in which the amplitude and position of the Moffat model were varied only. The offset of both [O III] core and wing component from the AGN position are smaller than 0.1 px.

In the text

Fig. 3. Origin of the forbidden [S II] emission in Mrk 1044’s centre. The left panel shows the Hα surface brightness map where the H0 outflow is located the central 0$\stackrel{″}{.}$5. The right panel shows the two component model for the inner 0$\stackrel{″}{.}$5 integrated spectrum where we fixed the host and outflow kinematics to what we found from the 3″ integrated spectrum. For the outflow, we do not find a significant contribution of [S II] and [N II] suggesting that the gas density lies above their critical densities.

In the text

Fig. 4. Mrk 1044’s nuclear region seen in the optical (left) and radio emission (right). Panel (a) shows the surface brightness map of the Hα narrow line emission extracted as described in Paper I. The panels on the right show the broad-band EVLA images of the radio emission in the C-band (b) and X-band (c) with the corresponding FWHM-contour of the beam shown in the lower left corner respectively. The contours of the star forming CNE encompassing the nucleus are shown as black contours (see Paper I). In both C- and X-band the nuclear radio emission is unresolved by VLA.

In the text

Fig. 5. Modelling Mrk 1044’s UV spectrum in the Ly-α window. We use a multi-Gaussian model for the Ly-α emission (blue), and the absorbing systems (orange). The best-fit model (red line) reproduces the spectrum (black line) well in the selected wavelength range. Since the Ly-α systems L1 and L2 are accompanied by high metallicity absorption lines at the same velocity (Fields et al. 2005), the absorbing clouds must be local to the galaxy.

In the text

Fig. 6. BPT diagnostics of the extended outflow H0. The grey data points correspond to the 8 × 8 co-added spaxels that belong to the star-forming CNE (see Paper I). The stars represent the two components that are present in the $0\stackrel{″}{.}5$-integrated spectrum that is shown in Fig. 3. The host component (red star) is consistent with SF-classification of the line ratios recovered in the CNE. For the outflowing component (blue star), both [O III] and [N II] emission might be suppressed due to collisional de-excitation in the high density gas. The BPT classification of the outflow as SF-ionised might therefore be misleading.

In the text

Fig. 7. Qualitative illustration of a possible alignment of the unresolved outflowing systems. Combining the dust torus model with the expanding shell model, the Ly-α absorption cannot be spatially linked with the [O III] emitting cloud. Hence, the inclination of the system as well as the geometry of the [O III] emitting structure is unconstrained. Compared to the extended [O III] emitting region, the Ly-α clouds are likely to be located closer to the AGN engine due to their higher velocity.

In the text

Fig. 8. Kinematical plot (k-plot) of the line broadening versus the line shift. It highlights the relative importance of chaotic/turbulent motions versus laminar/bulk motions and the connected physical processes. The yellow contours show the 1-3σ confidence intervals for the chaotic cold accretion (CCA) feeding mode found by high-resolution hydrodynamical simulations (Gaspari et al. 2018). The grey data points indicate measurements of warm and cold condensed gas in diverse galaxies (Gaspari et al. 2018). The ionised gas of Mrk 1044’s CNE is shown as black data points. The tail towards lower velocity is mainly driven by an inclined rotating disc, with several points overlapping with the main CCA area. In contrast, the four nuclear components (red and green points) have substantially higher bulk velocity indicating that they are not kinematically associated with the disc, but rather outflowing, in particular O2 and L2 (squares). O1 and L1 appear to interact more significantly with the surrounding medium, thus being more prone to turbulent condensation.

In the text

Fig. 9. Comparing Mrk 1044’s integrated ionised gas outflow with AGN multi-phase outflow scaling relations. The coloured data points are taken from Fiore et al. (2017) and show the molecular gas outflows traced by CO, OH (blue), ionised gas outflows (green) and X-ray outflows (red). Panel (a) shows the outflowing ionised gas mass. The dashed coloured lines correspond to the correlations for the respective gas phase reported by Fiore et al. (2017). Panel (b) shows the AGN wind kinetic power where the dashed, dotted and dash-dotted line correspond to the ${\dot{E}}_{\mathrm{kin}}=0.01,0.1,1L_{\mathrm{bol}}$. Panel (c) shows the mass loading factor η = ṀOF/SFR. Here, the brown star represents Mrk 1044’s SF in the immediate vicinity (< 160 pc) from the centre, the orange point to that of the SF in the CNE and the yellow point to that of Mrk 1044’s galaxy-wide SFR. Panel (d) shows the kinetic outflow power as a function of SFR. The solid red line indicates the predicted SFR from the AGN bolometric luminosity using the relation from Netzer (2009) together with the average ${\dot{E}}_{\mathrm{kin}}/L_{\mathrm{bol}}=0.025$ for molecular outflows (Fiore et al. 2017). Compared to the ionised gas outflows in the sample, Mrk 1044’s outflow carries a lot of mass in the ionised gas phase. Its mass injection into the ISM exceeds the current SFR in the immediate vicinity of the nucleus, but also compared to the galaxy-wide values.

In the text

Fig. A.1. Measuring the size of Mrk 1044’s spatially extended Hβ outflow. Panel (a) shows the flux map of the Hβ outflowing component which we extracted as described in Paper I. Panel (b) shows the best-fitting model that consists of the empirical PSF convolved with a Gaussian kernel and a constant background flux. The residuals in panel (c) show that broadened PSF describes the Hβ profile well. Panel (d) shows the χ2 distribution in the parameter plane where a global minimum is present around 0.55 px (4.6 pc).

In the text

	Fig. B.1. Illustration of the expanding shell geometry. The Ly-α occurs in the expanding shell, along the LOS towards the BLR with the velocity v. The [O IIII] emission originates from an intrinsically extended but spatially unresolved shell cap with a half-opening angle θ0, half-light radius α and inclination i. The observer measures the integrated flux and thus measures a smaller velocity.
In the text

	Fig. B.2. Dependence of the emission line flux maximum corresponding velocity on inclination. Independent of how concentrated the light profile is (α/θ0), the geometric model requires high inclinations of the [O III]-emitting shell to significantly reduce the peak velocity by a factor of approximately two.
In the text

Fig. B.3. Best-fitting model for the blue shoulder of Mrk 1044’s [O III]λλ4959,5007 emission line. The left panel shows the AGN spectrum (black line) around the [O III] doublet emission. The best-fit model (orange line) reproduces the line shape well. Both the slow (red line) and the fast shell (blue line) contribute to the asymmetric line shape. The panels on the right show percentile regions of the posterior distribution for the geometric parameters (left) and the flux contribution (right) respectively.

In the text