The multiphase gas structure and kinematics in the circumnuclear region of NGC 5728

TL;DR

This study combines multi-wavelength observations to analyze the complex gas structures and kinematics around the AGN in NGC 5728, revealing inflows, outflows, and star formation within the circumnuclear region.

Contribution

It provides a detailed multi-scale analysis of gas dynamics and AGN feedback in NGC 5728, highlighting the coexistence of inflows, outflows, and star formation in the circumnuclear environment.

Findings

Gas inflow rate of 1 M$_{\ m sun}$ yr$^{-1}$ at the ring

AGN-driven ionized gas outflow reaches 250 pc with 0.08 M$_{\rm sun}$ yr$^{-1}$

The circumnuclear region largely remains structurally intact despite AGN activity

Abstract

We report on our combined analysis of HST, VLT/MUSE, VLT/SINFONI, and ALMA observations of the local Seyfert 2 galaxy, NGC 5728 to investigate in detail the feeding and feedback of the AGN. The datasets simultaneously probe the morphology, excitation, and kinematics of the stars, ionized gas, and molecular gas over a large range of spatial scales (10 pc--10 kpc). NGC 5728 contains a large stellar bar which is driving gas along prominent dust lanes to the inner 1 kpc where the gas settles into a circumnuclear ring. The ring is strongly star forming and contains a substantial population of young stars as indicated by the lowered stellar velocity dispersion and gas excitation consistent with HII regions. We model the kinematics of the ring using the velocity field of the CO (2--1) emission and stars and find it is consistent with a rotating disk. The outer regions of the disk, where the dust lanes meet the ring, show signatures of inflow at a rate of 1 M$_{\sun}$ yr$^{-1}$. Inside the ring, we observe three molecular gas components corresponding to the circular rotation of the outer ring, a warped disk, and the nuclear stellar bar. The AGN is driving an ionized gas outflow that reaches a radius of 250 pc with a mass outflow rate of 0.08 M$_{\sun}$ yr$^{-1}$ consistent with its luminosity and scaling relations from previous studies. While we observe distinct holes in CO emission which could be signs of molecular gas removal, we find that largely the AGN is not disrupting the structure of the circumnuclear region.

Figures

Equations

$$R^i(A_{M/I})=\frac{R_{\mathrm{IB}}^i+0.568A_{M/I}{R}_{\mathrm{MB}}^i}{1+0.568A_{M/I}}$$

$$L_{\mathrm{CO}}^{\prime }=3.25\times 10^7\frac{S_{\mathrm{CO}}{R}_{12}{D}_L^2}{(1+z){\nu }_{\mathrm{rest}}^2}\mathrm{K}\mathrm{km}{\mathrm{s}}^{-2}\mathrm{p}{\mathrm{c}}^2$$

$$M_{\scriptscriptstyle\rm HII}=3\times 10^{19}\Biggl{(}\frac{F_{\scriptscriptstyle\rm Br\gamma}}{\rm erg\,cm^{-2}\,s^{-1}}\Biggr{)}\Biggl{(}\frac{D}{\rm Mpc}\Biggr{)}^{2}\Biggl{(}\frac{n_{e}}{\rm cm^{-3}}\Biggr{)}^{-1}$$

$${\dot{M}}_{\mathrm{out}}=\left(\frac{M_{\mathrm{HII}}{v}_{\mathrm{out}}}{\delta r}\right),$$

$${\dot{E}}_{\mathrm{out}}=\frac{1}{2}{\dot{M}}_{\mathrm{out}}{v}_{\mathrm{out}}^2,$$

$${\dot{p}}_{\mathrm{out}}={\dot{M}}_{\mathrm{out}}{v}_{\mathrm{out}},$$

$${\dot{M}}_{\mathrm{acc}}=\frac{L_{\mathrm{Bol}}}{c^2\eta },$$

Full text

††pubyear: 2019††pagerange: The multiphase gas structure and kinematics in the circumnuclear region of NGC 5728–B

1 Introduction

Beyond being the some of the most energetic objects in the Universe, active galactic nuclei (AGN) are thought to play an important role in the evolution of their host galaxies. In particular, large scale cosmological simulations require feedback from AGN to reproduce the galaxy population we observe today (e.g. Springel et al., 2005; Bower et al., 2006; Croton et al., 2006; Nelson et al., 2019). More complete knowledge of the processes that fuel AGN and the mechanisms through which they provide feedback is then necessary to understand the path galaxies take from being gas rich and star-forming to gas poor and quiescent.

To obtain this understanding however requires dedicated observations that cover a large range of wavelengths and spatial scales to reveal the dominant processes. AGN only occupy the very centres of galaxies but both streaming gas fuelling the AGN as well as outflowing gas from AGN driven winds or jets can reach at least several kpc away (e.g. Rupke & Veilleux, 2011; Cicone et al., 2014; Fischer et al., 2013; Bae et al., 2017; Baron et al., 2018; Fischer et al., 2018; Herrera-Camus et al., 2018; Kang & Woo, 2018; Mingozzi et al., 2018; Bischetti et al., 2019). Furthermore, this gas does not exist in a single phase but rather consists of a mixture of different phases from cold molecular gas to warm ionized gas to hot X-ray emitting gas. Thus to obtain a full accounting of the direct fuelling and feedback of the AGN will need multiple continuum and emission line tracers to probe each phase of the gas.

With the rise of integral field unit (IFU) instruments, large radio interferometers and adaptive optics (AO), it is now possible to study simultaneously the distribution and kinematics of both ionized and molecular gas at the same spatial scales in an efficient manner. For nearby galaxies, this translates to a wealth of information that can span tens of pc to hundreds of kpc, and the kinematics of the gas can further be compared to the stellar kinematics through the analysis of stellar absorption lines to reveal non-circular motions such as radial inflow or outflow. Indeed there has been an increased effort in combining mutliwavelength data to obtain a full view of the conditions and dynamics of gas around AGN. Recent examples include NGC 5643 (Alonso-Herrero et al., 2018), ESO 428-G14 (Feruglio et al., 2019), ESO 578-G009 (Husemann et al., 2019), NGC 3393 (Finlez et al., 2018), NGC 1566 (Slater et al., 2019), NGC 2110 (Rosario et al., 2019), and zC400528 (Herrera-Camus et al., 2018). These works primarily paired ALMA observations of CO emission with an optical or NIR IFU observation to measure inflow and outflow velocities and gas masses, gas excitation, and bar driven instabilities.

The importance of multiphase observations, in particular for studying AGN feedback, was outlined in Fiore et al. (2017) who compiled and reanalysed data probing AGN driven outflows in the molecular, ionized and X-ray emitting phase. Their findings revealed strong correlations between the AGN luminosity and outflow velocity, mass outflow rate, and wind momentum for all phases of the gas. However the scaling relations differ depending on the phase resulting in a changing fraction of outflowing gas in one phase that is dependent on the strength of the AGN. The Fiore et al. (2017) scaling relations also do not extend to lower luminosities that are more representative of the larger population of AGN. It remains to be seen whether they are valid for Seyfert-like AGN with recent work using dust as a tracer of AGN outflows suggests the relations break down (Baron & Netzer, 2019).

Assessing an AGN’s impact on its host galaxy requires measurements of important physical properties of the outflowing gas including the radial extent and the electron density. Without spatially resolved data and lines that accurately trace the electron density, estimates of the mass outflow rate and energetics can be off by orders of magnitude. Revalski et al. (2018a), using high spatial resolution HST imaging and long slit spectra, showed that while global estimates of outflow properties overall agree with spatially resolved estimates, they come with large uncertainties and highly depend on the assumed geometry and density of the system. Only with precise knowledge of the radial extent and gas density will globally averaged mass outflow rates, kinetic luminosities, and momentum rates come close to those obtained from spatially resolved data.

A parallel and equally important focus of studies of nearby AGN is not how gas is being removed but rather how gas is fuelling the AGN. Many mechanisms have been proposed to extract angular momentum and drive gas to the central supermassive black hole (SMBH) but so far no unique one has been found to explain all observations of AGN. Instead a range of mechanisms can likely fuel AGN that depends on the luminosity and morphology of the host galaxy. For low to moderate luminosity AGN, recent work suggests that gravitational torques produced by non-axisymmetric potentials are efficiently bringing gas down to the inner kiloparsec (e.g. Combes, 2003; García-Burillo et al., 2005; García-Burillo & Combes, 2012, ; however also see Sormani et al. (2018)). Hydrodynamic simulations show that gas flowing along a bar settles into a nuclear ring, spatially located in or near the Inner Lindblad Resonance (ILR), but further inflow inner to the ring is small without subsequent dynamical instabilities such as another nested nuclear bar (Combes & Gerin, 1985; Piner et al., 1995; Regan & Teuben, 2003; Maciejewski et al., 2002; Maciejewski, 2004). Instead of nuclear rings, nuclear spirals can also form and propagate down to the SMBH given a higher gas sound speed (Maciejewski et al., 2002; Maciejewski, 2004) and indeed nuclear spirals have been found in many galaxies hosting an AGN (Pogge & Martini, 2002; Martini et al., 2003; Prieto et al., 2005; Davies et al., 2009; Davies et al., 2014). Since cold molecular gas is the fuel for AGN, it makes sense to turn to observations of CO line emission, a primary tracer for molecular gas, to detect and measure the structure and instabilities bringing gas to the SMBH.

1.1 NGC 5728

NGC 5728 is a nearby SAB(r)a galaxy ($D=39$ Mpc, 1″= 190 pc), with a Compton thick, Seyfert 2 nucleus (Veron, 1981; Phillips et al., 1983). The bolometric luminosity of the AGN is $1.40\times 10^{44}$ ergs s*-1* where we have converted the absorption corrected 14–195 keV X-ray luminosity from the Swift Burst Alert Telescope, $L_{X}=2.15\times 10^{43}$ ergs s*-1* (Ricci et al., 2017), into $L_{\rm Bol}$ using the relation from Winter et al. (2012).

