The Close AGN Reference Survey (CARS) Tracing the circumnuclear star formation in the super-Eddington NLS1 Mrk 1044

Context. The host galaxy conditions for rapid supermassive black hole growth are poorly understood. Narrow-line Seyfert 1 (NLS1) galaxies often exhibit high accretion rates and are hypothesized to be prototypes of active galactic nuclei (AGN) at an early stage of their evolution. Aims. We present adaptive optics (AO) assisted VLT MUSE NFM observations of Mrk 1044, the nearest super-Eddington accreting NLS1. Together with archival MUSE WFM data, we aim to understand the host galaxy processes that drive Mrk 1044’s black hole accretion. Methods. We extracted the faint stellar continuum emission from the AGN-deblended host and performed spatially resolved emission line diagnostics with an unprecedented resolution. Combining both MUSE WFM and NFM-AO observations, we used a kinematic model of a thin rotating disk to trace the stellar and ionized gas motion from 10kpc galaxy scales down to ∼ 30pc around the nucleus. Results. Mrk 1044’s stellar kinematics follow circular rotation, whereas the ionized gas shows tenuous spiral features in the center. We resolve a compact star-forming circumnuclear ellipse (CNE) that has a semi-minor axis of 306pc. Within this CNE, the gas is metal-rich and its line ratios are entirely consistent with excitation by star formation. With an integrated star formation rate of 0 . 19 ± 0 . 05 M ⊙ yr − 1 , the CNE contributes 27% of the galaxy-wide star formation. Conclusions. We conclude that Mrk 1044’s nuclear activity has not yet a ﬀ ected the circumnuclear star formation. Thus, Mrk 1044 is consistent with the idea that NLS1s are young AGN. A simple mass budget consideration suggests that the circumnuclear star formation and AGN phase are connected and the patterns in the ionized gas velocity ﬁeld are a signature of the ongoing AGN feeding.


Introduction
AGN accreting with high Eddington ratios are often classified as Narrow-line Seyfert 1 (NLS1) galaxies.NLS1s offer an unobscured view of the broad-line region (BLR), such that their optical spectra exhibit both emission lines originating from the high-density BLR and from the low-density narrow-line region (NLR).As opposed to Broad-line Seyfert 1s (BLS1s), both permitted and forbidden emission lines are narrow (FWHM(Hβ) < 2000 km s −1 , Goodrich 1989).NLS1s often exhibit weak [O iii] emission ([O iii]/Hβ<3, Osterbrock & Pogge 1985) and strong Fe ii emission from the broad line region (BLR).This locates them at the extreme end of the Eigenvector 1 relation (EV1, Boroson & Green 1992;Sulentic et al. 2000) which is believed to be driven by the Eddington ratio.
The relative narrowness of the broad lines in NLS1 is usually interpreted as result of smaller rotational velocities of ion-ized gas clouds in the BLR.Thus, the inferred central black hole masses of NLS1 are small (10 6 -10 8 M , e.g.Zhou et al. 2006;Rakshit et al. 2017).Although 'narrow' broad lines could also be produced by inclined BLR discs in conventional broad line Seyferts 1s (BLS1s, Baldi et al. 2016), inclination-independent BH mass estimates (e.g.Du et al. 2014Du et al. , 2015;;Wang et al. 2014;Pan et al. 2018;Berton et al. 2021) and the host morphologies of NLS1s (Krongold et al. 2001;Järvelä et al. 2017Järvelä et al. , 2018;;Berton et al. 2019) suggest genuinely low BH masses.
NLS1s systematically lack large-scale diffuse radio emission (Komossa et al. 2006;Gliozzi et al. 2010;Doi et al. 2013;Lister et al. 2016), and AGN-ionized gas on kpc scales in contrast to broad-line Seyfert 1s (BLS1s) which suggests little to no AGN activity in their recent past (e.g.Husemann et al. 2008).It is often speculated that NLS1s represent the early stages of the AGN life cycle where the fuelling of the central super-massive black hole (SMBH) is not yet affected by feedback processes from the host galaxy (Mathur 2000;Collin & Kawaguchi 2004).
Although NLS1s share the same classification based on their optical spectra, they exhibit fairly heterogeneous properties.Only few are detected in the radio among which the dominating emission mechanism ranges from AGN-dominated, composite to host-dominated (Järvelä et al. 2021).Among the AGNdominated ones, some show traces of extended jets detected in the radio, reaching distances of tens of kpc from the nucleus.In individual NLS1 galaxies relativistic jets powered by the central engine have been found (Yuan et al. 2008;Foschini 2011;Foschini et al. 2015) as well as blazar-like phenomena such as high brightness temperature, a double-humped spectral energy distribution or gamma-ray emission (Romano et al. 2018;Paliya et al. 2019;Komossa 2018).X-ray observations have shown that NLS1s exhibit fast variability on typical time scales of less than one day (Boller 2000) which is expected for the small black hole masses and high accretion rates (Ponti et al. 2012).
One of the open questions regarding the place of NLS1s among the AGN population is the origin of the exceptionally bright Fe ii emission from their BLR.The high metallicity is difficult to explain from a galaxy evolution perspective as it requires substantial star formation (SF) to enrich the material in the accretion disk.At redshifts around z∼1 SF galaxies are more likely to host rapidly-growing, radiatively-efficient AGN despite having similar gas fractions and star formation efficiencies as their non-SF counterparts (Nandra et al. 2007;Goulding et al. 2014).This suggests that both galaxy-scale SF and AGN phase are triggered by the same secular processes.At intermediate and low redshifts however, the accretion becomes radiatively inefficient and mostly mechanical (Hickox et al. 2009).Furthermore, it is matter of discussion how gas is channeled from the galactic scales down to the sub-pc scales of the circumnuclear accretion disc, how the gas transport affects the SF and its spatial distribution in the host, and what is the main driver of such transport (e.g., turbulent mixing, inelastic collisions, gravitational torques, see Gaspari et al. 2020 for a review).
The host galaxies of NLS1s are predominantly classified as disc galaxies and often possess pseudo-bulges (Crenshaw et al. 2003).They often exhibit enhanced SF in their nuclear regions (Deo et al. 2006) and only a few of them show signs of major mergers or interaction (Ohta et al. 2007).Therefore, the gas transport towards the galaxy center requires internal process to prolong star formation and the AGN accretion.This idea is further supported by the frequent presence of pseudo-bulges in NLS1 host galaxies (Orban de Xivry et al. 2011) which are formed from secular processes (Kormendy & Kennicutt 2004).These asymmetries can develop circumnuclear spiral structures (Maciejewski 2004) and enable radial gas migration (Sakamoto et al. 1999;Sheth et al. 2005).
Local AGN in which circumnuclear star formation has been resolved have low Eddington ratios (e.g.Esquej et al. 2014;Ramos Almeida et al. 2014;Ruschel-Dutra et al. 2017;Esparza-Arredondo et al. 2018;Knapen et al. 2019) and may therefore have been caused by different fuelling mechanisms than accretion-mode NLS1s.The situation is further complicated by the difficulty to observe close to their unobscured nucleus.It requires an accurate deblending of the AGN emission that is broadened by the instrumental point spread function (PSF) and the faint host emission.Moreover, the NLS1 population is heterogeneous (Komossa et al. 2006;Järvelä et al. 2017) and low-redshift NLS1s are rare which makes a statistical evaluation of their host galaxy properties challenging.Given the diversity among NLS1s and their complex behaviour, their classification by their opti-cal classification alone does probably not reflect the underlying physical mechanisms.We therefore need in-detail studies of individual objects to explain their characteristics on the black hole -host galaxy interaction.
Mrk 1044 is a textbook example of a NLS1 and one of the nearest super-Eddington accreting AGN.It is a radio-quiet, luminous (L bol = 3.4 × 10 44 erg s −1 ) NLS1 at z = 0.0162 (Husemann et al. 2022).Its central engine is powered by a black hole with a reverberation mapped mass of M • = 2.8 × 10 6 M (Du et al. 2015).Numerous studies have confirmed Mrk 1044's high accretion rate with Eddington fractions λ Edd = L bol /L Edd that range from 1.2 (Husemann et al. 2022) up to 16 (Du et al. 2015).Although Mrk 1044 is a bright X-ray source (Dewangan et al. 2007), it shows no signs of AGN-driven hot outflows in the form of extended X-ray emission (Powell et al. 2018).The luminous central region, however, shows frequency-dependent variability and an X-ray excess that is likely caused by relativistic Compton scattering from a high-density accretion disc (Mallick et al. 2018).In a recent study, Krongold et al. (2021) used XMM-Newton observations to identify an ultra-fast outflow with a multi-velocity, multi-phase structure originating at the scales of the accretion disk.
On galaxy scales, Mrk 1044 has a barred spiral morphology (Deo et al. 2006) with an estimated H i mass of 2.6 × 10 9 M (König et al. 2009) and H 2 mass of 4 × 10 8 M (Bertram et al. 2007).Powell et al. (2018) have identified three concentric star forming rings.The outer ring at ∼ 7 kpc is likely to be associated with spiral arms in a rotating disk, whereas the inner SF region has not been resolved in earlier studies.Although Mrk 1044 is an extensively studied nearby NLS1, a complete picture of the host galaxy processes and how they affect its AGN activity is still missing.In this work we present adaptive-optics assisted optical integral field spectroscopic (IFU) observations of Mrk 1044 to study the ionized properties and explore the kinematic signatures of the black hole -host galaxy interaction at unprecedented spatial resolution.Throughout this paper we assume a flat ΛCDM cosmology with H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3, and Ω Λ = 0.7.In this, 1 corresponds to 0.333 kpc at the redshift of Mrk 1044 and the associated luminosity distance is 70.0 Mpc.