Its primary characteristics include a large stellar bar ($R\approx 11$ kpc, P.A. = 33°; Schommer et al., 1988; Prada & Gutiérrez, 1999) which is surrounded by a ring of young stars, a circumnuclear star forming ring ($R\approx 800$ pc; Schommer et al., 1988; Wilson et al., 1993), and extended ionization cones ($R\approx 1.5$ kpc, P.A. =118°; Schommer et al., 1988; Arribas & Mediavilla, 1993; Wilson et al., 1993; Mediavilla & Arribas, 1995). Rubin (1980) was the first to study the ionized gas kinematics finding that while at large radii the gas displayed normal rotation, the central regions showed strong non-circular motion as well as double peaked emission lines. Over the years much discussion ensued over the physical explanation of the non-circular motion and double peaked emission lines with radial outflow (Rubin, 1980; Wagner & Appenzeller, 1988; Arribas & Mediavilla, 1993; Mediavilla & Arribas, 1995) and radial inflow due to the bar (Schommer et al., 1988) the primary choices. van Gorkom et al. (1983) and Schommer et al. (1988) also found extended 6 and 20 cm radio emission spatially coincident with the ionization cones suggesting the presence of a jet. Further complicating interpretations was the finding of a strong isophotal twist ($\Delta$P.A. $\approx$ 60°) in the central $\approx 5$″ from NIR imaging that suggested the presence of a secondary nuclear bar.

This secondary nuclear stellar bar was thought to potentially be the driving mechanism to bring gas from the circumnuclear ring to the SMBH. However, Petitpas & Wilson (2002) found no significant nuclear molecular gas bar coincident with the stellar bar based on spatially resolved CO (1–0) data. Instead they note a double peaked morphology of the CO emission that straddles the centre but is perpendicular to the ionization cones and suggest that the NGC 5728 might have experienced a recent minor merger that could cause the off-center CO emission.

All of these intriguing components contained with NGC 5728 make it an ideal system to study the fuelling and feedback of AGN. In this Paper, we report on a combined analysis of recent Atacama Large Millimeter/submillimeter Array (ALMA), Very Large Telescope (VLT) MUSE, and VLT/SINFONI observations that jointly probe the morphology and kinematics of the stars, cold molecular gas, hot molecular gas, and warm ionized gas on both large and small spatial scales. Recently, Durré & Mould (2018) and Durré & Mould (2019) studied NGC 5728 using both the MUSE and SINFONI datasets. Much of our conclusions and interpretation are consistent, but we will note where we disagree and where their analysis enhances ours. Throughout this Paper, we assume a $\Lambda$CDM cosmology with $H_{0}=70$ km s*-1* Mpc*-1* and $\Omega_{\rm M}=0.3$.

2 Observations, Data Reduction, and Analysis Methods

2.1 SINFONI Observations

NGC 5728 was observed by VLT/SINFONI (Eisenhauer et al., 2003; Bonnet et al., 2004) as part of the Luminous Local AGN with Matched Analogues (LLAMA) program (Davies et al., 2015), a volume limited survey of nearby X-ray selected AGN with SINFONI and XSHOOTER. Within this program, all targets were observed by SINFONI using the H+K grating (1.4 – 2.5 µm) and AO mode providing integral field spectroscopy with spectral resolution $R\sim 1500$ over a roughly 3″$\times$3″ field of view (FOV). A standard Object-Sky-Object observation sequence was used with each integration lasting 300 seconds and the Object exposures dithered by 0.3″. Observing blocks for NGC 5728, under program ID 093.B-0057 (PI R. Davies), were carried out on the nights of February 23, March 06, May 09, June 04, and June 25, 2015 resulting in 25 Object integrations and 12 Sky integrations for a total on-source time of 125 minutes.

The raw SINFONI data were reduced using the custom reduction package spred (Abuter et al., 2006) which performs all the standard reduction steps needed for near-infrared spectra as well as additional routines necessary to reconstruct the data cube. We further applied the routines mxcor and skysub (Davies, 2007) to improve the subtraction of OH sky emission lines and lac3d, a 3D Laplacian Edge Detection algorithm based on lacosmic (van Dokkum, 2001), to identify and remove bad pixels and cosmic rays. Telluric correction and flux calibration were carried out using B-type stars. Finally, we also corrected for differential atmospheric refraction, a wavelength dependent spatial offset due to the atmosphere, using a custom procedure described in Lin et al. (2018).

The final SINFONI data cube covers a FOV of 3.5″$\times 3$.7″ with a pixel scale of 0.05″. With AO, we achieved a point spread function (PSF) full width at half maximum (FWHM) of 0.15″ based on Gaussian fitting of the unresolved continuum nuclei from the AGN in the LLAMA sample. This translates to a physical spatial resolution of 28 pc over a physical FOV of 660$\times$700 pc.

2.2 MUSE Observations

MUSE is a second generation VLT instrument consisting of 24 integral field units that span the wavelength range 4650–9300 Å (Bacon et al., 2010). We downloaded archival MUSE data of NGC 5728 from the MUSE-DEEP collection, an ESO Phase 3 data release that provides fully reduced, calibrated, and combined MUSE data for all targets with multiple observing blocks. The original observations were performed over two observing blocks on April 3 and June 3, 2016 for program ID 097.B-0640 (PI D. Gadotti), and were carried out in seeing limited Wide Field Mode that results in a 1′$\times$1′ FOV sampled with a pixel scale of 0.2″. The total exposure time for the combined MUSE cube is 79 minutes and Gaussian fitting of a foreground star within the FOV show a PSF FWHM of 0.62″. The MUSE data cube thus covers an 11.3 kpc square region with a spatial resolution of 177 pc at the distance of NGC 5728.

Figure 1 shows a three color gri image of NGC 5728 from the MUSE cube. As seen in the HST imaging in Figure 2, the MUSE observation covers a large portion of the large-scale bar which contains long dust filaments and a blue circumnuclear disk/ring. In the right panel, we highlight strong [OIII] and H$\alpha$ emission by using a simple narrow band filter for the blue and green colours of the image while the red colour is from the i-band.

2.3 ALMA Observations

ALMA observed NGC 5728 on May 15, 2016 as part of project 2015.1.00086.S (PI N. Nagar) for a total on-source time of 29 minutes in the C36-3 configuration resulting in an angular resolution of 0.56″$\times$0.49″(PA = -79.4°) in Band 6. This observation was part of an ALMA program to map CO (2–1) in the nuclear regions of 7 local AGN (Ramakrishnan et al., 2019). The spectral setup was designed to center one of the spectral windows on the redshifted CO (2–1) emission line ($\nu_{\rm rest}=230.538$ GHz) for NGC 5728.

We used the Common Astronomy Software Applications package (CASA; McMullin et al., 2007) to process the data. Within the tclean routine, we applied Briggs weighting with robust = 0.5 to reproduce the data cube at 2.54 km s*-1* velocity resolution and a peak residual (parameter threshold) of 2.5 mJy. All resulting data cubes were primary beam corrected and a 1.3 mm continuum image was extracted by combining and integrating the remaining three spectral windows without the CO (2–1) line.

2.4 HST Observations

The final dataset we primarily use for this work is archival multiband HST imaging. HST observed NGC 5728 in the F336W, F438W, F814W and F160W bands with the Wide Field Camera 3 instrument as part of program 13755 (PI J. Greene). We downloaded the standard pipeline reduced images from the Mikulski Archive for Space Telescopes111https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html.

Figure 2 shows the F814W (left) with a three colour zoom-in of the circumnuclear region that combined the F814W, F438W, and F336W images as red, green, and blue colours respectively. The larger scale F814W image also has overplotted the FOV’s of the ALMA (red circle) and MUSE (green box) datasets. Both ALMA and MUSE cover a large portion of the primary bar and the entire circumnuclear region. The FOV of the SINFONI cube (white box) is overplotted on the three colour zoom-in that shows we are probing the very inner regions of NGC 5728. The combination of all of these datasets allow us to consistently trace gas and structure from several kpcs down to tens of pc as well as the interaction between the host galaxy and the central AGN.

2.5 Astrometric Alignment and AGN Position

To ensure all of the data sets are registered to a common coordinate system, we implemented the following procedure. Luckily at least one bright foreground star lies within the FOV of every HST image and the MUSE cube. Starting with the F336W, F438W, and F814W HST images, we matched the bright stars to GAIA DR2 sources (Gaia Collaboration et al., 2018). To make corrections, we simply updated the WCS parameters by applying small pixel shifts until the positions matched between the images and the GAIA DR2 values. These corrected HST images were then used to correct both the NICMOS F160W image and the MUSE cube, again by applying pixel shifts until the locations of the stars in common matched.

For the SINFONI cube, no foreground stars appear in the FOV, therefore we first produced an H-band continuum image from the cube and used it to match distinct features between it and the HST F160W image. In particular, we aligned the brightest pixel in the SINFONI H-band continuum image with the brightest pixel in the HST F160 image and then checked to ensure other features such as the nuclear dust lane were aligned as well. This gives us an astrometric accuracy of $\approx 0.25$″(i.e. the angular resolution of HST at 1.6 µm) between SINFONI and the rest of the datasets. Finally, we did not apply any corrections to the ALMA data given the high astrometric accuracy inherent in radio interferometric observations.

For the position of the AGN, we used VLBI observations from The Megamaser Cosmology Project (C. Y. Kuo 2018, private communication). The coordinates of the AGN are RA = 14:42:23.872 and DEC = 17:15:11.016 with an accuracy of a few tens mas. The position is shown in the right panel of Figure 2 as a red cross and will be used as the reference point for all following figures. This position matches very well both the peak of 1.3 mm continuum emission from ALMA (see Figure 4) as well as the peak of the warm molecular gas traced by H2 (1–0) S(1) emission from SINFONI (see Figure 19).

2.6 Spectral Fitting

2.6.1 Continuum Fitting and Subtraction

The data reduction for the ALMA cube includes already a continuum subtraction however both the SINFONI and MUSE cube contain significant continuum signal. To remove the continuum as well as stellar absorption features, we use the Penalized Pixel-Fitting222http://www-astro.physics.ox.ac.uk/~mxc/software/#ppxf (pPXF) routine (Cappellari & Emsellem, 2004; Cappellari, 2017) together with the latest (v11.0) single stellar population (SSP) models based on the extended MILES stellar library (E-MILES; Vazdekis et al., 2016). We chose the models produced from the Padova 2000 stellar isochrones with a Chabrier initial mass function (Chabrier, 2003). The SSP models range from stellar population ages of 63 Myr–14.13 Gyr spaced logarithmically and 6 metallicities between -1.71 and 0.22.

For the MUSE cube, both to reduce the computation time and improve the signal-to-noise (S/N) ratio in the outer regions, we chose to utilize Voronoi binning before running pPXF. We used the method from Cappellari & Copin (2003) and the voronoi_2d_binning333http://www-astro.physics.ox.ac.uk/~mxc/software/#binning Python routine to bin the MUSE cube such that each bin contains a minimum S/N in the continuum of 20.

Each binned or single MUSE and SINFONI spaxel was fit with ppxf. For the MUSE spectra, we used the full range of metallicities and stellar population ages but we also incorporated regularization to better interpret the best fit metallicities and ages since the fitting of galaxy spectra is an ill-posed problem (Cappellari, 2017). For the SINFONI spectra, because they only cover the near-infrared and are relatively insensitive to young stellar population ages, we only fit using models with ages of 0.1, 0.5, 1.0, 5.0, 10.0, and 14.1 Gyr. We also did not use regularization because our only concern for the SINFONI cube is the removal of the continuum and determination of the stellar kinematics. Because both the MUSE and SINFONI spectral resolutions are only moderate, we only fit for the stellar velocity and velocity dispersion. We also included an additive fourth degree polynomial in the fit to model any emission from the central AGN. To create the continuum subtracted MUSE cube, we renormalised the best fit binned spectra to the median continuum level of each individual spectrum before subtraction. Figure 3 shows the results of our pPXF fits for the MUSE cube.

2.6.2 Emission Line Fitting

To measure the bulk distribution and kinematics of the molecular and ionized gas in NGC 5728, we fit single Gaussians to the relevant emission line tracers of the gas. In particular, we fit the CO (2–1) line ($\nu_{\rm rest}=230.538$ GHz) in the ALMA cube; the H$\beta$, [OIII] doublet, [OI], [NII] doublet, H$\alpha$, and [SII] doublet lines in the MUSE cube; and the [FeII]$\lambda 1.64$ µm, [SiVI], H2 (1–0) S(1), and Br$\gamma$ lines in the SINFONI cube. This set of lines represent the strongest emission lines observed in nearby galaxies and AGN and provide diagnostic power into the state of gas in a range of phases from cold molecular gas (CO (2–1)) to warm molecular gas (H2 (1–0) S(1)) to hydrogen recombination lines (H$\beta$, H$\alpha$, Br$\gamma$) and forbidden lines ([OIII], [OI], [NII], [SII], [FeII], [SiVI]) from ionized gas in a range of ionization levels.

For each line, we cut out the relevant sections of the continuum subtracted spectra that contain the single line or multiple lines in the case of the [OIII] doublet, H$\alpha$ and [NII] blended region, and [SII] doublet. We fixed the velocity and velocity dispersion to be the same for all doublets and fixed the intensity ratio to the expected theoretical value of 3 for the [OIII] and [NII] doublets. For the fit of the [SiVI] line, we also included Gaussian components for the nearby H2 (1–0) S(3) and Br$\delta$ lines. As an estimate of the local noise around each line, we measured the root-mean-square (RMS) of the line-free regions. Using the RMS, we produced 100 simulated spectra for every emission line by randomly perturbing the original spectra assuming a Gaussian distribution with a standard deviation equal to the RMS. Each of the simulated spectra were fit in the same way as the original spectra and uncertainties on the integrated flux, velocity, and velocity dispersion were calculated from the standard deviation of the best fit values from the simulated spectra. Finally, all velocity dispersions derived from MUSE and SINFONI data were corrected for instrumental broadening. For lines from MUSE, we assumed an instrumental dispersion of 49 km s*-1* for lines in the [OIII] region and 71 km s*-1* for lines in the H$\alpha$ region based on the spectral resolution curve as a function of wavelength in the MUSE User Manual444https://www.eso.org/sci/facilities/paranal/instruments/muse/inst.html. For SINFONI, we assumed an instrumental dispersion of 65 km s*-1* based on the width of OH sky lines. Figures 17, 18, 19, and 20 show the results of our single Gaussian fitting for all of the emission lines studied in this paper over the full FOV for each instrument. Figure LABEL:fig:key_maps_zoom shows the results for CO (2–1), H$\alpha$, and [OIII] zoomed into the central 10″x 10″.

3 Results

3.1 Large Scales and the Primary Bar

We begin our investigation of NGC 5728 with the largest scales. As mentioned in Section 1.1, NGC 5728 is characterised on large scales by a primary bar surrounded by a ring of young stars. Continuum emission in general is elongated from the NE to the SW and the MUSE three-colour image (Figure 2) shows bluer emission encircling the inner, redder bar. These bluer emission regions in the large ring are also spatially coincident with clumps of H$\alpha$ emission indicating recent star formation.

In the top left panel of Figure 5, we show in both colour and contours the HST F160W image which traces emission from older stars. Again, on the largest scales we observe primarily the large bar with the P.A. of the bar plotted as the white line, however towards smaller radii the stellar distribution becomes more and more spheroidal with the major axis shifting to more N–S. This is also the P.A. for the stellar kinematic major axis as shown in Figure 3 as well as the photometric major axis for the main disk that the bar lives within (Schommer et al., 1988). Prada & Gutiérrez (1999) observed this isophotal twisting as well in large I-band imaging and attributed it to the increasing influence of the central bulge. This would further make sense given the increasing stellar velocity dispersion which rises to $\sim 170$ km s*-1* before dropping at the circumnuclear ring (Figure 3, top right panel) as well as having a relatively old mean age of $\sim 5$ Gyr compared to both the larger main ring and circumnuclear ring which show ages of $\sim 1$ Gyr (Figure 3, bottom left panel).

Roughly straight dust lanes can also be seen in Figure 1 and 2 running the length of the bar with the more prominent dust lane occurring in the northern half and indicating that the west side of the galaxy is the near side (Schommer et al., 1988). To further enhance and highlight the dust lanes, we show in the top right panel of Figure 5 a $V-H$ map constructed from the F438W and F160W HST images. Areas with high extinction would have large $V-H$ values and appear dark in the map. Indeed the northern dust lane reveals itself as a dark band extending from the NE to the circumnuclear region. We also plot on top of the $V-H$ map contours showing the location of CO (2–1) emission. Interestingly, the northern dust lane lacks strong CO emission and it is instead the southern dust lane that appears as a single spiral arm connecting to the circumnuclear ring. This morphology is consistent with the observations of CO (1–0) emission from Petitpas & Wilson (2002) however they attribute the southern arm as potentially due to a larger ring structure or evidence of a past merger event. Instead it is quite clear within our datasets that it is associated with a fainter dust lane on the far side of the galaxy. It is unclear however why CO (2–1) is missing in the northern dust lane. It is possible that molecular gas here is more diffusely distributed over larger scales that were resolved out of the interferometric observations.

As hydrodynamic simulations show (e.g. Athanassoula, 1992; Maciejewski et al., 2002), two straight dust lanes are expected to appear along the bar major axis due to the formation of principal shocks which compress the gas and dust. This occurs due to the gradual shifting of the gas orbits from the $x_{1}$ family (i.e. parallel to the bar major axis) to the $x_{2}$ family (i.e. perpendicular to the bar major axis). As gas crosses these shocks, it loses angular momentum resulting in inflow towards the centre (e.g. Athanassoula, 1992; Maciejewski et al., 2002). As we will show in the next section, we do observe kinematic signatures of inflow at the locations where the dust lanes connect to the central circumnuclear ring strongly suggesting that gas is being driven from the outer disk to the circumnuclear region via the bar. Assuming the gas in the bar is rotation dominated it is expected to settle into a ring (e.g Piner et al., 1995) or nuclear spiral (Maciejewski et al., 2002; Maciejewski, 2004).

3.2 The Circumnuclear Ring

As many others have before we observe a central circumnuclear ring within the central $\approx 1.3$ kpc in projected radius. Multiple tracers highlight this structure through both their emission and kinematics.

Both the HST and MUSE continuum imaging show a distinct ring of clumpy, blue emission indicative of young hot stars. This is confirmed in the bottom left panel of Figure 3 which shows the mean light-weighted stellar age. The circumnuclear ring outlined by blue continuum emission also has a distinctly younger mean light-weighted stellar age (1 Gyr vs 5 Gyr) compared to the surrounding region. This fits with the LLAMA finding that on average AGN show a significant young ($<30$ Myr) stellar population within the central few hundred parsecs compared to inactive galaxies based on stellar population modelling of high spectral resolution VLT/X-Shooter spectra (Burtscher et al. in preparation). Further, Riffel et al. (2009b) found a significant population of intermediate age stars in stellar population modelling of NGC 5728 using a NIR IRTF spectrum. The stars in the ring further show evidence for lower metallicity which could be an indication that the increased nuclear activity is due to a recent minor merger.

The circumnuclear ring is further revealed kinematically in the stars via strong circular rotation and a drop in velocity dispersion (Figure 3, top row). While on large scales the stellar velocity peaks at the edge of the FOV and with a P.A.$\sim$0°, there is a clear kinematically decoupled component with velocity peaks at $\sim$4″ (760 pc) radius and a P.A.$\sim 10$°. The dispersion map shows the ring as a ring of decreased velocity dispersions compared to the surrounding regions (70 km s*-1* vs. 170 km s*-1*) which spatially corresponds to the ring of younger mean stellar ages. The reduced velocity dispersions can be explained as the primary stellar population not yet thermalising out of their birth clouds and still moving with regular rotation. These so-called dispersion “drops” have routinely been used in long slit observations of nearby galaxies as markers of nuclear star formation (e.g. Emsellem et al., 2001) .

All of the stellar tracers suggest recent star formation so there should be large amounts of both molecular gas from which the stars were born out of and ionized gas from photoionization by the young stars. Indeed, Figure 6 shows that exactly where we find the ring of lowered velocity dispersions, we also observe a strong ring of CO (2–1) emission as well as a ring of H$\alpha$ clumps outlining the star forming HII regions.