Observations and data reduction
For this study we use the integral field spectroscopic data acquired with the Multi Unit Spectroscopic Explorer (MUSE, Bacon et al. 2010Bacon et al. , 2014) ) at the Very Large Telescope (VLT).The unprecedented spatial resolution and spectral coverage of the adaptive-optics assisted narrow field mode (NFM-AO) allow us to map the host galaxy emission line properties and kinematics of both stellar and gas components down to a resolution of 89 mas (30 pc).In order to gain a complete picture of the physical processes on different spatial scales, we combine them with the IFU data from the Close AGN Reference Survey (CARS, Husemann et al. 2022) and a high resolution broad band image acquired with the UVIS instrument on the Hubble Space Telescope (HST).

Optical Imaging
Mrk 1044 has been imaged with HST under the program ID 12212 using WFC3/UVIS in the F547M filter.We collect the archival image from the Hubble Legacy Archive1 .The field of view (FOV) spans 3.66 ×3.56 with a spatial resolution of 0 .067 at 6000Å as reported in the WFC3 Instrument Handbook.

VLT MUSE WFM
As part of the Close AGN Reference Survey (CARS, Husemann et al. 2022)2 Mrk 1044 has been observed with the Wide Field Mode (WFM) of the MUSE instrument.The data were reduced using the standard ESO pipeline v.2.8.1 (Weilbacher et al. 2012(Weilbacher et al. , 2020) ) as described in Husemann et al. (2022).The resulting data cube has a FOV of 62 .7 × 63 .5 with a sampling of 0 . 2 per pixel and a seeing-limited resolution of 1 .03.Along the wavelengthaxis it extents from 4750 Å to 9300 Å with a spectral resolution of ∼ 2.5 Å that slightly varies over the spectral range (Bacon et al. 2017;Guérou et al. 2017).

VLT MUSE NFM-AO
We observed Mrk 1044 with the AO-assisted NFM of the MUSE Integral Field Unit under the program 0103.B-0349(A).The data were acquired on 23 Aug 2019, 24 Aug 2019 and 28 Sep 2019.The 12 exposures have an integration time of 550 s each and were dithered by a small angle to minimize the imprint of cosmic rays and flat-fielding artifacts.We use the MUSE pipeline v2.8.3-1 together with the graphical user interface ESO Reflex v2.11.0 to execute the EsoRex Common Pipeline Library reduction recipes.
For WFM data, the reduction pipeline performs a differential atmospheric refraction (DAR) correction that is based on a theoretical prediction (Weilbacher et al. 2020).However, the correction is poor for observations in the NFM since the performance of the hardware build-in atmospheric dispersion correction (ADC) is insufficient.In order to correct for the remaining wavelength-dependent atmospheric refraction, we use a 2D Moffat model3 to measure the AGN position of the individual exposures after the sky correction (muse_scipost).To remove the wavelength-dependence, we subtract the wavelength-dependent coordinate offset in the PIXTABLE files.Fig. 2 demonstrates the procedure for one example exposure.After the correction, the residual variation of the QSO position across the full wavelength range is of the same order of magnitude as its scatter (∼ 0 .04).
The spatial resolution of individual exposures depends on the wavelength-dependent performance of the AO system.The science observations taken on 24 Aug 2019 have inferior quality compared to the remaining ones due to the atmospheric conditions during the exposures.Especially at the blue end of the spectrum where we perform most of our diagnostics the difference in resolution is substantial.Here, the FWHM of the 2D Moffat to model the AGN PSF is larger by a factor of 3 − 4. In order to not degrade the resolution of the data cube, we therefore only combine the 8 exposures taken in the night of 23 Aug 2019 and 28 Sep 2019 with muse_exp_combine.
The final data cube consists of 369×378 spaxels, corresponding to a FOV of 9 .23 × 9 .45 with 127,741 spectra.From the width of the telluric emission lines we find a constant resolution of FWHM = 2.54 ± 0.10 Å across the full wavelength range 4750 Å -9350 Å, corresponding to 160.4 km s −1 and 81.5 km s −1 at the blue and the red end of the spectrum respectively.

Deblending the AGN and host emission
In 3D spectroscopic observations of type 1 AGN the light of the unresolved central QSO distributes over the host galaxy as dictated by the PSF.Due to the bright AGN and the low surfacebrightness host emission, the AGN contamination is particularly strong within the small FOV of MUSE NFM data cube (see right panel of Fig. 1).Disentangling the host emission from the QSO contamination requires a deblending in both spectral and spatial dimensions.
The AGN-host deblending process requires five steps: 1) estimation of the PSF at the broad lines by AGN-host decomposition, 2) reconstructing the AO-shaped PSF at the broad lines with a model, 3) generating a hybrid PSF that combines the empirical PSF with the model PSF, 4) interpolating the hybrid PSF along the wavelength axis, 5) applying an iterative AGN-host deblending procedure combining the wavelength-dependent hybrid PSF with the host galaxy surface brightness profile.In the following subsections, we describe each of the steps in detail.

PSF extraction at the broad emission lines
In order to separate the host emission from the QSO emission, we employ the dedicated software 4 QDeblend 3D (Husemann et al. 2013).It is based on the concept that the spectrum of each spaxel is a superposition of the host spectrum and a spectrum from the central QSO that is scaled in flux according to the PSF.Since broad lines from the BLR are spatially unresolved, i.e. point like, type 1 AGN allow to extract the empirical PSF at each broad line available.In the following, we give a brief outline of the concept.For a detailed description of the algorithm we refer the reader to Husemann et al. (2013).
As a first step, QDeblend 3D scales the broad-line wings to the flux of the pseudo-continuum to obtain a PSF.The high S/N QSO spectrum is then scaled to match the flux of the broad line wings in each spaxel (step 1), which corresponds to a tentative QSO data cube.As a next step, the emission from the host galaxy is obtained by subtracting the QSO data cube from the original data cube.Since the initially extracted QSO spectrum is inevitably contaminated by a small fraction of host galaxy light this leads to an over-subtraction of the host cube near the center.To estimate the host contribution, QDeblend 3D interpolates the host cube towards the central QSO position assuming a surface brightness (SB) profile.In our case we use a constant SB profile since the over-subtracted region is small (∼ 0 .5).The host spectrum in the central pixel is then subtracted from the tentative 4 https://git.io/qdeblend3dQSO spectrum in order to obtain an updated QSO spectrum that is used for the next iteration starting at step 1.
With this procedure we extracted the empirical PSF at Hβ, Hα and O iλ8446+Ca iiλ8498 (see Matsuoka et al. 2007) which have a FWHM of 89, 51 and 41 mas respectively.