The white ellipse plotted in Figure 6 was chosen to directly follow the ring of lowered stellar velocity dispersion. The parameters of the ellipse have a P.A. of 14°, a semimajor axis 822 pc and an axis ratio ($b/a$) of 0.8. If we make the assumption that the ring is intrinsically circular then the inclination is related to the axis ratio through $\cos(i)=b/a$. This leads to an estimate of the inclination of $\approx 37$° which is more face-on compared to the large scale photometrically determined inclination of 48–55° (Schommer et al., 1988; Davies et al., 2015). This could be due to our assumption of circularity, since its possible the stars in the ring are on elliptical $x2$ orbits associated with the large scale bar. These rings of young nuclear star formation associated with low velocity dispersions seem to be a common characteristic of local Seyfert galaxies as many studies using IFU data have found them within the central kpc (e.g. Riffel et al., 2010, 2011; Diniz et al., 2017; Dametto et al., 2019)

3.2.1 Gas Excitation

The MUSE and SINFONI IFU cubes allow for us to spatially resolve the distribution of emission line ratios and assess the dominant mechanism for gas excitation. Gas excitation mechanisms can be identified using popular line ratio diagrams with specific thresholds or demarcations that separate between the different mechanisms which include photoionization by a central AGN, ionization by UV photons from HII regions, and low-ionization nuclear emission line regions (LINER).

Figure 7 plots the three standard excitation diagrams (Veilleux & Osterbrock, 1987; Baldwin et al., 1981) for individual spaxels in NGC 5728. The circumnuclear region of NGC 5728 is dominated by gas photoionized by the AGN, as indicated by the yellow to red colours. The AGN-dominated region further follows the bicone shape of the strong [OIII] emitting region. HII regions instead follow a clear ring-like structure, encircling the bicone and tracing the ring-shaped emission seen in H$\alpha$, the three colour HST image, and the molecular gas distribution that is indicated by the white ellipse in the bottom row Figure 7. All of this definitively points to strong star formation occurring within the circumnuclear ring.

One interesting feature seen in the excitation diagram is a horizontal “spike” at high [OIII]/H$\beta$ ratios and pointing towards lower [NII]/H$\alpha$, [SII]/H$\alpha$, and [OI]/H$\alpha$ ratios. Within each of the diagnostic regions, we shaded the spaxels from light to dark based on the [NII]/H$\alpha$, [SII]/H$\alpha$, or [OI]/H$\alpha$ line ratio value. In the maps we can then see that the horizontal “spike,” indicated by yellow colours, is primarily occurring in a narrow region on the SE side of the circumnuclear ring and coinciding with where the ring intersects the AGN bicone.

This feature has been seen recently in other nearby AGN with spatially resolved line ratio maps. Mingozzi et al. (2018), using MUSE data for the MAGNUM survey, explained these line ratios as tracing the ratio of the solid angle subtended by matter bounded and ionization bounded gas clouds ($A_{\rm M/I}$) as first put forth in Binette et al. (1996). In short, matter bounded clouds are optically thin throughout the entire cloud to the ionizing radiation field while ionization bounded clouds are optically thick and absorb all of the incident ionizing radiation. If matter bounded clouds dominate the solid angle we would expect more emission from high ionization lines compared to low ionization ones and vice versa if ionization bounded clouds dominate.

Binette et al. (1996) using photoionization modelling produced predictions for the line ratios expected given an $A_{\rm M/I}$ with the following relation:

[TABLE]

where $R^{i}$ indicates a line ratio for a specific ion relative to H$\beta$ and $R_{\rm IB}^{i}$ and $R_{\rm MB}^{i}$ are constant ratios in the limiting cases for fully ionisation bounded clouds and matter bounded clouds respectively. We plot this sequence in the top panels of Figure 7 shown as a black line. Especially for the [NII]/H$\alpha$ plot, the increase of $A_{\rm M/I}$ perfectly explains the spike in the excitation diagram.

However, we can also explain this feature by an increase in the ionization parameter. (Baron & Netzer, 2019, ; BN19) developed a simple method for measuring the ionization parameter using the [OIII]/H$\beta$ and [NII]/H$\alpha$ line ratios (see Equation 2 of BN19). The main assumption is only that the gas is photoionized by an AGN. Therefore, we calculated the ionization parameter, $U$, as a function of radius towards the SW using a series of annular sections described Sec. B and shown in Fig. 23.

Fig. 8 shows the radial profile along the bicone of both [NII]/H$\alpha$ and $U$. Indeed, exactly where we observe a decrease of [NII]/H$\alpha$, the ionization parameter significantly increases. While both a change in $A_{\rm M/I}$ and a change in $U$ can explain the horizontal "spike" in the line ratio diagrams, we choose the simpler increase in $U$ model as our explanation. This would be consistent with an interpretation that this region is where the NLR bicone is intersecting the lower density circumnuclear ring since $U$ and gas density are inversely related.

3.2.2 Molecular Gas Mass

We estimate the total molecular gas mass in the circumnuclear region from the summed CO (2–1) spectrum resulting from integrating the ALMA cube within a circular aperture with a radius of 6.5″ (1.24 kpc). Because of the complex line profile, we simply integrated the spectrum between rest-frame velocities of $\pm 375$ km s*-1* to obtain a total CO (2–1) intensity ($S_{\rm CO}$) of 93.6 Jy km s*-1*. The CO luminosity is calculated following Solomon & Vanden Bout (2005):

[TABLE]

where $R_{12}$ is the conversion from CO (2–1) to CO (1–0) intensity, $D_{L}$ is the luminosity distance in Mpc, and $\nu_{rest}$ is the rest frequency of the CO (2–1) emission line in GHz. In this work we use $R_{12}=1.4$, the value found for nearby star forming galaxies (Sandstrom et al., 2013).

To calculate a molecular gas mass, we must assume a value for $\alpha_{\scriptscriptstyle\rm CO}$, the conversion factor from $L^{\prime}_{\rm CO}$ to $M_{H_{2}}$. For this work, we use $\alpha_{\scriptscriptstyle\rm CO}=1.1$ M*☉* pc*-2*/(K km s*-2*). This is the value found for the central regions of nearby galaxies, likely due to the increased pressure and/or turbulence (Sandstrom et al., 2013). Using this, we find a total molecular gas mass within the central $\sim$1 kpc of NGC 5728 to be 1.3$\times 10^{8}$ M*☉*.

This total mass matches well with that found in Rosario et al. (2018) who found a molecular gas mass of 108.53 M*☉*. Our value is 40% of theirs, however the Rosario et al. (2018) measurement is based on an JCMT observation with a beam size of 20.4″compared to the 6.5″aperture used here. This seems to be an indication that about 60% of the total molecular gas is on scales larger than 1 kpc given the largest angular scale of the ALMA observation was 5.34″. We tested whether convolving the ALMA cube to the JCMT beam and integrating over a larger aperture increased the measured molecular mass but did not find any change, an indication that diffuse flux was resolved out due to the lack of short spacings during the ALMA observation.

3.2.3 Disk Modelling

Kinematically, all three tracers (stars, molecular gas, and ionized gas) show strong signs of circular rotation (see Figures 3 and 4. Therefore, we can fit the velocity fields with a model of rotation to better constrain the inclination and dynamical mass driving the rotation.

We first fit both the CO (2–1) and stellar velocity fields using an updated Python version of dysmal (Davies et al., 2011), now called dysmalpy, an MPE developed software package for modelling the dynamics of galaxies. dysmalpy works by first simulating the mass distribution of a galaxy and populating a three dimensional cube with the corresponding physical velocities. After rotating and inclining the galaxy to match the preferred line of sight, the physical 3D cube is integrated along the line of sight to produce an observed “IFU” cube which is finally convolved with both a PSF and line spread function. In this way, we account for all instrumental effects including beam smearing. These simulated cubes can then be used to compare and fit against observed data using either standard non-linear least squares or a Monte Carlo Markov Chain (MCMC) algorithm with emcee (Foreman-Mackey et al., 2013).

We assumed a thin exponential disk mass distribution parameterised by the total dynamical mass ($M_{\rm dyn}$) and the effective radius ($r_{\rm eff}$) within which the enclosed mass is one half of $M_{\rm dyn}$. The inclination of the disk ($i_{\rm disk}$) and the position angle (P.A.disk) were also left as free parameters in the fitting. For the CO (2–1) model we fixed the centre by eye because we also masked out the inner 2.5″due to the presence of the non-circular motion described below. For the stellar velocity field we allowed the centre to vary between +/- 5 pixels in both the x and y direction. We ran MCMC with 100 walkers for 100 burn-in steps and 400 sampling steps to determine the best fit parameters and uncertainties from the sampled posterior distribution using the median and the 16th and 84th percentiles.

Figure 9 shows the results of our fitting. For the molecular gas, we find a best fitting $\log M_{\rm dyn}=10.372^{+0.002}_{-0.001}$ M*☉, $r_{\rm eff}=537^{+2}_{-2}$ pc, $i_{\rm disk}=43.3\degr^{+0.1}_{-0.1}$, and P.A.${}_{\rm disk}=14.01\degr^{+0.04}_{-0.03}$. We note that the very small errors on the model parameters only incorporate the statistical uncertainties and thus are more reflective of the high S/N and high spectral resolution of the ALMA data. These results are also consistent with those found using the same data in Ramakrishnan et al. (2019). For the stellar velocity field, we find a best fitting $\log M_{\rm dyn}=10.14^{+0.01}_{-0.01}$ M☉*, $r_{\rm eff}=640^{+10}_{-10}$ pc, $i_{\rm disk}=42.1\degr^{+0.4}_{-0.5}$, and P.A.${}_{\rm disk}=10.9\degr^{+0.1}_{-0.1}$.

The residuals in the stellar velocity field all have an absolute value less than 40 km s*-1* and the discrepancies are likely due to deviations from pure circular rotation since we expect stars in the circumnuclear regions to follow more elliptical $x1$ and $x2$ orbits. This can also be seen in measured velocity field where the line of nodes are not exactly perpendicular to the zero velocity line. These deviations are likely the reason for the discrepancy in the best fit model parameters between the molecular gas and stars. However, the good agreement between the models and the data as well as between the cold molecular gas and the stars indicates that our interpretation of a rotating circumnuclear ring is correct and that the young stars have been born out of the molecular gas and still exhibit the molecular gas kinematics. Finally we do note that the zero velocity line of both the stars and molecular gas are S-shaped which could be an indication of inflow of both stars and gas.

3.2.4 Molecular gas outflow and inflow?

The velocity residuals after removing circular rotation could be due to non-circular inflow and outflow. To investigate this, we placed circular apertures onto the residual velocity field spaced in increments of 0.56″ along the minor axis of the best fit disk model. Placement along the minor axis ensures the velocity residuals are primarily due to radial motion rather than differences in tangential motion compared to circular rotation.

Radial inflow will appear as redshifted residuals on the near-side of the disk and blueshifted residuals on the far-side. We know based on the large scale dust maps and H$\alpha$ and [OIII] flux maps as well as the excitation of the ionized gas that the near-side is West of the kinematic major axis. Radial outflow motion will have the opposite sense as radial inflow. We calculate the total mass in the aperture by summing the integrated flux of all pixels in the aperture and using Equation 2 and $\alpha_{\rm CO}$ as before. The mass rate through each aperture is then $Mv/\delta r$ where $\delta r$ is the diameter of the apertures.

Figure 10 plots both the radial velocity estimates (top) and the mass rate (bottom). Positive values of each correspond to outflow while negative values correspond to inflow. It is from this plot that the residuals suggest an outflow velocity of $\sim$80 km s*-1* and mass outflow rate of $\sim 0.5$ $M_{\sun}$ yr*-1* within 1 kpc. We compare these values later to those found for the ionized gas.

Figure 10 also suggests radial inflow outside of 1kpc with velocities reaching 50 km s*-1*. The mass inflow rate into the circumnuclear ring however is very small, only $\sim 0.1$ $M_{\sun}$ yr*-1*. The strongest inflow residuals though do occur in the NW and SE edges of the ring, exactly where the dust lanes meet the ring and where we would expect to find signatures of inflowing gas.

3.3 The Inner Warped Disk

Inside the ring, both the stars and gas dramatically change their distribution and kinematics. The bottom row of Figure 5 shows a zoom in of the inner 6.5″for the F160W image and $V-H$ dust map. We see inside from the contours of the continuum image and elongation along a P.A.$\sim$ 85°shown as a white line. This was first observed by Shaw et al. (1993) who suggested the presence of a stellar nuclear bar possibly supported by the $x_{2}$ orbits of the main large scale bar. Wilson et al. (1993) however also suggested it could just be scattered nuclear light off dust from the AGN.

The dust map shows heavy extinction along an arc extending from the NE to SW directly across the nucleus. This is also the location of two strong concentrations of cold molecular gas that connect two sides of the ring. These two clumps straddle the very centre such that there is a lack of CO (2–1) emission at the location of the AGN and correspond to the double peaks seen in lower resolution CO (1–0) maps Petitpas & Wilson (2002). Perpendicular to the two clumps and the stream are also two distinct holes of CO (2–1) emission. These holes are aligned in the direction of the NLR and ionized gas outflow, and each have a radius of roughly 165 pc and cover a surface area of 8.6$\times 10^{4}$ pc2.

Using a 2″aperture and Equation 2, we measure a molecular gas mass of $\sim 4\times 10^{7}$ M*☉*, or $\sim 25$% of the total molecular gas mass detected. Thus, while a large fraction of the molecular gas is contained within the rotating ring, a substantial amount seems to be experiencing irregular motion around the AGN.

The CO (2–1) line profiles near the centre contain multiple kinematic components which contribute to the increase of the measured velocity dispersion in the single Gaussian fits and are evidence for irregular motion. These kinematic components could be related to either the nuclear stellar bar or the ionized gas outflow so we performed a second fit of only the central 4″x 4″section of the ALMA cube. To fit each spaxel, we used the method described in Fischer et al. (2017) that employs the Bayesian MultiNest algorithm (Feroz et al., 2009; Buchner et al., 2014) to fit emission line profiles and calculate the evidence for a particular fit. In this way, by comparing the evidence from a single component fit to a double component fit, we can determine whether the second component is statistically necessary. This can be extended up to any number of components, and in this work we allow for a maximum of three components based on visual inspection of the ALMA cube.

A complication after fitting for multiple kinematic components is associating each velocity component within each spaxel with the ones in the other spaxels to form a coherent picture of the gas kinematics. To accomplish this, we first determined the likely independent components that existed in the region by fitting the integrated CO (2–1) spectrum within the central 2″. We found five distinct kinematic components that have velocities of roughly -200, -100, 0, 100, and 200 km s*-1* which are shown together with the observed spectrum in Figure 11. Then, for each of the components for a single spaxel fit, we associated it with one of the five components based on which one it was closest to in velocity. In this way, we constructed flux, velocity, and dispersion maps for each of the five components. Finally, based on inspection of each of the maps it was clear that the first and fifth components as well as the second and fourth components are the opposite sides of the same physical component so we plot them together. Figure 12 shows the results of our fits and reconstruction of each molecular gas component where Comp1 refers to the original first and fifth components ($\pm$200 km s*-1*), Comp2 refers to the original second and fourth components ($\pm 100$ km s*-1*), and Comp3 is the original third component (0 km s*-1*).

Our analysis therefore suggests the presence of three kinematically distinct cold molecular gas components in the nuclear region of NGC 5728. The first, shown in the top row of Figure 12, is characterised by LOS velocities near 200 km s*-1* and a P.A.$\sim$0°. The second, shown in the middle row, is characterised by LOS velocities around 100 km s*-1* and a P.A.$\sim$45°, while the third component has LOS velocities near systemic and a kinematic P.A. that is unclear since the velocity gradient is weak but a flux P.A. near 90°.

We can also determine the mass of molecular gas that is contained within each kinematic component using again Equation 2 and the conversions given above. We find a total mass within 2″ for Comp1, Comp2, and Comp3 of $1.6\times 10^{7}$, $1.6\times 10^{7}$, and $7.7\times 10^{6}$ M*☉* respectively. Thus about 40%, 40%, and 20% of the total nuclear molecular gas is contained in Comp1, Comp2, and Comp3.

On all of the flux maps, we also plot the contours of H2 (1–0) S(1)emission from the SINFONI cube. On the velocity maps, we plot the kinematic major axis for the circumnuclear rotating disk based on our fit in Section 3.2.3, the ionized gas bicone region based on the MUSE [OIII] map, and the major axis of the nuclear stellar bar based on Figure 5.

Immediately it is clear that Comp1 likely corresponds to gas still in circular rotation in the inner part of the disk. However, it seems that as the molecular gas is rotating towards the bicone, the gas is being accelerated giving rise to the enhanced velocities seen here and in the residuals of the disk modelling (see Fig. 9). However, as soon as the gas rotates into the bicone, all CO (2–1) emission disappears given none of the components we identified show any correspondence with the bicone. We discuss possible reasons for the disappearance of CO (2–1) emission in Section 3.7.

Comp2, however, seems to be aligned perpendicular to the bicone at least within the central 1″. Based on our later modelling of the ionized gas, we find an inclination for the bicone of 49°compared to an inclination for the disk of 43°. Since the bicone is directed almost straight into the disk, the molecular gas must warp by 90°in the inner nuclear region. Based on the P.A. shift for Comp2 compared to the larger scale disk and Comp1, we associate Comp2 with a warped inner rotating molecular disk that is contributing to shaping the bicone and likely providing the material necessary to fuel the AGN. Comp2 is further directly associated with the strong extinction seen in Figure 5 and is therefore also providing the obscuration necessary to hide the BLR. It also has a relatively high velocity dispersion of $\sim 80$ km s*-1* also indicating a thick disk which are commonly found in the central few 100 pc of AGN (Hicks et al., 2009; Riffel & Storchi-Bergmann, 2011a; Alonso-Herrero et al., 2018).

Comp3 perhaps is related to gas streaming along the nuclear bar. The velocities near 0 km s*-1* therefore are due to the orientation of the bar in the plane of the sky. Petitpas & Wilson (2002) were searching exactly for this molecular gas aligned with the nuclear stellar bar but were unable to detect it likely due to the low spatial resolution of their data. The thought was that nuclear gas in the bar would be an indication that the nuclear stellar bar would be driving gas further inward towards the AGN. However, we see here that while it does seem that there is molecular gas associated with the nuclear bar, the primary component that is driving gas to the centre is instead Comp2 which itself is arranged in a bar-like structure. This nuclear molecular gas bar then is significantly leading the stellar bar and in fact could be related to the $x_{2}$ orbits of the nuclear stellar bar which are perpendicular to the major axis of a bar. Hydrodynamical simulations do indicate that gaseous bars should lead stellar bars when both exist (Friedli & Martinet, 1993; Shaw et al., 1993) however note these simulations primarily consist of a single stellar and gaseous bar whereas for NGC 5728 we are observing a secondary nuclear and gaseous bar. The three molecular gas components make for a seemingly chaotic environment around the AGN and is likely not a dynamically stable setup without collisions between gas clouds. It’s possible though that this contributes to the variability of AGN where short-lived, unstable molecular gas environments quickly feed the AGN before breaking itself apart with the help of AGN feedback.

The morphology of H2 (1–0) S(1) emission seems to be a combination of all three components. The distinct S-shaped structure looks to follow the edges of Comp1 while the slight extension in the middle from NE to SW is either related to Comp2 or Comp3. The biggest difference between H2 (1–0) S(1) and CO (2–1) is that the peak of H2 (1–0) S(1)is directly located on the position of the AGN while the peaks of CO (2–1) straddle the centre. This combined with extended morphology following the edges of CO (2–1) emission suggests that the H2 (1–0) S(1) emission is produced by either X-ray heating from the AGN or shocks as the AGN driven outflow expands into the disk and is consistent with the conclusions of a study of NIR line emission in 62 AGN from long slit spectra (Rodríguez-Ardila et al., 2004; Rodríguez-Ardila et al., 2005; Riffel et al., 2013b). Indeed Durré & Mould (2018) show that this is the case through the use of H2 excitation diagrams.

3.4 AGN Driven Outflow

3.4.1 NLR Morphology

From the MUSE [OIII] maps (Figure 4) we find the NLR in a biconical shape that seems to be “piercing” through the centre of the ring and has an apex exactly at the the location of the AGN. Particularly noticeable is the strong dip in [OIII] emission towards the NW that indicates obscuration by the nuclear star-forming disk. Combined with the fact that the SE emission is systematically brighter than the NW emission, the NLR is inclined such that the bicone is in front of the nuclear disk towards the SE and behind the disk towards the NW until it breaks through at larger radii. This interpretation is consistent with previous studies of NGC 5728 using HST narrow band imaging (Wilson et al., 1993) and lower spatial resolution IFU observations (Davies et al., 2016).

We measure a biconical NLR extent of 1.7 kpc to the SE and 2.1 kpc to the NW. The central PA is roughly -60°and the full opening angle is $\sim$55°. The NLR emission, while extended out to kpc scales, is largely concentrated in two clumps on opposite side of the AGN. The brightest clump in the SE is only 130 pc away from the AGN while the NW one is about 230 pc. Separating the two clumps is another gap in emission that spatially corresponds to the Comp2 component of the cold molecular gas shown in Figure 12 and suggests that it is the component that is obscuring the AGN and producing the biconical shape of the NLR. The reason that the gap in [OIII] emission does not fully align with the Comp2 molecular gas component is due to the inclination of the bicone. The SE side of the bicone is in front of Comp2 while the NW side is behind it such that the gap in [OIII] emission appears to the NW of the location of Comp2. What is abundantly clear though is that the cavities of cold molecular gas are filled with warm ionized gas.

On smaller scales from the SINFONI data (Figure 19), we find that the ionized emission is largely confined to the NLR showing again a biconical structure and extending from the SE to the NW, similar to the [OIII] emission from MUSE. The brightest emission for all three emission lines ([FeII], [SiVI], and Br$\gamma$) is concentrated towards the SE, and with the higher spatial resolution of SINFONI, we locate it only 20 pc away from the AGN. Interestingly, [FeII] shows the same “hook” structure to the NW as the H2 (1–0) S(1) emission, except at locations inward from H2 (1–0) S(1). Both [SiVI] and Br$\gamma$ instead are distributed in distinct clouds and have an overall shorter extent compared to [FeII] and H2 (1–0) S(1). Figure 13 highlights the spatial differences between the different phases of the gas with warm molecular gas (red) encompassing the partially ionized gas (green; [FeII]) which then encompasses the fully ionized gas (blue; [SiVI]). This observed “stratification” of the ionization structure of the gas around the nucleus strongly suggests that the AGN is the source of the ionization. This is also further evidence that the so-called “Coronal Line Region” is simply the inner part of the NLR where the ionizing radiation field is stronger and can produce lines with higher ionization potential.

3.4.2 Ionized Gas Kinematics

In the nucleus, and spatially corresponding to the brightest clumps of ionized gas emission, we find elevated velocities reaching +/- 400 km s*-1*. The kinematic P.A. is along the same direction as the NLR and we see larger dispersions, although some of them are due to beam smearing and the presence of multiple kinematic components as in the cold molecular gas. All of this evidence points towards the presence of an AGN driven outflow, and the high spatial resolution SINFONI observations only further increase support for this interpretation. All of the SINFONI detected ionized emission lines show high velocities along the major axis of the NLR and dispersions reaching over 300 km s*-1*(see Figure 19). In particular, the [SiVI] line displays the largest velocities reaching 500 km s*-1*. Both Br$\gamma$ and [FeII] show in general the same velocity structure albeit with slightly lower absolute velocities. This is in contrast to the H2 (1–0) S(1)kinematics which show much lower velocities and dispersions and suggests there is a lack of molecular gas in the outflow due to the AGN inflating a bubble in the disk such as the one seen in NGC 1068 (May & Steiner, 2017).

3.4.3 Outflow Modeling

Building on the outflow interpretation of the ionized gas kinematics, we used DYSMALPY to model and fit for the outflow parameters. As we did for modelling the rotating disk, we used the 2D velocity map as our input and specifically decided on the [SiVI] velocity map since it shows the strongest outflow signatures. Our biconical outflow model is similar to those that have been extensively used in the literature for modelling the kinematics of gas around nearby AGN based on both long slit data (e.g. Crenshaw & Kraemer, 2000; Crenshaw et al., 2000; Crenshaw et al., 2010; Das et al., 2005, 2006; Fischer et al., 2013, 2014) as well as IFU data (e.g. Müller-Sánchez et al., 2011; Bae et al., 2017). The model consists of two axisymmetric cones that share an apex at the location of the AGN. The cones have two opening angles, an inner ($\theta_{\rm inner}$) and outer ($\theta_{\rm outer}$) opening angle that define the walls of the cone and where the emission originates. The bicone also has an inclination ($i_{\rm bicone}$ where 0° indicates the cones are oriented along the line of sight) and position angle (P.A.bicone) on the sky with respect to the line of sight.

As in Bae & Woo (2016), we assume the flux is exponentially decreasing with radius according to $F(r)=Ae^{-\tau r/r_{end}}$ where $r_{\rm end}$ is the radial edge of the bicone and $\tau$ controls the speed at which the flux falls off. Given we are only fitting the velocity map, we fixed $\tau$ to be 5 which is the value found for many nearby AGN. Finally, the model allows for different velocity profiles that change as a function of radius. We chose a linearly increasing profile that reaches a maximum ($v_{\rm max}$) at a specific radius ($r_{\rm turn}$) then decreases linearly until $r_{\rm end}$. This assumes a radially symmetric profile such that $r_{\rm end}=2r_{\rm turn}$. As with modelling the rotating disk, we determined the best fit parameters using MCMC with 800 walkers, 100 burn steps, and 100 sampling steps.

Figure 14 shows the results of our outflow fitting. We find best fit outflow parameters of $\theta_{\rm inner}=18\pm 1\degr$, \theta_{\rm outer}=23\pm 1\degr$$i_{\rm bicone}=49\pm 2\degr, P.A.${}_{\rm bicone}=43\pm 2\degr$, $r_{\rm end}=540\pm 20$ pc, and $v_{\rm max}=737\pm 25$ km s*-1*based on the median and standard deviations of the resulting posterior distributions.

Comparing the spatial orientations of the disk and outflow, we see that our initial interpretation that the NLR is in front of the disk towards the SE and behind the disk towards the NW is correct. The disk is inclined at an angle of 43° with respect to the LOS, while the outflow is inclined by 49°. However, an outer half opening angle of 23° places the front side of the NLR at only an angle of 26°, well in front of the disk. The similarity of the inclinations for both the disk and the outflow further shows that the NLR is completely intersecting the disk. The turnover radius for the outflow is only 270 pc which explains why we only observe enhanced velocities very near to the centre in the MUSE maps. Thus, while the NLR extends out to 2 kpc only the very central regions are outflowing and the rest is simply AGN excited gas within the disk.

3.5 Outflow Energetics

In this section, we explore the ionized gas energetics within the outflow. Because of the short extent of the outflow, we primarily use the diagnostics available from the SINFONI cube to spatially resolve the outflow properties. In particular, we want to measure the mass outflow rate, mechanical energy rate, and momentum rate as a function of radius and compare to the energy and momentum injection from the central AGN.

Similar to our analysis of the electron density radial profile in Appendix B, we binned the SINFONI cube using annular apertures along the SE redshifted half of the NLR given extinction is not as much of a factor. For the SINFONI data, we chose an annular width of 0.075″ which is half the PSF size and extended out to 1.65″. Within each aperture, we fit the [SiVI] and Br$\gamma$ line profile with either 1 or 2 Gaussian components, choosing between them based on $\log$(evidence) as we did for the very central regions of the ALMA data. If the profile needed two components, we calculated a flux-weighted average velocity since we assume that all of the emission in the inner regions is related to the outflow.

The upper panel of Figure 15 plots the resulting velocity profile for both [SiVI] and Br$\gamma$ with all radii and velocities deprojected based on the inclination from the best fit 2D model from Section 3.4.3. Both ionized gas lines show essentially the same profile, implying that both lines are tracing the same physical components of gas. The narrow and broad components are also both increasing with radius up to $\sim$250 pc, then decreasing toward larger radii just as we modelled previously. The general consistency of both the narrow and broad component suggests that these are also tracing the same gas and merely reflect the underlying velocity distribution of the gas. We also show the few inner bins of the [OIII] velocity profile from MUSE which traces roughly the same shape except a bit more flattened due to the larger beam size of the MUSE observation. We further show the velocity profiles from the stars and the best fit rotating disk model used to fit the CO (2–1) velocity map. These both show a nearly flat profile around 0 km s*-1* since the bicone axis is aligned nearly along the minor axis of the disk.

Because Br$\gamma$ and [SiVI] originate in the same gas we can measure the ionized gas mass using the total flux of the Br$\gamma$ line assuming Case B recombination gas conditions. Further since Br$\gamma$ occurs in the NIR, dust extinction is negligible. We convert from Br$\gamma$ flux to gas mass within each aperture using the formula given in Storchi-Bergmann et al. (2009):

[TABLE]

where $D$ is the luminosity distance in Mpc and $n_{e}$ is the electron density. $M_{\scriptscriptstyle\rm HII}$ then is the total ionized gas mass in solar units. For $n_{e}$, we assume a constant electron density within the SINFONI FOV and use the value obtained from our analysis of the auroral and trans-auroral lines (1550 cm*-3*; see Appendix B). We can then calculate the mass outflow rate, mechanical luminosity, and outflow momentum rate using the following standard equations:

[TABLE]

[TABLE]

[TABLE]

where $\delta r$ is the width of the annular apertures. The bottom three panels of Figure 15 show our calculations of $\dot{M}_{\rm{out}}$, $\dot{E}_{\rm{out}}$, and $\dot{p}_{\rm{out}}$ as a function of distance from the AGN. $\dot{M}_{\rm{out}}$, $\dot{E}_{\rm{out}}$, and $\dot{p}_{\rm{out}}$ all show peak values around 200 pc and peak values of 0.08 M*☉* yr*-1*, $4\times 10^{39}$ erg s*-1*, and $2\times 10^{32}$ Dyne. We note that our outflow measurements disagree from those presented in Durré & Mould (2019) who reported $\dot{M}_{\rm{out}}=38$ M*☉* yr*-1*, almost 500 times greater than ours. We attribute this difference to two factors: 1) they assumed an electron density of 100 cm*-3* compared to our measured value of 1000 cm*-3* and 2) they summed up the outflow rate determined within each single spaxel over the entire SINFONI FOV. For mass outflow rates, this is not correct as the rate needs to be calculated within a common radius.