Modelling the MUSE NFM-AO PSF
Our observations of Mrk 1044 with MUSE NFM were assisted by the Ground Atmospheric Layer Adaptive Optics for Spectroscopic Imaging (GALACSI) adaptive optics system (Ströbele et al. 2012) to improve spatial resolution.The AO-shaped PSF does not follow the Moffat profile suited for seeing-limited observations.Instead, it has a peculiar shape that depends on the performance of the AO system.Here, we describe the PSF modelling for bright point-like AGN emission required for the subsequent deblending process.
The NFM-AO PSF only slowly changes with wavelength.However, since we subtract the bright AGN emission to get the relatively faint host, already small inaccuracies in the PSF model could severely affect our analysis of the narrow host emission lines.Another problem is that the subtraction of an empirical PSF adds noise to the host data cube which dominates over the faint host emission already at small distances from the center.We are therefore interested in using a model that describes both the behaviour of the NFM-AO PSF towards the faint outskirts and its wavelength dependence.An analytic model for the long-exposure AO-corrected PSF has been presented by Fétick et al. (2019, hereafter PSFAO19).They describe the phase power spectral density (PSD) with a narrow Moffat core and a wide Kolmogorov halo.The PSFAO19 model consists of 5 parameters (+2 in the asymmetric case) which include the Fried parameter r 0 , the AO-corrected phase PSD background C, the Moffat amplitude A and the Moffat scale factors α and β.This parametrization is physically motivated by the design of the AO system.The actuators that deform the secondary mirror are separated by a pitch which sets the maximal spatial frequency of the phase f AO that can be corrected by the AO system.Beyond this frequency the residual phase PSD is not affected by the AO system which leaves a residual turbulence in the PSF outskirts.
To retrieve the best-fit parameters for the empirical PSFs available, we use the PSFAO19 model that uses a Levenberg-Marquardt algorithm to minimize the χ 2 -sum of the flux in the residual images.We fit the PSF over the full FOV, weighted by the inverse of the noise variance of the individual pixels.In Table 1 the best-fit parameters are listed.The Fried parameter r 0 increases with wavelength as λ 1.22±0.05and is consistent with Table 1.Best-fit PSFAO19 model parameters at three different wavelengths.We extracted the empirical PSFs from the data cube at each of the available broad lines.broad line  the theoretical prediction λ 6/5 .Since the atmospheric refraction is smaller in the red (i.e. the AO system performs better), the FWHM of the PSF is decreases with wavelength.

Generating the hybrid PSF
In contrast to the empirical PSF, the analytic model has the advantage that it does not contain noise.Thus, it does not add noise to the host cube when subtracting the PSF cube from the original data cube.This is particularly important at large distance from the center where the host emission is faint and the S/N of the empirical PSF is low.We find that the outskirts of the PSF are well reproduced by the PSFAO19 model (see Fig. 3).Near the center, however, the amplitude of the residuals is ∼ 1% of flux from the AGN and therefore of the same order of magnitude as the host signal.
In order to combine the advantages of both regimes, we create a hybrid PSF.For the high S/N central regions we use the empirical PSF, whereas the model PSF is superior at large distance from the center.We select the transition at the radius beyond which the azimuthally averaged S/N falls below a threshold value of 3.For Hβ, Hα and O i+Ca ii the transition radii are located at 39, 73 and 28 px respectively.

Interpolation of the PSF with wavelength
The PSF is wavelength dependent as can be seen from Fig. 4. Therefore, we need to interpolate the 2D PSFs that are available for the three broad lines listed in Table 1.
As a first approach, we use a spline interpolation of the parameters in Table 1 to generate a PSF cube across the full wavelength range.However, we find that this method leaves strong residual artifacts between the broad lines at which the PSF was extracted.This indicates that not all of the PSFAO19 parameters change smoothly with wavelength and only three broad lines are insufficient to constrain the exact wavelength-dependence.
An alternative approach involves a simple but effective pixelby-pixel interpolation of the hybrid PSF across the full wavelength range.Since Hβ and O iλ8446+Ca iiλ8498 are close to the edges of the MUSE wavelength range, no significant extrapolation is required.As a first step we normalize the hybrid PSFs extracted at Hβ, Hα and O iλ8446+Ca iiλ8498 to their peak flux.For each of the spaxels, we describe the wavelength dependence of the flux with a third order polynomial.Consequently, the PSF cube equals the available empirical PSF at the positions of the broad lines.The interpolated PSF cube is then scaled in the central spaxel to match the flux of the QSO spectrum.

Iterative AGN-host deblending from the IFU data
We now want to deblend the AGN from the host emission.By construction, the PSF cube contains the faint but non-negligible emission from the host in the central pixel.Therefore, the subtraction from the original datacube causes an over-subtraction in the center.We follow the AGN-host deblending approach described in Husemann et al. (2022) where the AGN spectrum is iteratively corrected for the host galaxy contribution.Thereby we assume a constant host galaxy surface brightness within the innermost 0 .5 (corresponding to 167 pc or 0.08 r e ).
The result of the iterative deblending process is displayed in Fig. 5.The 3 aperture spectra show that the point-like AGN emission contains the entire broad line emission and little to no contribution from narrow lines.This is because Mrk 1044 AGN is so luminous that it dominates over the narrow line emission on pc-scales.The host spectrum contains the stellar absorption and the spatially extended ionized gas emission in the form of narrow absorption and emission lines respectively.
Most of our diagnostics in the following are based on emission lines in the Hβ-[O iii] window.We therefore estimate the spatial resolution for the analysis from PSF cube at Hβ where the FWHM is 89 mas which corresponds to 30 pc in the galaxy system.

Analysis & Results
We aim to quantify the stellar and ionized gas properties from the IFU data and map them across the host galaxy.After the AGNhost deblending described in the previous section, the following analysis only employs the AGN-subtracted host emission.

Fitting the stellar continuum and the ionized gas emission
For the extraction of the galaxy properties we employ PyParadise (Husemann et al. 2016(Husemann et al. , 2022) ) which is a publicly available 5 updated version of the stellar population synthesis fitting code Paradise (Walcher et al. 2015).It is extended by an Markov chain Monte Carlo (MCMC) algorithm to find the parameters of the line-of-sight velocity distribution (LOSVD) and a routine to fit the emission lines.The main steps involved are the following.
(1) As a first step, both the input spectra and the template library spectra are normalized by a running mean where emission 5 https://git.io/pyparadiselines or strong sky line residuals are masked and linearly interpolated.This normalization has the advantage that the fitting result is less sensitive to non-physical continuum variations caused by the wavelength interpolation of the PSF, especially near the nucleus (see Fig. 5).
(2) In the next step, PyParadise uses an iterative scheme to independently determine the best-fit of the LOSVD and linear combination of simple stellar population (SSP) template spectra.Starting from an initial LOSVD guess obtained by fitting a single spectrum drawn from the template library, all spectra in the template library are convolved with the initial LOSVD guess to obtain the best-fit non-negative linear combination of template spectra.The resulting best-fit spectrum is used as the input for further iterations and the two step process of estimating the LOSVD and linear combination is repeated.For our analysis we select the CB09 SSP library, which is an updated version of the library presented in Bruzual & Charlot (2003).The template metallicities range from [Fe/H] = -1.44 to 1.44 and stellar ages from 1.7 Myr to 13 Gyr.The SSP spectra cover 3500 Å to 9500 Å at a spectral resolution of 2.5 Å.Our results do not change within the uncertainties if we employ the higher-resolution SSP library from González Delgado et al. (2005).
(3) The final best-fit spectrum is denormalized and subtracted from the original spectrum.
(4) To model the emission lines from the residual spectrum, PyParadise uses a set of Gaussian models with a common line-of-sight velocity, taking into account the spectral resolution that is approximately constant with wavelength (see Sect. 2.3).For the doublet emission lines [N ii]λλ6548,83 and [O iii]λλ4959,5007 we fix the flux ratio to the theoretical prediction of 2.96 (Storey & Zeippen 2000;Dimitrijević et al. 2007).Furthermore, we couple the emission lines in radial velocity and velocity dispersion in order to increase the robustness of the flux measurement.Our results do not change within the uncertainties if we do not kinematically tie the model parameters for emission lines with different ionization potentials.
Compared to the established pPXF software (Cappellari & Emsellem 2004;Cappellari 2017) which uses low-order polynomials to fit the stellar continuum, PyParadise allows us to model the stellar absorption lines irrespective of the details of the PSF subtraction.In this way we can robustly constrain the stellar kinematics close to the nucleus.In Fig. 6 we visualize an example spectrum from the MUSE data cube within a randomly selected aperture, together with the best-fit spectrum of the stellar continuum and the ionized gas emission.For a large fraction of the spaxels the S/N is too low to robustly model the low surface brightness emission lines and the faint stellar continuum emission.We therefore employ two spatial binning techniques that spatially co-add the spectra of the host data cube.For the stellar continuum we use the adaptive Voronoi tessellation routine of Cappellari & Copin (2003) to achieve a minimum signal-to-noise ratio of S/N=20 in the wavelength range 5080 Å < λ < 5260 Å.We then model the binned stellar continuum spectra with PyParadise using the updated CB09 version of the evolutionary synthesis model spectra from Bruzual & Charlot (2003) before projecting the stellar kinematics from the Voronoi-grid onto the initial MUSE sampling grid.This information is then used to model the stellar emission again for each spaxel of the host data cube, but this time keeping the stellar LOSVD parameters fixed.Finally, we spatially co-add 2×2 and 8×8 spaxels of the host data cube and fit the emission lines to the residual continuum spectra.The spatial binning allows us to extract the emission lines properties in regions that are more than one order of magnitude fainter than in the original data cube.We estimate the uncertainties for all emission line parameters with a Monte-Carlo approach.Both stellar continuum and emission line fitting are repeated 40 times after each spectrum is fluctuated randomly within the error of each spectral pixel.