The peak mechanical energy in the outflow is well below the AGN luminosity, $L_{\rm Bol}=1.4\times 10^{44}$ erg s*-1*, only reaching 0.003% of $L_{\rm Bol}$. Similarly, the ratio of the momentum rate to the radiation momentum rate from the AGN ($L_{\rm Bol}$/$c$) is only 0.043. These properties all suggest that the outflow currently is not energy driven but rather momentum driven given their extremely low values (e.g. Faucher-Giguère & Quataert, 2012; Zubovas & King, 2012).

While the values seem quite low, we can compare them to the properties of outflows found across literature. Fiore et al. (2017) recently compiled a comprehensive set of outflow energetics from various studies over the years spanning a large range of AGN luminosity. They used this data to establish clear AGN wind scaling relations, in particular between the mass outflow rate and outflow power with the AGN bolometric luminosity in all phases of the gas in the outflow.

In Figure 16, we compare our values of $\dot{M}_{\rm{out}}$ and $\dot{E}_{\rm{out}}$ for NGC 5728 with the ionized wind data from Fiore et al. (2017). NGC 5728 lies very close to the correlations and the low values of $\dot{M}_{\rm{out}}$ and $\dot{E}_{\rm{out}}$ seem to be simply due to the relatively modest AGN luminosity. With only one point, we cannot definitively determine if the correlations presented in Fiore et al. (2017) actually extend to lower luminosities. Therefore, we compiled more values for $\dot{M}_{\rm{out}}$ and $\dot{E}_{\rm{out}}$ with a focus on AGN with $L_{\rm Bol}\la 10^{46}$ erg s*-1*. We found data for 15 AGN: NGC 1068, NGC 2992, NGC 3783, NGC 6814, NGC 7469 (Müller-Sánchez et al., 2011), NGC 4151 (Crenshaw et al., 2015), Mrk 573 (Revalski et al., 2018a), Mrk 34 (Revalski et al., 2018b), NGC 7582 (Riffel et al., 2009a), Mrk 1066 (Riffel & Storchi-Bergmann, 2011a), Mrk 1157 (Riffel & Storchi-Bergmann, 2011b), Mrk 79 (Riffel et al., 2013a), NGC 5929 (Riffel et al., 2015), NGC 2110 (Schnorr-Müller et al., 2014), and NGC 3081 (Schnorr-Müller et al., 2016). Finally, we also include the sample of 234 type 2 AGN studied in BN19 using single integrated SDSS spectra and dust SED fitting.