Mapping the emission line properties
The surface brightness maps together with the velocity field and the dispersion are shown in Fig. 7.For the surface brightness maps we only select spaxels with S/N>3, and uncertainty of < 30 km s −1 for both velocity and dispersion.For the high S/N regions we overplot the 8×8 binned map with the 2×2 binned spaxels.The emission from the narrow lines is mostly concentrated in a nearly circular patch around the nucleus.Beyond ∼800 pc from the center, the signal from the emission lines is too weak to be robustly modelled.Within the patch, we detect a pronounced ring-like structure that is present in each of the Hβ, Hα and [N ii]λ6583 emission.The ring structure is most prominent in the Hα emission for which we extract the contour that is shown in the bottom left panel of Fig. 7. Since the inclination of the galaxy disk is 67 • (Powell et al. 2018), the deprojected ellipse has an axis ratio of 2.45.The size of the ring compares to the extent of the luminous inner region in HST WFC3/UVIS image, for which we show the contours in the top left panel of Fig. 1.
For each of the surface brightness maps the peak emission stems from near the nucleus.The associated high velocity and the high dispersion in that region indicate an outflowing component.We note that although the AGN-host deblending algorithm is capable of subtracting any compact emission, it is only a first order subtraction of the point-like emission.As demonstrated by Singha et al. (2022) the region near the nucleus may require a more complicated modelling with multiple narrow line emitting structures.Here and in the following analysis, we exclude the emission from the innermost 0 .5 (160 pc).A detailed analysis and discussion of Mrk 1044's outflow will instead be presented in Winkel et al. (in prep.).Ignoring the high velocity/dispersion feature in the center, the velocity field exhibits a smooth rotational field across the MUSE NFM FOV.The median dispersion of σ = 34.9km s −1 is fairly small, indicating dynamically cold ionized gas.

Excitation mechanism
From the emission line ratios the underlying excitation mechanism can be identified, i.e. star formation ionization in H ii regions, photo-ionization by the hard radiation field of an AGN and low-ionization nuclear emission-line regions (LINERs).A commonly used demarcation is made in the Baldwin-Phillips-Terlevic diagram (BPT, Baldwin et al. 1981) that uses the [O iii]/Hβ vs. [N ii]/Hα line ratios.The upper right panel of Fig. 7 shows the BPT diagram for the 8×8 binned spaxels.All spaxels are associated with a cloud that spreads over the SF domain but is elongated towards the AGN ionized region.The scatter originates from the finite spatial resolution of 30 pc which does not allow to resolve individual SF clouds.Together with the lineof-sight projection, different ionization conditions are inevitably captured within one spaxel.Furthermore, different ISM metallicities and densities lead to the star forming cloud occupying a large area in the BPT diagram (e.g.Smirnova-Pinchukova et al. 2022).Nonetheless, the vast majority of spaxels in the MUSE NFM FOV is consistent with star formation excitation.We therefore assume the AGN contribution in terms of ionizing radiation to be negligible.
The spatial distribution of the excitation mechanism is shown in the bottom right panel of Fig. 7. Interestingly, the spaxels classified as composite are located towards the outskirts of the bright Hα ring.This could either be caused by evolved stars during the short but energetic post-AGB phase which mimic the LINER emission (Binette et al. 1994;Stasińska et al. 2008), leaking photons from HII regions plus starburst driven shocks (Heckman 1980;Dopita & Sutherland 1996).More strikingly, it becomes evident that even in the very center of Mrk 1044, SF appears to be the dominant excitation mechanism in the host.Since the inner region is dominated by an outflow, we follow the procedure described in Singha et al. (2022) and use a two-component Gaussian model to fit the central 0 .5-integrated host spectrum.The resulting BPT classification for the host component is consistent with the single-component analysis of the individual spaxels (see Fig. 7).This confirms the success of our AGN-host deblending process, since an incorrect treatment of the PSF would generate residual spatially extended AGN ionization structures.On the other hand, it contradicts the findings in nearby AGN where the high excitation is typically found close to the nucleus (Davies et al. 2014;Richardson et al. 2014).This is particularly interesting given Mrk 1044's super-Eddington accretion rate.At such high luminosities quasar mode feedback is expected to drive galactic scale outflows through radiation pressure from the active nucleus (Nesvadba et al. 2007;Liu et al. 2013;Cicone et al. 2018).However, Mrk 1044's host does not exhibit any sign of interaction with the AGN photo-ionization field < 1 kpc despite the current AGN phase.If its BLR was obscured, Mrk 1044 would be identified as star forming galaxy.Moreover, the line ratios from the AGN-contaminated 3 -aperture spectrum clearly place Mrk 1044 in the star forming regime (top right panel of Fig. 7).There exist several examples of X-ray selected obscured AGN that appear optically inactive (XBONG galaxies, e.g.Severgnini et al. 2003;Caccianiga et al. 2007;Trump et al. 2009) and star forming in their diagnostic line ratios (Hornschemeier et al. 2005).Castelló-Mor et al. ( 2012) have shown X-ray bright galaxies optically classified as star forming are mostly NLS1s.This suggests that Mrk 1044 could be a prototype of many accretion-mode AGN at high redshift that are difficult to identify since their spatially unresolved narrow line ratios classify them as star forming galaxies.

Circumnuclear star formation
The star formation rate (SFR) is an important parameter for understanding the interaction between the central AGN and its host galaxy.The surface brightness maps derived in the previous section allow us to map the current star formation rate surface density Σ SFR with high spatial resolution.
As a first step we have to correct the Hα flux for foreground dust attenuation along the line-of-sight.Assuming case B recombination at an electron temperature T e = 10 4 K and electron density n e = 100 cm −3 the intrinsic Balmer ratio is 2.86 (Osterbrock 1989).The nebular colour excess can be estimated from the Balmer decrement as Here, κ(λ Hβ ), κ(λ Hα ) is the reddening curve for star forming galaxies from (Calzetti et al. 2000, hereafter C00) evaluated at Hβ and Hα, respectively.Using the C00 reddening curve, the dust extinction at wavelength λ is connected to the colour excess by With the above relation we correct the Hα luminosity as L intr (Hα) = L obs (Hα) × 10 A Hα 2.5 (5) Following the argument from the previous section, we assume that all of the extinction-corrected Hα luminosity is generated by SF ionization.We then compute the SFR using the relation from Calzetti et al. (2007, C07) which was calibrated for nearby star-forming galaxies assuming a constant SFR over 100 Myr

SFR [M yr
This calibration applies for the default two power law stellar IMF from Starburst99 (Leitherer et al. 1999).Moreoever, it contains the underlying assumption that the ISM has solar metallicity.For a robust measure of the Σ SFR we impose a minimum S/N for both Hβ and Hα line fluxes.Further, we select E(B − V) = 0 for spaxels with F Hα /F Hβ < 2.86.After applying these criteria the selected 8 × 8 binned spaxels sample the entire star forming ring, i.e. we do not lose flux when computing integrated properties of the SF ring.
The left panel of Fig. 8 shows that the extinction is particularly high in the ring of enhanced Hα emission.Within the Hα-contour (and excluding the innermost 0 .5, see Sect.3.2) the median extinction is 0.65 mag, indicating significant abundance of dust along the line-of-sight.This coincides with the presence of dusty spiral absorption features in the HST WFC3 image (Fig. 1).As shown in the right panel of Fig. 8, the Hα-ring shows particularly high SFR surface density.The azimuthally averaged peak SFR located at a projected radius of r = (306 ± 43) pc.Taking into account the galaxy inclination of 67 • , the deprojected structure is an ellipse with an axis ratio of 2.45.We therefore refer to the circumnuclear region of enhanced SFR as the circumnuclear ellipse (CNE).As before, we exclude the innermost 0 .5 when computing integrated properties of the CNE.
The average Σ SFR within the CNE is 0.16 ± 0.04 M yr −1 kpc −1 , as computed from the integrated extinction-corrected Hα flux.This value lies significantly above the minimum of 10 −3 M yr −1 kpc −2 for driving galactic winds (Veilleux et al. 2005) and is typical of circumnuclear star formation bursts (Kennicutt 1998;Kennicutt & De Los Reyes 2021).It is comparable to the lower end of the Σ SFR in ultra-luminous infrared galaxies (ULIRGS) that host the most powerful starbursts known (Dopita et al. 2002;Genzel et al. 2010).
For typical star forming spiral galaxies with Σ gas = 10 M pc −2 the Kennicutt-Schmidt law predicts Σ SFR = 6 × 10 −3 M yr −1 kpc −1 .Following the same procedure as for the CNE, we use the MUSE WFM data cube to compute the galaxy-wide SFR surface density of Σ SFR = (5.9± 1.0) × 10 −3 M yr −1 kpc −1 .We estimate an integrated value of SFR tot = 0.70 ± 0.17 M yr −1 .This value is consistent with SFR = 0.6 ± 0.2 M yr −1 estimated by Smirnova-Pinchukova et al. ( 2022) who used the integrated IR luminosity.While the galaxy-wide value SFR places Mrk 1044 among the star forming population, it is evident that the inner star forming ellipse with SFR CNE = 0.19 ± 0.05 M yr −1 accounts for 27% of the absolute galaxy star formation rate.

Metallicity of the star forming CNE
Another important characteristic of the ionized gas is the metallicity.In contrast to the star formation rate, which traces shortlived (∼5 Myr) emission from H ii regions, the gas-phase metallicity provides a cumulative measure for the enrichment history of the ISM.Furthermore, in-and outflows can affect the chemical composition of the gas on short timescales.The gas-phase metallicity therefore provides an important diagnostic to constrain its origin and motion.
There exist several metallicity indicators that are based on the line ratios emitted by H ii regions.Since Mrk 1044's star forming CNE is completely consistent with excitation by star formation, we can use strong line ratios covered by the MUSE wavelength range to estimate the gas-phase metallicity.The N2S2 index defined as N2S2 = log ([N ii]λ6583/[S ii]λλ6717,31) has proven to be a reliable indicator for the O/H abundance ratio (Dopita et al. 2016).
We use N2S2 index calibration established in Husemann et al. (2019).Since the O/H abundance is significantly dependent on the empirical calibration adopted (Kewley & Ellison 2008), we use the emission line ratios of star forming galaxies in the SDSS DR7 value-added catalog Brinchmann et al. (2004) to determine a N2S2-calibration We thus have a self-consistent measure based on the metallicities from Tremonti et al. (2004).Our results do The dynamical time estimate compares to the timescale of ISM enrichment by massive stars, indicating that the ISM is well mixed on < 1 kpc scales.However, the CNE may be affected from gas in-and outflows occurring on longer timescales.We therefore interpret the constant gas-phase metallicity to be a signature of a homogeneous enrichment process through the ongoing starburst.

Kinematic modelling
To investigate whether the gas and stellar motion are consistent with dynamically cold rotation on a circular orbit, we use a simple model to fit the 2D velocity maps extracted in Section 3.1.Let's consider the first moment of the LOSVD, the intensityweighted mean velocity, which we assume to trace the material in a thin disk.If the inclination of the system i is unknown, the projected velocity of a circular orbit in the plane of the galaxy can be expressed through the azimuthal angle φ measured from the projected major axis in the galaxy plane: Here, v sys is the systemic velocity, PA is the position angle of the rotation axis and (x, y) are the Cartesian coordinates in the observed plane with the galaxy center located at (x c , y c ).We assume the AGN position to be the same as the kinematic center of the circular motion.Furthermore, we adopt the inclination of the disk i = 67.4± 4.5 • computed by Powell et al. (2018).Thus, our analytic model consists of three free parameters (v sys , PA, v) for which we determine the best-fit values using the Bayesian MCMC software GAStimator6 .This algorithm maximises the log-likelihood of the residuals between the observed and the model LOSVD.In order to trace the rotation curve v(r), we sample the FOV with concentric rings within which the model is evaluated.We chose their radial width to equal the FWHM of the PSF, i.e. 1 .03 for MUSE WFM and 89 mas for NFM-AO respectively.Finally, we assume a smooth rotation curve by radially interpolating the parameters (PA, v) to compute the model velocity in each pixel.

Kinematics of the stellar component
Although the stars are likely to have significant pressure support, we conduct a simple experiment here by assuming purely circular motions.The stellar emission is spatially binned such that the stellar continuum reaches S/N≥20 as described in Sect.3.1.Since the bin size varies across the FOV, we discretize the model by projecting it onto the Voronoi grid and taking the average velocity within each bin.We then use this information to maximize the log-likelihood of the LOSVD.The top panels of Fig. 10 shows the stellar velocity field together with the best-fit model and the residuals.Both observations from MUSE WFM and NFM-AO show a smooth rotational pattern that ranges from the smallest resolved scales of ∼100 pc to galaxy scales ∼20 kpc.The thin rotating disk model reproduces this pattern well as we do not see systematic extended structures in the residual map.This indicates that the stellar velocity field is largely unperturbed, justifying our initial (simplistic) assumption of circular orbits for the stellar rotation.Approaching the center, the stars have decreasing rotational velocity v, while the PA stays constant.This is compatible with the quiescent rotating velocity field of a disk galaxy that is dominated by the gravitational potential of the host galaxy alone.

Kinematics of the ionized gas phase
Compared to the stellar rotation, the ionized gas rotates significantly faster than the stellar component as shown in Fig. 9. Since the stars experience a pressure support from their velocity dispersion, their motion in a disk system is not described by circular orbits.The higher rotational velocity of the ionized gas component is a reflection of this effect, the so-called asymmetric drift.
We model the ionized gas velocity on the native MUSE pixel sampling since the velocities measured from the high S/N emission lines are well constrained.Since the high velocity gas in Fig. 10.Modelling the stellar (top) and ionized gas velocity field (bottom) with a thin rotating disk.For each of the two the panels are organized as following: From left to right the panels show (1) the contours of the HST WFC3/UVIS image (2) Mrk 1044's stellar/ionized gas velocity field (3) the kinematic model of a rotating disk and (4) the difference between the observed and model velocity field.For each stars and ionized gas respectively, the top panels show the MUSE WFM data and the bottom panels the zoom-in onto the NFM-AO data.For the ionized gas velocity field we exclude the innermost 0 .5 where the gas does not co-rotate with the disk.Both the stellar and ionized gas are well reproduced by the rotating disk model.Only the ionized gas shows tenuous clumpy structures in the velocity residuals which are highlighted in Fig. 11. the center is clearly not associated with the rotating disk of the host, we use a circular aperture to mask the pixels in the innermost 0 .5. Furthermore, the ionized gas velocity field captured by MUSE WFM exhibits an arc-shaped stream perpendicular to the galaxy disk that stretches from 1.3-4 kpc.Powell et al. (2018) suggested that this 'bloom' was first ejected by the high star formation surface density in Mrk 1044's CNE.Since it does not co-rotate with the galaxy disk, we also exclude this feature from the model.
Apart from the two features that we manually excluded, the ionized gas velocity field appears quiescent, indicating an undisturbed rotational motion.The bottom right panels of Fig. 10 show that on galaxy scales the ionized gas is surprisingly well described by the disk rotation.Only within the innermost 1 kpc we detect spatially extended structures in the velocity residuals that coincide with the location of the star forming CNE.We highlight these extended structures in Fig. 11, where the ionized gas velocity residuals show both receding and approaching components with an spiral-like shape.We select two prominent patches  which were obtained from smoothing the map with a Gaussian kernel of 3 px and imposing a minimum velocity of ±3 km s −1 .The contours of the smoothed map together with the original map are shown in Fig. 11.The size and orientation of the clumps coincide with the direction and extent of the spiral pattern of the dust absorption seen in the HST WFC3 image (see left panel of Fig. 1).The arcs reach from 400 pc to the 150 pc towards the central region as highlighted in Fig. 11.Both regions are ionized by star formation and have symmetric line shapes which indicate that the velocity residuals are genuine rather than artefacts from fitting asymmetric emission lines.The approaching region has a median velocity offset of ∆v r = 5.9 km s −1 and the receding region ∆v r = −7.0km s −1 .Clump-like structures are typical sig-natures of radial streaming motions within the disc, likely due to slow inflow in the disc plane.
Fig. 12 shows the radial dependence of the kinematic parameters measured with the rotating disk model.As expected for an inclined rotating disk, the rotational velocity increases with increasing distance from the center up to a maximum value of v gas ≈ 100 km s −1 and v ≈ 80 km s −1 respectively before the rotation curve flattens.Also the model confirms that the ionized gas rotates faster at all radii which is expected due to the asymmetric drift (see Fig. 9).Gadotti et al. (2020) report that the nuclear discs have larger rotational support as compared to the underlying main disc.For Mrk 1044 we recognize a velocity jump around ∼ 200 pc where both the ionized gas and the stellar rotation appear to rotate faster in the CNE.However, the seeing-limited resolution of the MUSE WFM data of FWHM=1 .03, the rotation velocity in the innermost ∼ 330 pc is underestimated due to beam smearing.Only MUSE NFM-AO allows to resolve the stellar motion on scales of the CNE, but doesn't have the coverage to compare with larger scales.Thus, the velocity jump could also be caused by the different resolution of MUSE WFM and NFM-AO.Within the uncertainties, the PA of the rotational axis is constant across the whole galaxy.The median values of PA = 95.2 ± 2.4 • and PA gas = 99.0 ± 16.2 • indicate that both stellar and gas rotation are aligned.Altogether, Mrk 1044's kinematics are quiescent.The smooth radial behaviour of the kinematic parameters indicates that Mrk 1044 does not currently experience major interactions with a companion or external disturbances.

Discussion
Despite its super-Eddington accretion rate, Mrk 1044 shows an unperturbed, rotation dominated velocity field.Moreover, there is no close companion visible in continuum light that could trigger a nuclear starbursts or AGN activity.Alternatively, intrinsic processes could be responsible for the enhanced star formation activity near the nucleus.

Star Formation as potential driver for the AGN fuelling
Contradicting results regarding the connection between star formation and AGN have been reported in the literature.While some authors claim a correlation between global SFR and AGN luminosity (e.g.Bonfield et al. 2011;Rosario et al. 2012;Chen et al. 2013;Xu et al. 2015;Lanzuisi et al. 2017;Stemo et al. 2020) others find only weak or no correlation (Baum et al. 2010;Rosario et al. 2013;Azadi et al. 2015;Stanley et al. 2015Stanley et al. , 2017)).This discrepancy can partially be explained by the different timescales of SF and AGN variability and the fact that the time-delay to trigger SF is typically longer than the AGN duty cycle (Silk 2013;Zubovas et al. 2013;Harrison et al. 2021).In addition, different spatial scales over which galaxy integrated properties are measured (∼10 kpc) compared to the BH accretion disk (<1 pc) add further scatter to the correlation.
The importance of the distance on which AGN affect their hosts has been pointed out by Volonteri et al. (2015) who used a suite of hydrodynamical simulations to compare the BHAR with the galaxy SFR on different spatial scales.We show their results together with Mrk 1044's position in Fig. 13.In their simulations, Volonteri et al. (2015) describe galaxies as ioslated ('stochastic' phase), during a 'merger' that lasts ∼ 0.2 − 0.3 Gyr or in the 'remnant' phase.Their simulations predict that, in contrast to works from Mullaney et al. (2012) and Chen et al. (2013), the galaxy-wide SFR is uncorrelated with BHAR.Instead, only the nuclear (∼100 pc) SFR correlates with BHAR.Mrk 1044 is in-line with this prediction.We have shown that on galaxy scales, the star forming CNE contributes a significant part of the total SFR which pushes Mrk 1044's into the merger regime.On nuclear scales, however, the CNE in Mrk 1044 follows the SFR 100 pc -BHAR correlation.Volonteri et al. (2015) further find that nuclear SFR and BHAR variations have similar amplitudes and both vary on similar timescales.
These results are supported by starbursts observed in the inner regions of nearby Seyfert 1s where the starburst intensities increase with AGN luminosity and Eddington ratio (Imanishi & Wada 2004;Davies et al. 2007;Watabe et al. 2008;Diamond-Stanic & Rieke 2012).In this regard, Mrk 1044 represents an extreme example where the high star formation rate surface density coincides with a phase of strong nuclear activity.The critical question is whether both processes are physically connected?SFR vs. BHAR Since a substantial fraction of Mrk 1044's galaxy-wide star formation is concentrated in a circumnuclear ellipse, we want to find out whether it can explain the enhanced BHAR.Assuming that in Mrk 1044's star forming CNE massive stars quickly return their material to the ISM through stellar winds and SNe, we can compare the mass fluxes from star formation and BH accretion.
Hα nebular emission arises from the recombination of gas ionized by O-and early-type B-stars with M 20 M (Peters et al. 2017).Due to their short lifetime, the presence of Hα emission implies ongoing star formation within the last 5 Myr.We assume that the lifetime of the massive stars is short compared to the total duration of the present day star formation activity, i.e. they instantaneously return their material to the ISM except for a remnant mass M r .We adopt the metallicity dependent compact remnant mass function from Fryer et al. (2012, eq. 5,7 and 9).Due to the vast abundance of molecular hydrogen in the host galaxy (4 × 10 8 M , Bertram et al. 2007), we assume a constant SFR that equals the present-day value of 0.19 M yr −1 estimated in Sect.3.2.2.We can express the rate of enriched material that Under the above assumptions the present-day net mass flux from stellar winds and SNe in the CNE exceeds the current BHAR by a factor of 1.5.Since we have not made a geometric assumption, comparing both quantities requires that all stellar ejecta remain within the CNE and are channeled towards the accretion disk which may not be the case in reality.We also note that choosing a higher η that is typical of high accretion rates (Ohsuga et al. 2005;McKinney et al. 2015) reduces the BHAR.Nonetheless, our simple model supports the scenario where the star formation in the circumnuclear ellipse drives the high accretion rate of Mrk 1044.