The inclusion of both single source literature compilation and the BN19 sample shows that while NGC 5728 lies close to the Fiore et al. (2017) correlations, there is likely substantial scatter at lower luminosities. That said, combining all of these samples together can be dangerous given the widely varying methods and data quality used to measure these properties. In particular, while Fiore et al. (2017) renormalized all outflow properties to a constant electron density, BN19 showed that object by object determined densities can change the mass outflow rate considerably and thus weaken the mass outflow rate correlation. What is needed to confirm these correlations is a large, spatially resolved sample of AGN where all of these properties can be properly and systematically measured.

3.6 AGN Feedback in NGC 5728

With such a low mass outflow rate and outflow power, it is unlikely that the AGN driven wind is effectively suppressing star formation in NGC 5728 and affecting its evolution. With a total molecular gas mass in the circumnuclear disk of 108.1 M*☉, it would take 217 Myr to completely clear the disk of its gas even including the potential molecular gas outflow. Using Herschel photometry and SED fitting, Shimizu et al. (2017) estimated a global SFR of 2.1 M☉* yr*-1*. Combining this with $\dot{M}_{\rm{out}}$, the total gas depletion time still 47 Myr. However, this is under the assumption that all of the outflowing gas can escape the gravitational potential of the galaxy. Following Rupke et al. (2002), we can estimate the escape velocity under the assumption of a truncated isothermal sphere, $v_{\rm esc}(r)=\sqrt{(}2)v_{\rm circ}[1+\ln(r_{\rm max}/r)]^{1/2}$. We assume an $r_{\rm max}=100$ kpc and use $v_{\rm circ}=320$ km s*-1* from Rubin (1980). At 250 pc, the radius of peak outflow velocity, $v_{\rm esc}=1200$ km s*-1*, far above the outflow velocities of either the ionized or molecular gas. Outflows in moderate luminosity AGN do not seem to have the energy to escape their host galaxies as this has been observed in many other systems (e.g. Davies et al., 2014; Fischer et al., 2017; Herrera-Camus et al., 2018).

We can further compare the mass outflow rate with the mass flow rate through the bicone. Using the best fit parameters for the outflow, the total surface area covered by the bicone is 0.19 kpc2. From our disk modelling we find an average rotational velocity of the molecular gas of 273 km s*-1* through the bicone. Combined with the average surface mass density of $55$ M*☉* pc*-2*, we calculate that molecular gas flows through the bicone at a rate of $\sim$ 9.5 M*☉* yr*-1*, about 16x higher than the mass outflow rate. Essentially, the gas that is entering the bicone is not being ejected at rates large enough that would substantially disturb the structure and dynamics of the disk. This is especially true when we take into account the fact that the molecular gas is only moving radially at speeds of 40 km s*-1* which means that by the time it exits the bicone it has only moved about 55 pc.

Finally, we can compare the mass outflow rate to the current mass accretion rate onto the SMBH using the following equation:

[TABLE]

where $\eta$ is the efficiency of accretion to convert rest mass energy into radiation and $c$ is the speed of light. Assuming $\eta\approx 0.1$ for a standard geometrically thin, optically thick accretion disk, we find a current $\dot{M}_{\rm{acc}}\approx 0.025$ M*☉* yr*-1*. Therefore, even the ionized gas mass outflow rate is three times the current mass accretion rate onto the SMBH implying that the majority of gas near the AGN does not end up fuelling it.

3.7 Missing molecular gas?

The most striking evidence of AGN feedback is the presence of what appears to be large “holes” of cold molecular gas as probed by CO (2–1) emission that are co-located exactly along the outflow and NLR. Therefore, it’s reasonable to interpret that the AGN driven outflow has completely evacuated these regions of the disk of molecular gas. However, these holes do not exist in maps of the ionized gas and in fact it is these regions where emission from ionized gas is strongest. If an AGN driven outflow had indeed evacuated these regions, we should expect it to drive out all phases of gas.

A plausible explanation for these holes is instead that CO (2–1) is not accurately tracing the presence molecular gas in the NLR. Under standard ISM conditions, CO (2–1) would be a reliable tracer of the molecular gas mass, however in the presence of the strong radiation field of an AGN, it could be that the CO (2–1) line in particular or all CO lines are suppressed. Whereas in a normal star-forming disk, molecular gas is primarily heated by nearby OB stars with FUV photons dominating the process, near an AGN, X-rays will dominate over FUV photons and control the heating and chemical composition of the gas producing a so-called X-ray dissociation region (XDR).

This is being increasingly observed in local AGN. Both Rosario et al. (2019) and Feruglio et al. (2019) reported recently on a lack of CO (2–1) emission near the central AGN that is instead filled with H2 (1–0) S(1) and ionized gas emission. These cavities of CO emission are further cospatial with hard X-ray emission, confirming that X-ray radiation is suppressing the CO emission. While we observe the same effect at the very centre of the galaxy corresponding to the location of the AGN, we do not observe H2 (1–0) S(1) emission filling the larger scale CO (2–1) holes that correspond to the bicone. Instead, H2 (1–0) S(1) is only lining the edges of the bicone where gas is entering the bicone. It’s possible then that the inner NLR is devoid of molecular gas and a continuation of the stratification we observe in the SINFONI line emission continues on larger scales, i.e. that cold molecular gas only exists outside the bicone region with hotter and hotter gas becoming more prevalent towards the inner NLR. To confirm the lack of molecular gas in the NLR, however, observations of higher-J CO lines are needed.

4 Summary & Conclusions

We have presented a comprehensive analysis of multiwavelength datasets for the nearby Seyfert 2 galaxy, NGC 5728. Using primarily HST imaging and MUSE, SINFONI, and ALMA cubes, we measured and analysed the distribution and kinematics of the stars, ionized gas, and molecular gas. Our results and conclusions are summarised as follows:

• •

Prominent dust lanes are observed extending along the major axis of the large stellar bar from the outer disk to the circumnuclear ring. These are likely coincident with the primary shocks produced by the axisymmetric bar potential and driving inflow to the circumnuclear region.

• •

The circumnuclear ring is observed as a ring of low stellar velocity dispersion, bright H$\alpha$ clumps, and CO (2–1) emission indicating the presence of on ongoing star formation and young stars. This star formation is likely induced by the build up of gas at the Inner Lindblad Resonance of the primary bar. The kinematics of the ring are well fit by circular rotation, however we also seen signs of inflow into the ring where the ring and dust lanes meet at a rate of 1 M*☉* yr*-1*.

• •

Inside the ring, we find a three distinct kinematic components of the molecular gas corresponding to gas rotating with the ring, gas inflowing to the AGN, and gas following the nuclear stellar bar. The gas inflowing to the AGN is distributed across the nucleus at the location of heavy extinction and is the source of obscuration for the AGN due to its perpendicular orientation compared to the NLR.

• •

The AGN is driving a weak outflow primarily seen in ionized gas with a mass outflow rate of 0.08 M*☉* yr*-1* that only reaches to radii of 250 pc before decelerating down to stellar rotation velocities. There are hints of molecular outflow in residuals of the disk modelling indicating a molecular mass outflow rate of 1 M*☉* yr*-1* but only 40 km s*-1*velocity. We determine that the outflow is unlikely to be largely disturbing the structure of the circumnuclear disk.

• •