Lifetime of the CNE star formation
The abundance of dense molecular gas from which the stars are formed correlates with the BHAR (Izumi et al. 2016;Shangguan et al. 2020).Theoretical models suggests a physical connection between molecular gas abundance and BHAR on scales of a few 100 pc down to the accretion disc, independent from galaxy-scale processes (Umemura et al. 1997;Thompson et al. 2005;Kumar & Johnson 2010).We can estimate the maximum duration of the BH growth at the current rate by assuming a closed box model.The mass evolution of the gaseous disk with the cold gas mass M H 2 is given by Assuming that the cold gas content M H 2 = 4 × 10 8 M (Bertram et al. 2007) will be entirely accreted onto the BH, the maximum duration of the gas accretion is With this simplistic estimate, the central BH of Mrk 1044 will have consumed its entire pristine cold gas reservoir after 1.7 Gyr with a final mass M • = 4 × 10 8 that equals the presentday HI mass.However, we already see a spatially resolved outflowing component in the ionized gas kinematics (see Fig. 7).Neglecting in-and outflows from the CNE might therefore not be a robust assumption.Using the BH mass -extended narrow line region size relation from (Husemann et al. 2022), the expected timescale for an AGN-driven outflow to propagate through the CNE is ∼ 10 4.5 yr, well below the time that the BH requires for significant mass growth.Furthermore, it is questionable whether the strong AGN phase will affect the host galaxy star formation activity through AGN feedback.Due to the proximity of the CNE near the active nucleus, the energy input irrespective of the coupling mechanism would first affect the CNE star formation before it would propagate through the entire host galaxy.The presence of ongoing star formation so close to the nucleus means that Mrk 1044's BH fuelling is a recent event.A detailed discussion on Mrk 1044's AGN-driven outflow and its future impact on the host galaxy SFR will be presented in Winkel et al. (in prep.)

Chemical Enrichment of the CNE
Using UV spectra obtained with HST Imaging Spectrometer Fields et al. (2005a) have identified two outflowing absorbing systems that escape from Mrk 1044's center.From the column densities of O iv, N v and HI measured with the Far Ultraviolet Spectroscopic Explorer Fields et al. (2005b) estimated their metallicity to be at least five times solar.Such high metal abundances are difficult to explain in a galaxy evolution context.Fig. 14 shows that, compared to the overall star forming galaxy population Mrk 1044's star forming CNE has a gas-phase metallicity that is higher-than-average by 0.11 dex.This result is in agreement with the unperturbed velocity field, as we do not see signs of an external gas reservoir that could fuel Mrk 1044's center with fresh metal-poor gas.Instead, the high SFR suggests that Mrk 1044 is currently diverging from the mass-metallicity relation.The CNE size and its enrichment state are well in agreement with the TIMER studies of nuclear discs and nuclear rings (Bittner et al. 2020).With an approximate bar size of ∼ 5 kpc bar (see Fig. 1), Mrk 1044 fits well into the correlation between the nuclear discs size and bar size (Gadotti et al. 2020), which suggests a secular bar-driven formation of the CNE.

On the black hole feeding mechanism
We have shown above that the mass ejection rate from dying stars in Mrk 1044's circumnuclear star formation region is comparable to the current BHAR.Furthermore, the high metal abundances in the CNE and the accretion disk suggests that the gas has the same origin.To provide both the ongoing SF and the AGN with material, there must be efficient gas channelling from galaxy scales towards the center.In the following we discuss the secular host galaxy processes that could trigger the radial gas migration.
Since the ionized gas co-rotates with the galaxy disk, the gas transport towards the center must be accompanied by a significant loss of angular momentum.One ingredient to achieve this is to increase the turbulence in the ISM which allows the gas to dissipate its angular momentum.There have been several disk-internal processes suggested that induce turbulence.For example, SN explosions and stellar winds deposit energy into to a pressure-supported circumnuclear gas disk.The higher turbulence increases the effective kinetic viscosity, i.e. enables a more efficient angular momentum transfer (Kawakatu & Wada 2008;Wutschik et al. 2013).For Mrk 1044, the high SFR in the circumnuclear ellipse could be a reflection of this process.
Another process has been suggested for collisional disks where the cold gas could radially migrate towards the galaxy center via chaotic cold accretion (CCA, e.g.Gaspari et al. 2019).The circular motion is perturbed by recurrent collisions between cold gas cells which reduces their angular momentum and thus enable radial mass transport towards the center.In this way CCA can rapidly boost inflow rates from the galaxy meso-(∼ 100 pc) down to the micro-scale (<1 pc, Gaspari et al. 2017).Independent of how it is injected into the rotating disk, turbulence allows the gas to dissipate its angular momentum and radially migrate towards the galaxy center.The enhanced star formation seen in Mrk 1044's CNE could either be a result of gravitational instabilities within the self-gravitating eccentric disk or vice-versa.
Other astrophysical processes may further boost the angular momentum transport.The strong eccentricity of Mrk 1044's CNE suggests that its formation is driven by a galaxy scale process.Furthermore, Mrk 1044's host morphology is not that of a bona fide barred spiral.As visible in the left panels of Fig. 1, the inner ∼500 pc are dominated by dusty spiral arms and are slightly elongated from the south-east to the north-east.The spiral arms are surrounded by a tenuous bar with a semi-major axis of ∼2 kpc that in contrast stretches in the east-west direction.On the largest scales up to ∼15 kpc, the galaxy outskirts exhibit once more a spiral structure that is elongated in the north-south direction.Due to its peculiar host morphology, the gas migration towards Mrk 1044's central engine could be driven by angular momentum transport through gravitational torques (Shlosman et al. 1990).The minor axis of the CNE coincides with the orientation of the bar and is therefore likely a resonant ring that typically forms at the inner Lindblad resonances (ILR).At that location of the ILR, the gas-flow towards the center is slowed down and it accumulates in a reservoir (Combes 1996;Comerón et al. 2010).This is is supported by the enhanced SF within Mrk 1044's circumnuclear ellipse which is likely the reflection of a gas repository.
At smaller distances than 100 pc from the center, large-scale torques become inefficient and disk instabilities are required to transport gas down to scales of the accretion disk ("bars within bars", Shlosman et al. 1989).Hydrodynamic simulations from Hopkins & Quataert (2010) show that this process would explain the high star formation rate through the formation of circumnuclear molecular gas clumps.This means that gas accretion and SF are self-regulating and can coexist in an inhomogeneous, gravitationally unstable disk (Vollmer et al. 2018;Schartmann et al. 2018).In this scenario, the gas may only be partially converted into stars since the advection timescale is shorter than the star formation timescale (Thompson et al. 2005).Therefore, the radial gas transport allows for efficient channelling to the accretion disk, which would explain Mrk 1044's high accretion rate.As the motion of the stellar component is driven by the gravitational potential alone, the angular momentum transport only affects the gas component.
For Mrk 1044, we see a possible signature of this process in the spiral trails in the ionized gas velocity field (Fig. 11).The majority of galaxy spiral arms are trailing (e.g.Pasha & Smirnov 1982;Puerari & Dottori 1992) which can be explained with swing amplification of density waves that direct inwards (Toomre 1981;Athanassoula et al. 2009).Thus, the deviation from circular rotation likely traces the radial gas migration towards Mrk 1044's galaxy center.On the other hand, the ionized gas contributes a small fraction of the supply for the SMBH compared to the cold molecular gas (e.g.Gaspari et al. 2017).While the ionized gas is accelerated by momentum injection through radiation and winds form SNe, the dense inner regions of molecular clouds might be shielded from both and its motion remain unperturbed.However, in a CCA scenario, the molecular gas is expected to inherit part of the galaxy-scale kinematics via the top-down condensation, hence leading to tighter phase correlations (Gaspari et al. 2018).Overall, we require high angular resolution of the molecular gas kinematics in Mrk 1044's central region to understand the detailed gas feeding towards the accretion disk, in particular to unveil the cooler phases.

Summary and conclusions
In this work we have presented new MUSE NFM-AO observations of the nearby unobscured NLS1 galaxy Mrk 1044.The deblending of the AGN from the host emission enabled us to extract the host galaxy properties with unprecedented resolution.We have mapped the host galaxy stellar velocity field and the emission line properties in Mrk 1044's central kpc.Furthermore, we have used a kinematic model to identify structures that deviate from a quiescent rotating disk.Our key results are the following.
-A large fraction (27%) of the galaxy-wide star formation is concentrated in a compact star forming circumnuclear ellipse (CNE) with Σ SFR = 0.16 ± 0.04 M yr −1 kpc −1 , a semi-minor axis of 306 ± 43 pc and an axis ratio of 2.45.Within the CNE the optical emission line ratios are entirely consistent with excitation by star formation.Surprisingly, we do not find signs of extended photo-ionized gas emission, despite Mrk 1044's high accretion rate.-From 30 kpc down to ∼100 pc both stellar and ionized gas follow the quiescent kinematic pattern of a rotating thin disk.
In the velocity fields, there are no signs of a major interactions with a companion.-Within the compact CNE, the ionized gas exhibits clumpy streams with ∼ 100 pc extent that deviate from rotation with the disk, potentially contributing to the SMBH feeding.Their projected line-of-sight velocity is of the order of 10 km s −1 .-The gas-phase metallicity within the star forming CNE is above the expected value from the mass-metallicity relation.
-The estimated mass ejection rate from massive stars exceeds the current BHAR.
Our results suggest that the circumnuclear star formation is connected to the high accretion rate of Mrk 1044's central black hole.It remains unclear, however, if the star formation drives Mrk 1044's black hole accretion or if both share an underlying process that triggers them at the same time.Very long baseline interferometric observations of Mrk 1044's nuclear region have been conducted recently with the Atacama Large Millimeter/submillimeter Array (ALMA).The information on the molecular gas abundance and distribution will help to understand the fuelling of Mrk 1044's central engine.However, to constrain the initial process driving the rapid black hole growth in NLS1, a comparative study of nearby BLS1s and NLS1s with high spatial resolution is required.