We observe cavities of CO emission co-spatial with the AGN outflow that could either be explained by a deficiency of molecular gas in the NLR or a suppression of CO emission by the hard X-ray radiation of the AGN.

Overall, for NGC 5728, we find signatures of both feeding and feedback of the AGN. However, feeding of both the AGN and central region of the galaxy is currently the dominant process. NGC 5728 would need a large increase to its accretion rate to boost the energetics of its outflow and severely disrupt its evolution.

Acknowledgements

We thank the anonymous referee for helpful comments and suggestions which improved the paper. We thank Cheng-Yu Kuo from the Megamaser Cosmology Project for sharing the coordinates of the AGN in NGC 5728. TTS thanks Travis Fischer for sharing and discussing his line decomposition software package as well as discussing the results of this paper. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2015.1.00086.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan) and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This research has made use of NASA’s Astrophysics Data System Bibliographic Services. This research made use of Astropy,555http://www.astropy.org a community-developed core Python package for Astronomy (Astropy Collaboration et al., 2013, 2018).

Appendix A Full FOV Emission Line Fits

Appendix B Electron Densities

An important diagnostic IFU studies can provide is the electron density ($n_{\rm e}$) of the ionized gas. The electron density is needed in particular for the determination of mass outflow rates for both star formation and AGN driven outflows. Historically, these have been measured from single spectra using the [SII]$\lambda\lambda$6716,6731 line ratio with high uncertainties or assumed based on local studies with usual estimates of $n_{\rm e}<100\,\rm{cm}^{-3}$. Recently, IFU studies have begun to spatially resolve the electron density distribution within individual galaxies at low redshift (e.g. Kakkad et al., 2018; Nascimento et al., 2019). At high redshift, stacks of the nuclear spectra from IFU surveys have produced estimates of the electron density for both HII regions and broad outflows (Förster Schreiber et al., 2019). These studies have all indicated that the assumed electron densities have been underestimated which in turn can greatly overestimate the mass outflow rate.

We use three independent methods to determine the electron density in the circumnuclear region of NGC 5728. The first and most popular method in the literature uses the [SII] doublet ratio as the density diagnostic due to the weak dependence on electron temperature. With a critical density around 1000 cm*-3* the [SII] line ratio provides reliable electron densities within the range of 50–2000 cm*-3*. Above or below these densities however, the [SII] doublet ratio saturates and becomes unresponsive to higher or lower densities.

The second method, introduced by Holt et al. (2011) and used extensively in Rose et al. (2018) and Santoro et al. (2018), incorporates the ratios of four auroral or trans-auroral emission lines: the [OII] doublets $\lambda\lambda$3726,3729 and $\lambda\lambda$7319,7330 and the [SII] doublets $\lambda\lambda$4069,4076 and $\lambda\lambda$6716,6731. This method constructs line ratios from the total flux between the doublets instead of the flux in individual lines and is therefore less susceptible to degeneracies in decomposing the doublets. The large wavelength range further provides an estimate of the gas reddening. As shown by Holt et al. (2011), this diagnostic is also sensitive to a much larger range of electron densities than the single [SII]$\lambda\lambda$6716,6731 line ratio (102–106.5 cm*-3*). However, both [OII] doublets and the [SII]$\lambda\lambda$4069, 4076 doublet are usually weak and to simultaneously measure all four doublets requires either multiple spectra or a spectrum that covers a large wavelength range.

The third method was developed in Baron & Netzer (2019, hereafter BN19) and we provide a brief summary here. Instead of using the weak auroral and trans-auroral lines, they calibrated a method of measuring $n_{e}$ using the strong [OIII], H$\alpha$, H$\beta$, and [NII] lines under the assumption that the gas is photoionized by an AGN. They found that the [OIII]/H$\beta$ and [NII]/H$\alpha$ line ratios can provide a reliable measure of the ionization parameter following a simple equation which can then be converted into an electron density assuming the radial distance to the AGN is known.

MUSE does not extend to short enough wavelengths for the trans-auroral analysis, but as part of the LLAMA program, we obtained X-Shooter IFU spectra for all of the targets. For NGC 5728, we detected all four of the auroral or trans-auroral doublets in the spectrum extracted from a central 1.8″ box aperture. In Figure 21, we show cutouts of the X-Shooter spectrum around the four doublets. Unfortunately the [OII]$\lambda\lambda$3726,3729 doublet is next to the edge of one of the UVB echelle orders which caused the increased noise blue-ward of the doublet.

To determine the necessary kinematic components to model the doublets, we first fit the [OIII]$\lambda\lambda$4959,5007 emission lines to use as a robust template for the presence, location, and width of the existing kinematic components comprising the ionized gas as we recently did for this same galaxy in Shimizu et al. (2018). In fitting the [OIII] doublet, we fixed the velocity and dispersion of the components to be the same between each of the lines and fixed the intensity ratio (5007/4959) to 3, the theoretical value. The aperture used to produce the X-Shooter spectrum encompasses both sides of the outflow so we used a combination of red and blue-shifted Gaussian components in our model. To properly reproduce the [OIII] line profile, we found we needed 6 kinematic components, 3 red-shifted and 3 blue-shifted.

In fitting the auroral and trans-auroral lines, we fixed the velocity and dispersions of all of these components but left the intensities to vary. However, we found while fitting the [SII]$\lambda\lambda$6716,6731 doublet that we needed to include a seventh component with 0 km s*-1* velocity and $168\pm 7$ km s*-1* dispersion. The remaining 3 doublets were then well-fit with these 7 components. Figure 21 shows the best fit model for each of the doublets as well as the combined red-shifted, blue-shifted, and systemic components of the model.

While the fits overall are satisfactory, the low S/N and relatively small separation of three of the doublets creates strong degeneracies between the amplitudes of the individual components. However, as mentioned previously, the electron density measurement from the trans-auroral lines is only dependent on the total flux of each of the doublets. Therefore, for NGC 5728 the combined auroral and trans-auroral lines can provide a robust estimate of the flux-weighted electron density but cannot provide separate electron densities for each of the kinematic components.

Figure 22 plots the location of the [SII] and [OII] doublet ratios along with the expected values for a grid of models with different $E(B-V)$ and $n_{\rm e}$. The models are the same as those presented in Rose et al. (2018) and were produced using the cloudy photoionization code under the assumptions of solar metallicity gas photoionized by a central source with a power law continuum ($\alpha=-1.5$ for $F_{\nu}\propto\nu^{\alpha}$). They varied the electron density of the gas to model intrinsic trans-auroral line ratios and then applied a Calzetti et al. (2000) reddening law for a range of $E(B-V)$.

Using 2D spline interpolation, we estimated the best-fit $\log\,n_{\rm e}$ and $E(B-V)$ for NGC 5728 assuming they vary smoothly as a function of the trans-auroral line ratios. Taking into account only the uncertainties on the line ratios, we estimate a flux-weighted $\log n_{\rm e}=3.19\pm 0.12$ (1550 cm*-3*) and $E(B-V)=0.18\pm 0.05$.

Within the same X-Shooter spectrum, we then fit the H$\alpha$+[NII] complex and H$\beta$ line profile using the same kinematic components to determine $n_{e}$ via the BN19 method. As a comparison to the trans-auroral estimate, we only calculated [OIII]/H$\beta$ and [NII]/H$\alpha$ line ratios for the total line profile. From the line ratios, we measure a $\log U=-2.59\pm 0.11$. The biggest uncertainty here is what radial position to choose to convert $\log U$ into an electron density. We chose 171 pc which corresponds to half the size of the X-Shooter aperture. We find $\log n_{e}=3.20$ which agrees well with the value from the auroral and trans-auroral lines.

We can compare this to using only the [SII]$\lambda\lambda$6716,6731 line ratio. We find an [SII] 6716/6731 ratio of $0.65\pm 0.34$. Using the pyneb666http://www.iac.es/proyecto/PyNeb/ software package and assuming an electron temperature of 104 K, this translates to a $\log n_{\rm e}=3.36\pm 0.95$ (2300 cm*-3*). Clearly the large amount of blending between the [SII] doublet lines causes severe degeneracies and leads to almost a factor of ten uncertainty even though the value agrees quite well with the trans-auroral value.

We also compared the X-Shooter derived $n_{\rm e}$ with ones obtained from the MUSE cube based on both the [SII] 6716/6731 ratio and the BN19 method. We measured the line fluxes using annular slices along a PA of -65° West of North, ie along the bicone axis. The annular slices spanned angles +/- 20° away from the centre, had a width of 0.5″, and ranged from 0 to 11.5″ away from the AGN. We chose to only use the redshifted SE side of the NLR since it is most unaffected by extinction and has higher S/N. Figure 23 plots the annular apertures on top of the [OIII]$\lambda$5007 flux map from our single Gaussian fits to the MUSE cube and shows how the annuli cover the majority of the biconical structure of the NLR.

Spectra were extracted from each of the apertures and the [OIII] doublet was again fit first to determine the necessary number of kinematic components in the model since the spatial location, aperture size, spectral resolution, and S/N is different from the X-Shooter spectrum. The H$\beta$, H$\alpha$+[NII] complex, and [SII] doublet were then fit with the same components.

Figure 24 plots $n_{\rm e}$ as a function of radius as determined by the [SII] 6716/6731 ratio and the BN19 method from the MUSE cube. Also plotted in the centre is our estimate from the X-Shooter spectrum. Inside of 1kpc, there is a strong increase of $n_{e}$ based on BN19, whereas the [SII] based $n_{e}$ slightly decreases toward the centre. This is due to the strong blending of the [SII] 6716/6731 doublet in the centre where the velocity gradient is large and there are multiple kinematic components. Therefore a robust [SII] doublet measurement of the electron density becomes very difficult without high spectral resolution. The BN19 based values are also very well aligned with the auroral and trans-auroral based value which is the most robust estimate. It’s clear then that within 1 kpc, $n_{e}\approx 1000$ cm*-3*, therefore for all later analysis we will use the auroral and trans-auroral value as the electron density of the gas. At larger radii, we note there is better consistency between the BN19 and [SII] doublet methods but still some discrepancies which could be explained by the different emission regions for the various lines. It’s possible that while [OIII] and H$\alpha$ are being emitted throughout a whole ionized cloud, [SII] instead is only emitted at the ionization front.