Fig. 1 .
Fig. 1.Comparison of the FOV between Mrk 1044's observations with HST WFC3/UVIS image in the F547M filter (left), VLT MUSE WFM (center) and VLT MUSE NFM-AO (right).The two images from MUSE are white-light images which were created by integration over the wavelength axis.Both WFC3/UVIS and MUSE WFM have a FOV that covers the host galaxy beyond its effective radius of its disk r e = 7.14 ± 0.04 kpc (Wang et al. 2014), whereas MUSE NFM-AO only captures the innermost 3.1×3.2kpc.

Fig. 2 .
Fig. 2. Wavelength-dependent AGN position on the detector plane before (transparent) and after (opaque) correcting for differential atmospheric refraction.The wavelength range around 5890 Å is blocked by a dichroic in the optical path to avoid contamination and saturation of the detector by the strong sodium emission from the AO laser.

Fig. 3 .
Fig. 3. Fitting the MUSE NFM-AO PSF.From left to right the panels show the empirical PSF for the broad Hβ line extracted as described in Sect.2.4.1, the corresponding PSFAO19 model and the residual flux map.The rightmost panel shows the cross section of the PSF at the QSO position and fixed δ.While the turbulent halo of the AO-induced PSF is well reproduced by the PSFAO19 model, the systematic errors near the center are significant compared to the relatively faint host signal.

Fig. 4 .
Fig. 4. Construction of hybrid PSF from the MUSE data cube.The left panels show the continuum (blue) and broad line (red) spectral regions from which the broad line intensity is extracted around Hβ, Hα and O iλ8446+Ca iiλ8498 (from top to bottom).The resulting empirical PSFs with logarithmic intensity scaling are shown in the 2 nd column.The corresponding best-fit PSFAO19 models are shown in the 3 rd column.The fourth column shows the hybrid PSFs which we generate by replacing the empirical PSFs with the modeled PSFs beyond a the radius where the S/N drops below a threshold value (dashed red circle).

Fig. 5 .
Fig. 5. Result of the iterative AGN-host deblending for Mrk 1044 in spectral (left) and spatial (right) dimensions.The panel on the left shows spectrum the original blended spectrum (black) within a 3 aperture centered on the QSO position.The AGN emission (red) contains the broad lines, whereas in the the host galaxy spectrum (blue) only narrow emission lines are present.The bright point-like AGN emission distributes over the entire FOV of the data cube which is captured by the Hα narrow-band images on the right.Both original (left) and deblended host (right) images share the same colour-scaling, which demonstrates the successful deblending.

Fig. 6 .
Fig. 6.Fitting of an example host spectrum with PyParadise.The panel on the left shows an Hα narrow-band image of the AGN-subtracted host emission.An example 8×8-binned spaxel is highlighted with a white square north of the nucleus.The corresponding host spectrum extracted from this aperture is shown in the right panel as a black line in the Hβ (left) and Hα window (right).Furthermore, we show the best-fit spectrum of the stellar continuum (blue) and the ionized gas contribution (red) which we obtained with PyParadise.The residuals in the bottom panels show that the model reproduces the spectrum within the 3σ confidence region.

Fig. 7 .
Fig. 7. Mapping the ionized gas properties of Mrk 1044.The top left panel shows the flux contours of the HST WFC3/UVIS image which capture the spiral dust absorption features.The bottom left panel contains the outer contours where Hα surface brightness exceeds 3 × 10 −18 erg s −1 cm −2 which indicates the core region of enhanced narrow Hα emission.Here and in the following we show both contours as a reference for the size of the structures identified.The surface brightness maps for the prominent narrow emission lines after subtracting the AGN emission are shown in (a) -(d).The maps are spatially co-added by 2×2 pixels and 8×8 pixels for regions where S/N of the emission line is < 3. The LOS velocity (e) and velocity dispersion (f) of Mrk 1044's ionized gas appear quiescent.Panel (g) shows the BPT diagram for Mrk 1044's central region for 8×8 spatially co-added spaxels.The grey star shows Mrk 1044's position based on its AGN-contaminated 3 aperture spectrum.The empirical demarcation lines from Kauffmann et al. 2003 (continuous), Kewley et al. 2001 (dashed) and Cid Fernandes et al. 2010 (dotted) define classifications into star forming (SF), composite (inter), Low-Ionization Nuclear Emission Region (LINER) and AGN ionized regions.Already the AGN-blended 3 spectrum is located in the SF-regime.For the AGN-subtracted cube, all spaxels as well as the 0 .5-integrated host component are consistent with excitation by star formation.(h) shows the spatial distribution of the excitation mechanism.Mrk 1044 is ionized by star formation, even in the immediate vicinity of its nucleus.

Fig. 8 .
Fig. 8. Mapping the star formation rate in the vicinity Mrk 1044's nucleus.From left to right the panels show (1) the HST/UVIS contours, (2) Mrk 1044's extinction map, (3) the star formation rate surface density and (4) the gas-phase metallicity.The contour of the Hα ring is overplotted in white.While both extinction and star formation are concentrated in a circumnuclear ring, the chemical composition of the ionized gas is homogeneous within the inner 1 kpc.

Fig. 9 .
Fig. 9. Difference between Mrk 1044's gas and stellar velocity field from MUSE WFM (top) and NFM-AO (bottom).Except for the high velocity bloom that stretches ∼5 -15 north from the center, the gas rotates faster than the gravitationally driven motion of the stellar component.

Fig. 11 .
Fig. 11.Zoom-in onto the clumpy structures detected in the residual gas velocity ∆v r after subtracting the rotating thin disk model.The contours highlight spaxels of the smoothed map (see text) and correspond to an approaching (blue) and receding (red) clump respectively.The panels on the right show the rest-frame spectra of the host in the Hα+[N ii] region, integrated over the respective patch.The emission lines have a symmetric shape and are well reproduced by the fit, indicating a genuine shift of the velocity residuals ∆v r .

Fig. 12 .
Fig. 12. Radial dependence of the rotating disk model parameters.Both model parameters of the galaxy stellar (blue) and ionized gas (blue) velocity field are shown.The data points acquired from the MUSE NFM-AO observations are shown as stars and have errors smaller than the size of the symbol.The data points from MUSE WFM are shown as errorbars.From 0 .8 to 25 the PA of the disk-like rotation is constant and matches between ionized gas and stars.For both, the rotational velocity increases from the center up to ∼6 before it flattens out.

Fig. 13 .
Fig. 13.Comparison of Mrk 1044's BHAR and SFR on different spatial scales with theoretical predictions.The contours are taken from Volonteri et al. (2015) and describe the BHAR during merger (blue), stochastic (red) and remnant (yellow) phase respectively.The lines represent the relations between galaxy-scale SFR and BHAR from Mullaney et al. 2012 (mass-selected galaxies, dashed line) and Chen et al. 2013 (star forming galaxies, dashed).The dash-dotted line marks the line that separates AGN and SF dominated regions.On galaxy scales Mrk 1044 has a high star formation rate that compares to that of merging galaxies.In the center, the high SFR of Mrk 1044's CNE follows the predicted correlation with BHAR.
fading stars return into the ISM as ṀSNe = m max m min [m − M r ] × SFR × ξ(m)dm (11) where ξ(m) is the initial mass function (IMF) and SFR × ξ(m) corresponds to the birth rate of stars with mass m.We chose a minimum mass of m min = 8 M for stars that enrich the ISM on the relevant timescale of ∼ 10 Myr.Further, we assume a default Starburst99 IMF (Leitherer et al. 1999) from 0.1 M to 100 M and metal-rich gas (Z/Z =2) to estimate the mass ejection rate from stellar winds with Eq. 11 ṀSNe = 3.6 × 10 −2 M yr −1 = (1.4 ± 0.2) × 10 44 erg s −1 reported in Husemann et al. (2022) and a radiative efficiency of η = 0.1, we estimate the black hole accretion rate (BHAR) Ṁ• = L bol ηc 2 = 2.5 × 10 −2 M yr −1 (13)