© ESO 2021

1. Introduction

Active galactic nuclei (AGN) are compact luminous sources at the very centre of many giant galaxies. They are powered by the potential energy loss of material falling into a supermassive black hole (SMBH; Salpeter 1964; Lynden-Bell 1969). AGN come in a multitude of varieties distinguishable by their spectral properties (see e.g. Table 1 in Padovani et al. 2017). AGN are interesting objects by themselves, and also because of their coupling with their host galaxies. An example of this is the fairly tight correlation between the SMBH mass and the velocity dispersion of the stellar spheroid (the so-called ℳ_BH − σ_⋆ relation; Ferrarese & Merritt 2000; Gebhardt et al. 2000). AGN feedback is thought to be one of the mechanisms that limit the growth of massive galaxies (e.g. Harrison 2017, and references therein) and that contribute to the transformation of dark matter cusps into cores (e.g. Peirani et al. 2008).

Active galactic nuclei feedback mechanisms include outflows (for a review on AGN feedback, see Morganti 2017). The first outflows were detected in ionised gas (see the historical discussion in Veilleux et al. 2005), but they are nowadays known to be multi-phase and also carry H I (Morganti et al. 2005) and molecular gas (Feruglio et al. 2010). The kind of feature studied in this paper, ionised outflows, is sometimes complex and might require a multi-component approach to be accurately described (e.g. McElroy et al. 2015; Lena et al. 2015; Mingozzi et al. 2019). Outflows are thought to be part of the self-regulation mechanism for the growth of the SMBH and to contribute to the quenching of star formation in the nuclear regions of the host galaxy (for a recent review, see Veilleux et al. 2020).

The feeding of AGN is a matter of controversy, since it requires the existence of a mechanism for the gas to lose its angular momentum to reach a galaxy centre. Galaxy-galaxy interactions can generate torques to that effect (Negroponte & White 1983). In non-interacting galaxies, and at scales larger than 1 kpc, inwards gas transportation can be efficiently triggered by energy dissipation at shocks and gravitational torques associated with bars (Schwarz 1984; Athanassoula 1992) and spiral arms (Lubow et al. 1986; Kim & Kim 2014). Strong observational evidence of bars funnelling material towards galactic centres is provided by the detection of enhancements in the star formation, gas concentration, and central mass concentration in barred galaxies (e.g. Heckman 1980; Hummel 1981; Hawarden et al. 1986; Devereux 1987; Sakamoto et al. 1999; Sheth et al. 2005; Díaz-García et al. 2016; Lin et al. 2017; Díaz-García et al. 2020). Large-scale bars and spirals bring the gas to the inner Lindblad resonance (ILR; Schwarz 1984) region, which is located at roughly one kiloparsec from the centre and sometimes traced by spectacular star-forming nuclear rings (Knapen et al. 1995; Comerón et al. 2010). Once near the ILR, it is unclear how the gas loses its remaining angular momentum to move further in, but it has been proposed that this can be achieved by shocks and gravitational torques in a ‘bar-within-bar’ scenario (Shlosman et al. 1989; Hunt et al. 2008; Querejeta et al. 2016) or in a nuclear spiral scenario (Combes et al. 2014; Kim & Elmegreen 2017). The same fuel that feeds the central engine can also ignite intense circumnuclear star formation episodes, or ‘nuclear starbursts’. Galaxies hosting both a Seyfert AGN and a nuclear starburst are referred to as ‘composite’ (Telesco 1988).

The study of the innermost parts of galaxies is crucial to understanding how AGN are fed (inflows; Storchi-Bergmann & Schnorr-Müller 2019) and how they affect their surroundings (through, e.g. outflows Ramos Almeida & Ricci 2017; Hönig 2019). Spectroscopic data are necessary to study both the kinematics and the physical conditions of the circumnuclear medium. Because of the relevant scales (a few hundred parsecs or less) sub-arcsecond angular resolution is required even for the closest galaxies. Hence, the advent of the laser guide star, GALACSI laser adaptive optics (AO) module (Stuik et al. 2006), in the Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) integral field spectrograph at the VLT, provides a new tool to improve our understanding of galaxy-AGN interplay. In narrow field mode (MUSE-NFM), MUSE + GALACSI combine the wide MUSE wavelength range (4750–9350 Å with a gap at 5780–6050 Å to prevent contamination from the laser guide stars) and an extraordinary angular resolution that is nominally below $0\stackrel{″}{.}1$ over the whole wavelength range. As a consequence, MUSE is able to obtain integral field data at an angular resolution comparable to that of the Hubble Space Telescope (HST) over a field of view of about $7\stackrel{″}{.}5\times 7\stackrel{″}{.}5$. Such angular resolutions were achievable in the near-infrared (see e.g. the works by Davies et al. 2009; Riffel et al. 2009, 2010, done with VLT SINFONI and Gemini NIFS data, respectively), but MUSE has expanded these capabilities to optical wavelengths, where the lines necessary to build, for example, Baldwin, Phillips, and Terlevich (BPT) diagnostics are found.

In Knapen et al. (2019), we published the first ever MUSE-NFM AO-assisted study of the circumnuclear medium in an AGN-hosting galaxy, NGC 7130. Now, we build upon our previous work and provide a detailed analysis of the data in order to unveil the complex physics of the circumnuclear medium, including the outflow. In Sect. 2, we summarise the properties of the target galaxy and the findings reported in the literature. In Sect. 3, we describe the data processing, including the reduction and the spectral analysis. In Sect. 4, we describe our results, which are then discussed in Sect. 5, where we also present a simple model to explain the observations. We summarise our findings in Sect. 6.

Throughout this paper, we assume the cosmology derived from the five-year WMAP mission combined with Type Ia supernovae and baryonic acoustic oscillation data (Hinshaw et al. 2009), that is a Hubble-Lemaître constant of H₀ = 70.5 km s⁻¹ Mpc⁻¹, a matter density parameter of Ω_m, 0 = 0.27, and a cosmological constant density parameter of Ω_Λ, 0 = 0.73.

2. Previous studies of NGC 7130

The galaxy NGC 7130, also known as IC 5135, is a southern galaxy found at right ascension $\mathrm{RA}={21}^{\mathrm{h}}{48}^{\mathrm{m}}19\stackrel{\text{s}}{.}520$ and declination $\delta =-34°{57}^{\prime }04\stackrel{″}{.}48$ (Epoch J2000.0) with a redshift z = 0.016151, according to the NED¹. It is a peculiar Sa galaxy (de Vaucouleurs et al. 1991) where infrared observations reveal a bar (Mulchaey et al. 1997). An inspection of the HST images presented in Malkan et al. (1998) and Elias-Rosa et al. (2018) also shows the bar in optical, albeit partially obscured by conspicuous dust lanes. The bar is surrounded by a star-forming inner pseudo-ring (Dopita et al. 2002; Muñoz Marín et al. 2007). The proper, luminosity, and angular-diameter distances are D_p = 64.8 Mpc, D_L = 65.8 Mpc, and D_A = 63.9 Mpc, respectively (based on the velocity with respect to the cosmic microwave background provided by the NED, 4586 ± 23 km s⁻¹). At that distance, one arcsecond corresponds to 310 pc.

The infrared luminosity of NGC 7130 is log (L_IR/L_⊙) = 11.35 (Sanders et al. 2003, who assumed a value of H₀ = 75 km s⁻¹ Mpc⁻¹) so it is classified as a luminous infrared galaxy (LIRG). In this kind of galaxy, the intense infrared emission is usually due to an intense star formation episode often linked to interactions between spiral galaxies (Sanders & Mirabel 1996). NGC 7130 forms a pair with IC 5131 (Sandage & Bedke 1994), which is located at a distance of 12′, or 220 kpc in projection. No obvious signs of interaction between the two galaxies are seen, but the distorted appearance of the outskirts of NGC 7130 (already reported in de Vaucouleurs et al. 1964, 1976) might indicate a past close encounter between them or with a smaller unidentified member of the group. The asymmetric velocity and velocity dispersion maps of the ionised gas in the galaxy (Bellocchi et al. 2012) are further indicators of a likely recent interaction.

The galaxy NGC 7130 is nearly face-on, with an axial ratio q = 0.88 and position angle PA = 160° measured at the K_s-band 20 mag arcsec⁻² isophote (from 2MASS; Skrutskie et al. 2006). Orientation parameters obtained from optical images are very similar (Lauberts & Valentijn 1989).

The AGN of NGC 7130 was originally classified as a Seyfert 2 (Phillips et al. 1983). This was later refined to Seyfert 1.9 (Véron-Cetty & Véron 2010), but see Sect. 5 for further details on the Seyfert type. The composite H II + Seyfert nature of the nucleus of NGC 7130 was independently found by Véron (1981) and Phillips et al. (1983) and confirmed by Thuan (1984) and Shields & Filippenko (1990), but Radovich et al. (1997) claimed that nuclear star formation is not required to explain the spectrum of the inner kpc. The core of NGC 7130 has been found to emit in radio (Norris et al. 1990), and its optical spectrum has two kinematic components, narrow and broad (Busko & Steiner 1988), interpreted to be associated with H II regions and the AGN, respectively (Shields & Filippenko 1990). The broad component is blueshifted and was later hypothesised to be associated with an outflow (González Delgado et al. 1998; Bellocchi et al. 2012; Davies et al. 2014).

High-resolution UV continuum images obtained by the HST reveal a tiny circumnuclear ring 1″ in size (major axis) and a few UV knots along the spiral arms associated with the bar (González Delgado et al. 1998). They also found a UV knot, presumably highly obscured, at the suspected location of the AGN engine. The ring is reminiscent of the ultra-compact nuclear rings (UCNRs) presented in Comerón et al. (2008). The fact that the AGN of NGC 7130 is highly obscured was confirmed by the modelling of the spectral energy distribution (Contini et al. 2002) and by X-ray observations (Levenson et al. 2005). The latter authors found that the AGN emits most of the hardest X-rays in NGC 7130 (> 2 keV), but that two thirds of the total X-ray emission can be ascribed to an extended component associated with star-forming regions.

The centre of NGC 7130 has two dusty spiral arms within the bar. They coincide with molecular gas as traced by the CO(6–5) transition (Zhao et al. 2016). Extended CO emission, maybe partially correlated with the star-forming UCNR, is also found in the central 1″. Zhao et al. (2016) found no traces of a molecular gas outflow using this CO transition (which traces very dense molecular gas). The spectral line energy distribution of CO in the central regions of NGC 7130 requires both star formation and X-Ray emission to be modelled (Pozzi et al. 2017).

In Knapen et al. (2019), we presented the first optical high-angular resolution integral field study of the centre of NGC 7130. We found a tiny bipolar pattern in the velocity maps that we interpreted as a kinematically decoupled core $0\stackrel{″}{.}2$ in radius. We confirmed the presence of an outflow whose line ratios indicate AGN ionisation. There, we assumed the location of the nucleus of NGC 7130 to be at the spaxel where a single-component fit of the ionised gas kinematics yields the largest velocity dispersion. This coincided with the centre of the small bipolar structure and with a bright knot located within the UCNR. Here, we assume the same position for the engine of the AGN.

3. Observations and data processing

3.1. Data obtention and reduction

We obtained MUSE-NFM AO data of the central region of NGC 7130 as part of the science verification programme for this observation mode. Our proposal aimed to obtain 4 × 600 s on-target exposures intertwined with 180 s off-target exposures to model the sky. Ten 600 s exposures were taken on the nights of September 15, 16, and 18 in 2018. Unfortunately, only in the two exposures taken on September 18 did the AO work well enough to bring the seeing below $0\stackrel{″}{.}5$ (full width at half maximum (FWHM) of $0\stackrel{″}{.}16-0\stackrel{″}{.}17$ as measured using the imexamine tool in IRAF; Tody et al. 1986). After publishing our letter (Knapen et al. 2019), we checked the ESO archive for additional data. We were pleasantly surprised to find that eight 300 s MUSE-NFM AO exposures of the same region had been taken during the commissioning of the instrument mode on June 19 and 21 2018. Three of these exposures have an angular resolution comparable to those used in Knapen et al. (2019) and are included in the present work, bringing the total exposure time to 2100 s.

The raw MUSE data were processed using the standard MUSE pipeline (version 2.8.1; Weilbacher et al. 2010, 2014) run under the version 2.9.1 of the EsoReflex environment (Freudling et al. 2013). The five exposures were manually aligned and then combined using the muse_exp_combine recipe. The processed data cube has an angular resolution of $0\stackrel{″}{.}17$, comparable to the UVIS HST images of the same region presented by Elias-Rosa et al. (2018, angular resolution $0\stackrel{″}{.}11$). The combined data cube has two extensions, namely one with the signal and another one with variances (error estimates). In Fig. 1, we show a white-light and continuum-subtracted Hα and [O III] λ5007 images (Sect. 3.3) produced from the final reduced data cube to illustrate the quality of the data.

Fig. 1. Left panel: white-light image of the centre of NGC 7130 obtained from integrating the reduced MUSE data cube along the spectral direction. Middle panel: continuum-subtracted Hα image obtained from the same data cube. The red ellipses indicate the outline of the nuclear rings (see Sect. 4.3). Right panel: continuum-subtracted [O III] λ5007 image. The intensities are scaled logarithmically. In this and all subsequent maps, the plus sign indicates the inferred galaxy centre and the insets show an enlarged version of the centremost region of the galaxy. The axes are in arcseconds. North is up and east is left.

3.2. Use of GIST and pyGandALF

We processed the data cube using a modified version of the Galaxy IFU Spectroscopy Tool² (GIST, version 2.0.0; Bittner et al. 2019) and obtained the emission line properties using the Python implementation of GandALF (Sarzi et al. 2006; Falcón-Barroso et al. 2006), which is called pyGandALF (Bittner et al. 2019), included in GIST. GIST is an all-in-one Python pipeline that comprises many features to easily extract physical information and data cubes. The pipeline is complemented by Mapviewer, an extremely powerful interactive visualisation tool that has been key to understanding our data.

The GIST functions that we were particularly interested in were (1) the Voronoi binning of the data and (2) the extraction of stellar kinematics. How these are implemented in GIST and how these functions were modified for our purposes is explained in detail in Sects. 3.3 and 3.4. In Sects. 3.5–3.7, we explain how we used pyGandALF to obtain the emission line properties.

3.3. Voronoi binning

Voronoi binning was performed with the code written by Cappellari & Copin (2003). The user provides a signal and a noise value for each spaxel that are used to produce bins with a chosen signal-to-noise ratio (S/N). GIST’s original implementation allows for only a single Voronoi binning to be used both for the stellar and the gas emission. This was not convenient for our purposes because our data have a combination of poor S/N stellar emission and high S/N line emission. We thus modified the code to support two binnings.

The stellar emission binning was made using GIST’s original binning procedure, using the rest-frame wavelength range of 5500 Å − 5680 Å. We found the median signal and noise over the chosen wavelength range on a pixel-by-pixel basis before being fed into the binning code. We required a signal-to-noise ratio of S/N = 50 per stellar bin, which resulted in 111 bins.

The emission line or ionised gas binning was made based on a Hα continuum-subtracted image. We built the image by integrating the data cube flux in a 20 Å window centred on the restframe Hα wavelength of NGC 7130, and by subtracting the integral flux of a window of equal width but 50 Å redwards. The result is shown in Fig. 1. The noise was estimated by quadratically summing the noise values from both the line and the continuum windows. For this binning, we required S/N = 100, which resulted in 2689 bins.

Both stellar and emission line binnings required a minimal single spaxel S/N to be considered for binning. We set this threshold as S/N = 0.1 for stars and S/N = 0.5 for ionised gas. Because some isolated spaxels or small spaxel clusters can be above the threshold without being connected to the main body of regions with a signal above the threshold, we modified GIST not to consider clusters smaller than 500 spaxels for the binning.

For each binned spectrum we also computed a variance spectrum by summing the variances at each wavelength over all the individual spaxels in a bin. The procedure developed to produce Hα continuum-subtracted images was also used to obtain continuum-subtracted images in other lines, such as [O III] λ5007 (Fig. 1).

3.4. Stellar kinematics

GIST recovers the stellar kinematics using pPXF. The latter code fits the spectra for each bin (stellar emission binning in this case) with a linear combination of spectral energy distribution templates convolved with the line-of-sight velocity distribution (LOSVD), which is what one ultimately desires to obtain. We only fitted the two lowest LOSVD momenta, namely the velocity, V_⋆, and the velocity dispersion, σ_⋆.

We used the templates from the E-Miles library³ (Vazdekis et al. 2016) with BaSTI isochrones (Pietrinferni et al. 2004), a Kroupa Universal stellar initial mass function (Kroupa 2001), and the ‘base’ abundances. E-Miles has a spectral resolution of 2.51 Å (FWHM) in the spectral range of interest. The spectral resolution of MUSE as a function of wavelength has been modelled to be (Bacon et al. 2017)

$$\begin{array}{c}\hfill \mathit{FWHM}(\lambda )=5.866\times {10}^{-8}{\lambda }^2-9.187\times {10}^{-4}\lambda +6.040,\end{array}$$(1)

where λ is expressed in Å. The spectral resolution is worse than that of E-Miles on its blue side (λ < 6762 Å) and better on its red side (λ > 6762 Å). GIST supports a wavelength-dependent Gauss convolution of the templates so they match the instrumental spectral resolution. We have implemented a Gauss-convolution of the data so their spectral resolution matches that of the templates at wavelengths where the instrumental resolution is better than that of the templates.

Stellar kinematics were obtained using the wavelength range 4750 Å−7113 Å (4675 Å−7000 Å in the restframe of NGC 7130). We ignored the reddest wavelengths because of the presence of many sky lines, and even though that region contains deep absorption lines that in principle can be used to characterise the stellar kinematics (the Ca II triplet), our experiments showed that using those regions caused the solutions to be noisier. The stellar continuum shape was modelled with an eight-degree additive Legendre polynomial.

The wavelength range affected by the AO laser guide stars (5774 Å−6054 Å) was masked. Sky lines with a peak flux larger than 6 × 10⁻¹⁹ W m⁻² Å⁻¹ arcsec⁻² in the sky line atlas by Hanuschik (2003) were masked by windows with 20 Å in width. We also masked the O₂ B telluric feature at 6862.1 Å − 6964.6 Å (wavelength range from Buton et al. 2013).

Based on our single-component fits in Knapen et al. (2019), we know that emission lines from the galaxy display high velocity dispersion values of up to 400 km s⁻¹, so we masked emission lines with windows of 3000 km s⁻¹ in width. The MUSE observations are so deep that we detected more lines in the innermost arcsecond of NGC 7130 than previous observers, and we resorted to deep studies of another Seyfert galaxy, NGC 1068, for the identification of some of those lines (Koski 1978, Osterbrock & Fulbright 1996; Kraemer & Crenshaw 2000). The masked lines are listed in Table 1.

Our initial guesses for the fit were V_⋆ = cz = 4842 km s⁻¹ (where c is the speed of light in vacuum) and σ_⋆ = 30 km s⁻¹, respectively. The stellar kinematic maps are shown in Fig. 2. The velocity amplitude of the butterfly pattern is between 40 and 50 km s⁻¹, which is similar to what is found for the ionised gas disc (Sect. 4.1). Hence, the circumnuclear stellar population is rotation-supported.

Fig. 2. Left panel: velocity map of the stellar component. The zero in velocity corresponds to z = 0.016221 (cz = 4863 km s−1; see Sect. 4.2). Right panel: velocity dispersion of the stellar component. In this and subsequent maps, the white contour delineates the field of view of the observations.

3.5. One- to six-component fits to emission lines

To study the ionised gas, we started by producing single-component fits of the gas using GIST (akin to those made for Knapen et al. 2019). While examining them with Mapviewer, we discovered that a single-component description does not work well for many regions, especially in the innermost arcsecond. Hence, we wrote our own code that uses pyGandALF to perform multi-component fits with a series of criteria to decide the number of Gaussian components required for a given spectrum.

Increasing the number of gas components comes at the cost of increasing the number of free parameters, so we decided to reduce the complexity of the fit by tying the flux ratios of several doublets to their predicted values as calculated from the Einstein coefficients theoretically estimated by Storey & Zeippen (2000). Those flux ratios are shown in Table 2. For all the lines in a component the kinematics were tied to those of Hα.

We fitted the lines that are used in BPT diagnostics (Sect. 4.6), namely Hβ, [O III] λ4959, [O III] λ5007, [O I] λ6300, [O I] λ6364, [N II] λ6548, Hα, [N II] λ6583, [S II] λ6716, and [S IIλ6731]. Additionally, we fitted the [S III] λ9069 line, which, although in the infrared, is found in a window with no prominent sky lines. We only fitted windows with a width of 5000 km s⁻¹ centred in the lines of interest. These windows contain two sky lines that were masked in the stellar kinematic fit, namely at λ = 6364 Å and λ = 6864 Å, which were left unmasked in here.

Since the ionised gas binning is much finer that the stellar one, the impact of stellar emission is very small and easily close to the noise level. We therefore ignored stellar emission altogether and simply assumed an underlying continuum modelled by a multiplicative eight-order Legendre polynomial (hence ignoring stellar lines). This is further justified by the fact that we only fit a narrow wavelength interval around the emission lines.

After a careful eyeball examination of hundreds of spectra, we decided that up to six Gaussian components per spectrum were required. To help the minimisation within pyGandALF to find the global minimum in multi-component fits, we established a set of initial guesses for the kinematics of the components (Table 3). The need for this is illustrated in Ho et al. (2016). The initial guesses were selected based on experience after manually fitting several hundred spectra representative of the whole data cube. For each fit, we calculated the chi-square:

$$\begin{array}{c}\hfill {\chi }_{j,n,k}^2={\displaystyle \sum _i\frac{{(O_j^i-F_{j,n,k}^i)}^2}{{\left(s_j^i\right)}^2},}\end{array}$$(2)

where O corresponds to the observed spectrum, F corresponds to the fit, s² corresponds to the variances of the spectrum, i is the running index over the unmasked elements of the spectrum, n is the number of components included in the fit, and k labels the fits with the same number of components but different initial guesses. For a bin j and a given number of components, the best fit F_j, n is defined to be the one that yields the smallest chi-square:

$$\begin{array}{c}\hfill {\chi }_{j,n}^2=\underset{k}{min}\left({\chi }_{j,n,k}^2\right).\end{array}$$(3)

We also calculated the chi-square for two restricted ranges in wavelength, namely restframe 4985 Å−5035 Å covering the brightest of the fitted [O III] lines (blue chi-square),

$$\begin{array}{c}\hfill {\chi }_{j,n}^{\mathrm{b}}{}^2={\displaystyle \sum _i\frac{{(O_j^i-F_{j,n}^i)}^2}{{\left(s_j^i\right)}^2}\forall i\mid 4985\mathrm{\AA }<{\lambda }_i<5035\mathrm{\AA },}\end{array}$$(4)

and 6520 Å−6610 Å covering the complex of lines including the [N II] doublet and Hα (red chi-square),

$$\begin{array}{c}\hfill {\chi }_{j,n}^{\mathrm{r}}{}^2={\displaystyle \sum _i\frac{{(O_j^i-F_{j,n}^i)}^2}{{\left(s_j^i\right)}^2}\forall i\mid 6520\mathrm{\AA }<{\lambda }_i<6610\mathrm{\AA }.}\end{array}$$(5)

We used the chi-square values to decide whether adding extra components to a fit improved it enough to justify the growth in complexity. Originally, this was done by calculating the ratio of the chi-square values for n + 1 and n components, ${\chi }_{j,n+1}^2/{\chi }_{j,n}^2$ and comparing it to a threshold ratio of chi-squared values Y_{n + 1, n} chosen so the automatic procedure would, in general, choose as many components as we would if fitting the spectra manually. This method was, for example, implemented by Davis et al. (2012) to choose between single- and two-component fits in IFU data of NGC 1266. However, occasionally, adding an extra component to a fit did not improve it significantly, but adding two caused a large leap in the quality of the fit. Therefore, our final implementation takes into account both Y_{n + 1, n} and Y_{n + 2, n + 1} before considering whether extra components are necessary. Also, since some fit improvements only show in very narrow wavelength ranges that are not well described by the global chi-square value ${\chi }_{j,n}^2$, we also took into account the improvements in the blue and red restricted ranges as quantified by ${\chi }_{j,n+1}^{\mathrm{b}}{}^2/{\chi }_{j,n}^{\mathrm{b}}{}^2$ and ${\chi }_{j,n+1}^{\mathrm{r}}{}^2/{\chi }_{j,n}^{\mathrm{r}}{}^2$ through the thresholds $Y_{n+1,n}^{\mathrm{b}}$, $Y_{n+2,n+1}^{\mathrm{b}}$, $Y_{n+1,n}^{\mathrm{r}}$, and $Y_{n+2,n+1}^{\mathrm{r}}$. In Appendix A.1, we show examples of fits that justify considering these choices. Our complete multi-component fitting strategy is depicted in the flow diagram in Fig. 3.

Fig. 3. Flow diagram illustrating how our code decides the number of components that a bin j requires to have its ionised gas lines fitted.

After comparing automatic fits with manual ones, we found that if we choose the same Y_{n + 1, n}, $Y_{n+1,n}^{\mathrm{b}}$, and $Y_{n+1,n}^{\mathrm{r}}$ values for all bins, we would either be overfitting the bins with a large surface, or underfitting the bins with a small one. A possible cause for that is that the variances of the spectra calculated by the MUSE pipeline might be underestimated in the low-surface-brightness regime, which impacts the χ² determinations. Therefore, we chose different thresholds for bins with fewer than N = 200 spaxels, bins between N = 200 spaxels and N = 999 spaxels, and bins with N = 1000 spaxels or more. The selected values are listed in Table 4. Although these values are to some degree arbitrary, they constitute a reasonable compromise between the need to describe as many visually identified components as possible and the necessity to avoid spurious components that would be overfitting the data.

3.6. Identification of the kinematic components

We found nine distinct kinematic components based on their spatial location, their kinematics, and BPT ratios (see more on the latter criterion in Sect. 4.6). One corresponds to the disc, five are narrow components (σ ≲ 250 km s⁻¹), and three are broad components (σ ≳ 250 km s⁻¹). Since we were only fitting up to six of them per bin, not all are co-spatial. The kinematic properties of these components (as measured after the refining process explained in Sect. 3.7 and Fig. 4) are summarised in Table 5. The refined kinematic maps and the contours of the components are shown in Figs. 5 and 6, respectively. The labels in the top-left corners of the velocity maps identify each component and are assigned a colour that is used consistently throughout this paper. Because the colour scale in Fig. 5 does not allow one to see the details in the kinematic maps of the disc component, we redisplay it in Fig. 7 with an adapted colour scale.

Fig. 4. Flow diagram that illustrates the procedure for rerunning the fit of a spectrum j. This process was done three times after the initial run.

Fig. 5. Kinematic maps of kinematically defined components. They are obtained from fits made accounting for all the emission lines of interest (Sect. 3.6) simultaneously. Nine sets of panels are shown, one for each of the kinematic components identified in the circumnuclear medium of NGC 7130. For each set, the left-hand panel corresponds to the velocity map, and the right-hand panel corresponds to the velocity dispersion map. The velocity colour bar is set so the spacing between velocity steps for |Vgas| > 200 km s−1 is three times larger than for |Vgas| < 200 km s−1. The zero in velocity corresponds to z = 0.016221 (cz = 4863 km s−1; see Sect. 4.2). The labels identifying components are given colours that are used consistently in all the plots in this article.

The disc component has a butterfly pattern, a low velocity dispersion (σ <  100 km s⁻¹), and BPT ratios compatible with star formation (Sect. 4.6). All of the remaining kinematic components, except the zero-velocity narrow component, have line ratios compatible with AGN excitation. The locations of the spiral arms within the bar (Fig. 1) have velocities that are slightly different to those of their surroundings: the north-eastern arm is blueshifted, and the south-western arm is redshifted. We argue that this is evidence for gas inflow through the arms (Sect. 5).

The blueshifted broad component occupies a 3″ − 4″ (900 pc–1200 pc) wide region that runs in the south-east to north-west direction. All of the other blue- and redshifted components are found within the area covered by this component.

The blueshifted narrow component nearly overlaps (in projection) with the blueshifted broad component in the region north-west of the nucleus. The region with the highest velocity relative to the disc of NGC 7130 (darker shades of blue in Fig. 5, or V ∼ −200 km s⁻¹) runs first to the north and then to the north-west. Other regions of this component, especially those to the north-east and to the west of the nucleus, are not as blueshifted and might be to some degree confused with the zero-velocity narrow component. The η parameter, which describes the ionisation mechanism, was used to distinguish them (Sect. 4.6).

The crescent narrow component is found $0\stackrel{″}{.}7$ (200 pc) west of the nucleus. The most redshifted regions of the crescent narrow component have velocities comparable to those in the redshifted narrow component 1. We considered them as two separate components because we see a narrow gap between them (Fig. 6). Also, the ionisation mechanism of this component is different to that of other outflow components, as it has LINER line ratios rather than Seyfert ones (Sect. 4.6).

Fig. 6. Outlines of the nine kinematic components shown in Fig. 5 overlaid on the continuum-subtracted Hα image from Fig. 1. The contours are colour-coded according to the labels identifying each component in Fig. 5. To ease the readability of the figure, we omitted holes within components and regions smaller than N = 50 spaxels disconnected from the main body of the component. The purple x symbols and numbers indicate the bins whose spectra are discussed in Appendix A. Those corresponding to bins in the inner $1\stackrel{″}{.}5\times 1\stackrel{″}{.}5$ are indicated only in the inset to avoid clutter.

The redshifted narrow components 1 and 2 overlap in projection throughout most of their extent. In regions around and south of the nucleus, they correspond to well-defined shoulders in [O III]. At 1″ (300 pc) north of the nucleus, these components are less well defined and might actually be describing the wings of the disc and the blueshifted components. The redshifted broad component is associated with the regions south of the nucleus where the redshifted narrow component 2 is located, and it could be interpreted as a non-Gaussian red wing.

The zero-velocity narrow component is seen in the inner arcsec in clumps north-east, north-west, and south-west of the nucleus (in locations overlapping with the UCNR, see Fig. 1) as well as in the north-eastern arm. Thus, this component correlates with star formation (the line ratios indicate ionisation from star-forming regions; Sect. 4.6). Because it is slightly blueshifted with respect to the disc (V >  −100 km s⁻¹) it could correspond to star-formation driven outflows.

The zero-velocity broad component is found in regions where the redshifted narrow components are also present. As discussed in Appendix A, although a visual examination of the spectra does not reveal an obvious need for this component, some regions require a low-amplitude very broad line for the redshifted narrow components to be properly characterised. Thus, this component might hold no physical meaning, and instead, in many bins, may describe the sum of the effects of non-Gaussian wings of the remaining components. In a few bins it fits with a single Gaussian what a human observer would distinguish as the blueshifted and the redshifted narrow components.

3.7. Refinement of the fits

The velocity and velocity dispersion maps of the kinematic components obtained after the first run were not smooth (they showed large jumps in V and σ between neighbouring bins). Inspired by Davis et al. (2012) and Ho et al. (2016), we refitted the spectra using the median of the previously fitted values for the closest bins (including the bin of interest itself) as initial guesses for V and σ. The distance between bins was calculated as that between geometric centres. The refitting algorithm is described in the flow diagram in Fig. 4. This procedure was done three times to obtain the final maps (Fig. 5). All the results, the discussion, and the conclusions are based on those refined fits.

The median absolute differences in V and σ between the third and the fourth iteration are considered to be representative of the fitting uncertainties. For line fluxes, we estimated median relative uncertainties as

$$\begin{array}{c}\hfill \frac{\delta f_n}{f_n}\equiv \mathrm{median}\left(\frac{|f_{n,i,4}-f_{n,i,3}|}{f_{n,i,4}}\right),\end{array}$$(6)

where the sub-indices 3 and 4 refer to the third and fourth iterations of the fits, respectively, and i refers to the bins where the kinematic component n is defined. The estimated uncertainties are listed in Table 6.

In Appendix A, we show the spectra of a few representative bins indicated by the purple numbers and x signs in Fig. 6. We also discuss their fits.

3.8. The [Fe VII] λ6087 coronal line

We also fitted the [Fe VII] λ6087 coronal line kinematics. This line has a high ionisation potential (99.1 eV; De Robertis & Osterbrock 1984), which unequivocally associates it with AGN photoionisation. Because of its faintness, we treated this line differently from the emission lines discussed above, by giving it its own binning constructed with a methodology similar to that employed for the Hα binning. Here, we used a 20 Å window centred at λ = 6087 Å to characterise the emission, and an equally wide window 50 Å bluewards for the continuum. We required the bins to have S/N = 10. This resulted in 137 bins.

The [Fe VII] line was fitted with a single Gaussian component using GIST applied over the spectral range 6000 Å − 6200 Å. The stellar continuum shape was modelled with a second-degree multiplicative Legendre polynomial. The V and σ maps are shown in Fig. 8.

Within a radius of $0\stackrel{″}{.}2$, we find a bipolar structure centred at the nucleus and oriented in the north-south direction. The northern lobe is blueshifted, whereas the southern one is redshifted. The bipolar structure coincides with that in the single-component fits in Knapen et al. (2019). Beyond the central $0\stackrel{″}{.}2$, the maps are rather noisy, but there are also signs of a blueshifted northern component and a redshifted southern one, possibly corresponding to the blueshifted narrow component and redshifted narrow component 1, respectively.

4. Results

4.1. The orientation of NGC 7130

Disc orientation determinations for NGC 7130 are prone to uncertainties due to the distorted shape of its outer isophotes. The position angle PA = 160° obtained from 2MASS photometry (Skrutskie et al. 2006) is clearly not compatible with the kinematic axes of the stellar and gas discs as seen in Figs. 2, 5 and 7.

Fig. 7. Velocity (left) and velocity dispersion (right) maps of disc component. The data are the same as in the top-left panels in Fig. 5, but the colour scales are adapted to highlight details.

We first attempted to estimate the orientation of NGC 7130 with kinemetry (Krajnović et al. 2006), performing a tilted-ring fit of the velocity map of the disc component (Figs. 5 and 7). To constrain the parameter space, we limited the axis ratio range to 0.7 <  q <  1.0. We found that kinemetry yields a stable position angle at PA = 58 ± 8° (obtained from averaging the fitted orientations for annuli with a semi-axis major axis larger than 2″; uncertainty from the standard deviation), but not a stable axis ratio, which does not allow us to estimate the disc inclination (Fig. 9).

To obtain an estimate of the inclination of NGC 7130, we assumed that the plateau in velocities seen in Fig. 9 corresponds to the maximum circular velocity. By averaging over radii larger than 2″, we found the projected circular velocity to be 42 ± 2 km s⁻¹. Combining this information with the intrinsic luminosity of the galaxy, the Tully-Fisher relation can be used to estimate the inclination.

Because NGC 7130 is rather dusty, we measured its luminosity in the mid-infrared using 3.6 μm science-ready images from the Spitzer Heritage Archive⁴ (Programme ID 90031; PI O. D. Fox). We calculated the total magnitude of the galaxy using the curve of growth technique, following Muñoz-Mateos et al. (2015). For this, we hand-masked stars and galaxies overlapping NGC 7130, then we measured the curve of growth as the cumulative sum of the flux of the galaxy using circular apertures that increase in size by 2″ steps. The magnitude at infinite aperture is where the slope of the curve of growth reaches exactly zero, hence we calculated the gradient of the curve of growth as a function of radius, then fit a line between this gradient and the enclosed magnitude in the outer regions of the curve (for an example, see Fig. 7 of Muñoz-Mateos et al. 2015). The intercept of this line is, by definition, the asymptotic magnitude.

We calculated the uncertainty as both the Poisson error, and the systematic uncertainty due to the local sky subtraction. We estimated the latter by refitting the asymptotic magnitude after subtracting randomised local sky values drawn from a Gaussian distribution with a mean and standard deviation matching those of the original sky measurement. The total magnitude we calculated for this galaxy is m_3.6 μm = 10.82 ± 0.013, where the error is the quadrature sum of the aforementioned random (0.009 mag) and systematic (0.01 mag) errors. Additional uncertainty (not included) comes from our flux interpolation across the relatively large masks, and we also included the AGN flux in this measurement, making this magnitude likely an overestimate.

Assuming a luminosity distance modulus m − M = 34.09 we computed an absolute magnitude M_3.6 μm(AB) = − 23.27. Using the Tully-Fisher relation from Eq. (1) in Sorce et al. (2014), this corresponds to a circular velocity of rotation of v_c = 318 ± 32 km s⁻¹, which is probably an upper limit to the real velocity due to the bias introduced by the AGN. The error bars were calculated using the scatter of 0.43 mag found for the galaxies in Sorce et al. (2014). If 50% of the light came from the AGN, the circular velocity would be v_c = 265 km s⁻¹. As a sanity check, we estimated v_c from the Tully-Fisher relations for Cousins B and R in Pierce & Tully (1992) and adopting the luminosity estimates from Lauberts & Valentijn (1989) of m_B = 12.86 and m_R = 11.57. We obtain v_c = 265 km s⁻¹ and v_c = 281 km s⁻¹, respectively.

If we adopt a round number of v_c ≈ 300 km s⁻¹, this would imply an inclination of i ≈ 8°. Even if the AGN were strongly biasing the Tully-Fisher results and the actual velocity were v_c = 250 km s⁻¹, the resulting inclination would be i ≈ 10°. We can therefore conclude that NGC 7130 is almost face-on.

The galaxy NGC 7130 has an inner pseudo-ring made of clockwise outward-winding spiral fragments encircling the bar (see HST images in Elias-Rosa et al. 2018). Assuming that the spiral fragments are trailing, the near side of NGC 7130 is the southern one.

4.2. The systemic velocity of NGC 7130

The third panel in Fig. 9 shows the systemic velocity as a function of radius for NGC 7130. The inner arcsecond of the galaxy has a velocity of ≈−40 km s⁻¹ with respect to the outer regions. Those outer regions have an average systemic velocity that is 21 ± 2 km s⁻¹ more than what is reported in the NED (cz = 4842 km s⁻¹). We therefore reevaluate the recession velocity to be cz = 4863 km s⁻¹ with respect to the barycentre of the Solar System, which corresponds to a redshift of z = 0.016221. This systemic velocity has been used as the zero in velocity in the velocity maps in Figs. 2, 5, 7, and 8.

Fig. 8. NGC 7130 velocity map (left panel) and velocity dispersion map (right panel) derived from the [Fe VII] λ6087 line. The zero in velocity corresponds to z = 0.016221 (cz = 4863 km s−1; see Sect. 4.2). The colour scales are as in Fig. 5.

Fig. 9. From the upper to the lower panel: position angle, axis ratio, systemic velocity (with respect to z = 0.016151), and circular velocity as functions of radius for the disc component of NGC 7130, as obtained with kinemetry. The dashed horizontal lines show the values averaged over radii larger than 2″.

4.3. Morphology

The inner regions of NGC 7130 host two nuclear rings (Fig. 1), the shape of which has been obtained as described in Comerón et al. (2014). The first one is the tiny elongated UCNR reported by González Delgado et al. (1998) and then confirmed by Knapen et al. (2019). We remeasured it using our Hα continuum-subtracted image and determined a semi-major axis of $0\stackrel{″}{.}68$ (210 pc), an axis ratio of q = 0.59, and a position angle of PA = 163°. The continuum-subtracted Hα image reveals a larger and rounder ring that is not seen in the UV images in González Delgado et al. (1998). This ring has a projected semi-major axis of $3\stackrel{″}{.}0$ (930 pc), an axis ratio of q = 0.77, and a position angle of PA = 119°. The largest of the rings is crossed by Hα emission and dusty spiral arms that reach the innermost region where the UCNR is found. Most of the western side of the ring is missing, which makes it hard to constrain its shape.

We characterised the bar with the 3.6 μm image that we used to estimate the luminosity of NGC 7130 (Sect. 4.1). By running Python’s implementation of ellipse (Jedrzejewski 1987), we found that the bar has its smallest axis ratio, q = 0.48, at a radius of $7\stackrel{″}{.}6$ (2.4 kpc), where the position angle is PA = 6°, which is representative of the whole bar.

We decided against deprojecting the ring parameters, because NGC 7130 is nearly face-on (Sect. 4.1). We found that the difference in position angle between the largest of the rings and the bar is |ΔPA| = 68°. This does not quite match the theoretical picture of rings made of x₂ orbits perpendicular to the x₁ orbits constituting the backbone of the bar (see e.g. Sect. 4.2.1 in Knapen et al. 1995). We note that this expectation may not hold in statistical samples of nuclear rings (Comerón et al. 2010, although the conclusions are also critically dependent on the reliability of the deprojection parameters of the galaxies).

The outermost nuclear ring in NGC 7130 is exceptionally large: its semi-major axis is only 2.5 times smaller than that of the bar (as measured by the radius where it has its maximum ellipticity). This violates the empirical law in Comerón et al. (2010): ‘for barred galaxies, the maximum radius that a nuclear ring can reach is a quarter of the bar radius’. To the best of our knowledge, the only other similar ring is that in ESO 565-11, where the peculiarities have been suggested to be caused by the ring being in an early and fast-evolving stage before it settles into a configuration with a smaller radius (Buta et al. 1999). We note, however, that the other peculiarity of the nuclear ring of ESO 565-11, its large ellipticity, is not shared by that in NGC 7130.

4.4. The [O III] λ5007 surface brightness of the kinematical components

Figure 10 shows the [O III] λ5007 surface brightness maps for each of the nine kinematic components. We display [O III] because it is a tracer of AGN ionisation.

Fig. 10. Surface brightness of [O III] λ5007 emission of the nine kinematically identified components.

The disc component map is very noisy. This is because the [O III] λ5007 disc emission is, in many regions, sub-dominant compared to that from the outflow as seen, for example, in the spectra in Figs. A.1, A.2, A.5, and A.6. Hence, it is poorly constrained by the fits. Hα images of the disc component (not shown) display many details of the spiral arms and the star formation knots in the UCNR.

The [O III] emission of the blueshifted broad component is rather featureless. Its distribution has a nearly circular symmetry in the inner arcsecond (although the maximum emission is not at the galaxy centre, but $0\stackrel{″}{.}1$ or 30 pc to its north). At larger radii, the emission is enhanced in the region north-west of the nucleus where the narrow blueshifted component emission is also enhanced.

The blueshifted narrow component emission is less smooth than that of its broad counterpart. Both share the location of the maximum in emission. The blueshifted narrow component exhibits a well-collimated feature running to the north of the nucleus. This feature coincides with the regions where the component has its largest negative velocities (Sect. 3.6 and Fig. 5) and is also seen in the continuum-subtracted [O III] λ5007 image (Fig. 1).

In Fig. 10, we see further evidence for the separation between the crescent narrow component and the redshifted narrow component 1. Indeed, whereas the redshifted narrow component 1 is centrally concentrated and fades as we move away from the central regions, we see that the crescent narrow component has a sharp border on its western side, where it is the closest to redshifted narrow component 1.

The redshifted narrow components 1 and 2 and the redshifted broad component have their peak in emissions slightly less than $0\stackrel{″}{.}1$ south of the nucleus (roughly symmetric with the peak in the blueshifted components). The maxima in emission of the blue- and redshifted components coincide with the bipolar structure revealed by the coronal gas kinematics (Fig. 8). The redshifted narrow components 1 and 2 were fitted both south and north of the nucleus. We see, however, that the extension of those components to the north-west has a very low surface brightness. Therefore, most of the light in redshifted components comes from regions south of the nucleus (as opposed to the blueshifted light, especially for the narrow component, that comes mostly from north of the nucleus).

4.5. The radio jet

A radio jet as traced by 8.4 GHz continuum emission was observed with the VLA by Thean et al. (2000). Y. Zhao kindly provided us with the re-reduced data presented in Zhao et al. (2016). The radio emission south of the nucleus roughly coincides with the position of the redshifted narrow components 1 and 2. The detached blob of radio emission north of the nucleus is, strikingly, located just north of the collimated blueshifted feature seen in [O III] (Fig. 11).

Fig. 11. [O III] λ5007 surface brightness of blueshifted narrow component (same map as in the corresponding panel in Fig. 10). The black line overlay shows the VLA 8.4 GHz continuum data obtained for Thean et al. (2000) as processed by Zhao et al. (2016). The contour levels cover a range from 0.001 to 0.007 mJy beam−1 in steps of 0.001 mJy beam−1. The angular resolution in radio is $0\stackrel{″}{.}60\times 0\stackrel{″}{.}19$ and the position angle of the point spread function ellipse is ∼10° (Zhao et al. 2016).

4.6. Resolved BPT diagrams: The ionisation mechanisms of the kinematical components

To study the ionisation mechanisms of the kinematic components, we used three different types of the BPT line diagnostics (Baldwin et al. 1981) and produced resolved BPT maps. The main ionisation mechanisms of the components are listed in Table 7.

A caveat of the diagnostics is that they rely on the Hα line. Unfortunately, Hα is often blended with [N II], which makes it very hard to constrain. For example, in many bins the Hα amplitude of outflow components (in particular for the redshifted ones) is fitted as zero, which results in missing bins in the BPT maps. Even when a non-zero Hα amplitude is fitted, its value for components with a low relative Hα surface brightness is prone to large uncertainties (Table 6), which increases the scatter in the horizontal axis in the line diagnostics shown in Figs. 12 and 13.

Fig. 12. Nine sets of two panels are shown, one for each of the kinematic components. For each set, the left-hand panel corresponds to the [O III] λ5007/Hβ versus [N II] λ6563/Hα BPT diagram used to build the resolved BPT map in the right-hand panel. The colour coding indicates the η parameter defined in Erroz-Ferrer et al. (2019). The lines delimiting the regions in the BPT diagram are from Kewley et al. (2001, solid line) and Kauffmann et al. (2003, dashed line).

Fig. 13. Nine sets of two panels are shown, one for each of the kinematic components. For each set, the left-hand panel corresponds to the [O III] λ5007/Hβ versus ([S II] λ6716+[S II] λ6371)/Hα BPT diagram used to build the resolved BPT map in the right-hand panel. The lines delimiting the regions in the BPT diagram are from Kewley et al. (2006).

4.6.1. The [O III] λ5007/Hβ versus [N II] λ6563/Hα line diagnostic

We first discuss the BPT diagnostic that uses the [O III] λ5007/Hβ and the [N II] λ6563/Hα line ratios (Fig. 12). We follow the division of the [O III] λ5007/Hβ-[N II] λ6563/Hα plane by Kewley et al. (2006), where areas below the dashed line (derived semi-empirically by Kauffmann et al. 2003) are considered to be ionised by star formation, and areas above the solid line (derived theoretically by Kewley et al. 2001) are considered to be ionised by the AGN. The region between the two lines is considered to be ionised by a mix of the two. A quantification of the strength of each mechanism is obtained through the parameter η defined by Erroz-Ferrer et al. (2019), so η = −0.5 at the transition between the star-forming and the mixed regimes, and η = 0.5 at the limit between the mixed and the AGN regimes. In the star-forming regime, it is defined to be η <  0.5 and measures the orthogonal distance to the Kauffmann et al. (2003) line. In the AGN regime, η is defined to be η >  0.5 and measures the orthogonal distance to the Kewley et al. (2001) line. In the mixed region, it measures the orthogonal distance to the bisector of the two above-mentioned lines, where η = 0. The concrete formulae to compute η are detailed in Erroz-Ferrer et al. (2019). As indicated in Sect. 3.6, the parameter η was used to decide whether a fitted component corresponded to the blueshifted narrow component (η >  0.8) or the zero-velocity narrow component (η <  0.8). In some bins, both components were present. There, we assigned the component with the largest negative velocity to the blueshifted narrow component, irrespective of η.

We show the BPT maps for the kinematic components in Fig. 12. As in Knapen et al. (2019), the colour-coding traces η, with blue denoting star formation and red indicating AGN ionisation.

The disc is mostly ionised by star formation. In the south-east to north-west axis, there are regions with η >  −0.5, which might indicate the effect of the AGN. We cannot, however, discard the effects of some confusion, especially between the disc and the blueshifted narrow components. Indeed, one of the disc regions with the largest η corresponds to the position of the collimated feature seen in [O III] for the blueshifted narrow component (Sect. 4.4).

The zero-velocity narrow component is also partially ionised by star formation (due partly to how we defined the component). The fraction of a spiral arm delineated by this component (Sect. 3.6) has a very low η, mostly compatible with a pure star-forming ionisation. The innermost clumps corresponding to the UCNR are probably contaminated by the blueshifted narrow component, which explains why η >  −0.5. The remaining components have been ionised by the AGN, as indicated by their orange and red hues in Fig. 12.

4.6.2. The [O III] λ5007/Hβ versus ([S II] λ6716+[S II] λ6371)/Hα line diagnostic

For this line diagnostic, we used the criteria of Kewley et al. (2006). The results are displayed in Fig. 13 and are very similar to those found with the above BPT diagnostic (Fig. 12), namely that the disc and the zero-velocity narrow components are ionised by star formation, whereas the other components are ionised by the AGN. We find that the crescent narrow component straddles the Seyfert and LINER regimes, albeit weighted more strongly towards the LINER side.

We also studied the [O III] λ5007/Hβ versus [O I] λ6300/Hα line diagnostic. Because the results are nearly identical to those in Fig. 13, we do not show them here.

4.7. The electron density, the reddening, and the [S III]/[S II] line ratio

We estimated the mass outflow rates and kinetic powers (Sect. 4.8). To calculate these values the electron density, n_e, and the reddening, E(B − V), of the components must be measured. In our discussion (Sect. 5), we use the [S III]/[S II] line ratio as a proxy for the ionisation parameter (e.g. Díaz et al. 2000). We attempted to produce maps of these three magnitudes for the nine kinematic components. Unfortunately, these maps were very noisy, so we opted to present a single integrated value for each component obtained from the total Hα and Hβ fluxes (for the reddening), the total [S II] λ6716 and [S II] λ6731 fluxes (for the electron density), and the total [S II] and [S III] fluxes (for the ionisation parameter). The uncertainties in the integrated line fluxes are adopted to be those in Table 6. Our estimated n_e, E(B − V), and [S III]/[S II] values are presented in Table 7.

We estimated n_e using the [S II] λ6716/[S II] λ6731 ratio and the parametrisation from Sanders et al. (2016). Unfortunately, the redshifted narrow and broad components are barely detected for these low-amplitude lines (see example spectra in Appendix A). We therefore assumed that they have the same n_e as their blueshifted counterparts.

Extinctions were estimated using the Balmer decrement and assuming a Calzetti et al. (2000) extinction law. Unfortunately, Hα is often blended with the [N II] lines, which makes it very hard to measure its flux, especially for the components with the lowest surface brightness. Thus, for the redshifted broad component, we assumed the extinction to be the same as for the redshifted narrow component 2, with which it overlaps in projection. As expected, the outflow components seen through the disc (the redshifted ones) are more extincted than the ones in front of the disc.

When estimating the [S III]/[S II] line ratio, the flux in [S III] is defined to be the sum of the λ = 9069 Å and the λ = 9532 Å lines, whereas that in [S II] is that of the sum of the λ = 6717 Å and the λ = 6731 Å lines. Since MUSE does not cover [S III] λ9532, we followed Mingozzi et al. (2019) and assumed [S III] λ9532/[S III] λ9069 = 2.5 as theoretically determined by Vilchez & Esteban (1996). This value is similar to the line ratios found by Ramos Almeida et al. (2009) for a sample of five Seyfert galaxies (four of their galaxies have ratios in the 2.2–2.7 range). Because of the large separation in wavelength between the lines, we corrected their fluxes for extinction. We were not able to detect the redshifted broad component in [S III], so we assume that its [S III]/[S II] line ratio is similar to that of its blueshifted counterpart.

In the outflow components, the ionisation parameter is usually larger than in the disc. We find that the crescent narrow component has an extremely low ionisation parameter when compared to the other outflow components ([S III]/[S II] ≈0.1 versus at least 0.9).

4.8. The mass outflow rate and the kinetic power

Here, we followed the procedure from Rose et al. (2018) to calculate the mass outflow rate and the kinetic power of the outflow. We first estimated the mass in each of the kinematic components:

$$\begin{array}{c}\hfill M=\frac{L{(\mathrm{H}\beta )}_{\mathrm{corr}}{m}_p}{{\alpha }_{\mathrm{H}\beta }^{\mathrm{eff}}h{\nu }_{\mathrm{H}\beta }{n}_e},\end{array}$$(7)

where L(Hβ) is the Hβ luminosity corrected for extinction (Table 7), α_Hβ is the effective Case B coefficient (we adopted 3.03 × 10⁻¹⁴ cm³ s⁻¹, which corresponds to an electron temperature T_e = 10⁴ K; Osterbrock & Ferland 2006), m_p is the mass of the proton, and hν_Hβ is the energy of an Hβ photon. To obtain the mass-loss rate, we needed the timescale of the outflow, which can be estimated as

$$\begin{array}{c}\hfill {\tau }_{\mathrm{o}}\sim \frac{R}{V_{\mathrm{o}}},\end{array}$$(8)

where R and V_o are the size and the velocity of the outflow, respectively. For the error propagation, we estimated the uncertainty in τ_o to be of 50%.

If we divide Eq. (7) by Eq. (8), we obtain the mass loss rate

$$\begin{array}{c}\hfill \dot{M}=\frac{L{(\mathrm{H}\beta )}_{\mathrm{corr}}{m}_{\mathrm{p}}{V}_{\mathrm{o}}}{{\alpha }_{\mathrm{H}\beta }^{\mathrm{eff}}h{\nu }_{\mathrm{H}\beta }{n}_{e}R}.\end{array}$$(9)

For V_o, we adopted the velocity of each kinematic component averaged over all bins. The values of R were obtained from the size of the components in Fig. 6. Both are listed in Table 8.

We estimated the kinetic power by accounting for both the net velocity of the outflow, V_o, and its velocity dispersion, σ_o:

$$\begin{array}{c}\hfill {\dot{E}}_{\mathrm{kin}}=\frac{\dot{M}}{2}(V_{\mathrm{o}}^2+3{\sigma }_{\mathrm{o}}^2).\end{array}$$(10)

For σ_o, we adopted the velocity dispersion of each kinematic component averaged over all bins. To gauge the uncertainties in the adopted V_o and σ_o, we estimated the inhomogeneities in each outflow component by finding the dispersion of the velocity and the velocity dispersion, ς(V) and ς(σ), over a kinematic component (Table 8).

The mass-loss rates and kinetic power estimates for each of the outflow components are shown in Table 7. We find a total mass loss rate of Ṁ = 0.55 ± 0.33 M_⊙ yr⁻¹ and a total kinetic power of Ė_kin = (8.8 ± 5.9) × 10⁴⁰ erg s⁻¹ for the blueshifted components. If we also account for the more poorly constrained redshifted components, these values are Ṁ = 1.5 ± 0.9 M_⊙ yr⁻¹ and Ė_kin = (3.4 ± 2.5) × 10⁴¹ erg s⁻¹, respectively. The latter values have to be taken with caution because a third of the mass-loss rate and two thirds of the kinetic power come from the hard-to-characterise redshifted broad component. Indeed, we might only be detecting it over a small fraction of its true extent, which can cause a large overestimate in its associated mass-loss rate and kinetic power (through Eq. (9)). Assuming a bolometric AGN luminosity of L_bol = 2.2 × 10⁴⁴ erg s⁻¹ (Esquej et al. 2014, we corrected the luminosity for the different distance estimates), the fraction of power emitted in kinetic energy is F_kin = 0.040 ± 0.027% for the blueshifted components and F_kin = 0.15 ± 0.11% if accounting for all components.

5. Discussion

In this section, we discuss a plausible interpretation of the complex circumnuclear medium in NGC 7130, in terms of an AGN torus almost perpendicular to the galaxy disc, various red- and blueshifted components of AGN-driven outflow, and star-forming regions. A cartoon depicting our model is presented in Fig. 14.

Fig. 14. Visual representation of our toy model for the ionised circumnuclear gas in NGC 7130. The colours of the lines follow the general colour scheme used throughout this article. Dashed lines indicate features that we hypothesise are obscured by disc dust. The torus is not drawn to scale. The observer’s point of view is indicated by the eye graphic at the top of the image. North is to the right of the cartoon and south to the left.

Whereas NGC 7130 is nearly face-on (Sect. 4.1), the torus of the AGN is probably close to edge-on, as indicated by its Seyfert 1.9 (Véron-Cetty & Véron 2010) or Seyfert 2 (Phillips et al. 1983) classification. In fact, the spectrum of the central spaxel does not show any signature of a broad component associated with the broad-line region in the Hα line, which implies that the broad components are associated with the outflow. We therefore favour the Seyfert 2 classification.

The north-eastern arm of the disc is more blueshifted than its surroundings, and the south-western arm is more redshifted than its surroundings (see Fig. 1 for the location of the arms and Fig. 7 for the kinematics). Given that the southern side of the galaxy is the closest one (Sect. 4.1), this indicates inward streaming motions bringing the gas to the inner few hundreds of parsecs. The disc component velocity map resulting from gas inflow through spiral arms is reminiscent of that for NGC 7213 (Schnorr-Müller et al. 2014).

We find that both the outflow (Fig. 10) and the synchrotron emission (Fig. 11) follow an almost north–south axis. An outflowing bicone close to the plane of the disc and perpendicular to the torus would explain the presence of blueshifted gas on both sides of the nucleus. This would imply that the AGN-driven wind is sweeping material from the disc outwards, including molecular gas. This could be potentially detected with ALMA observations aiming for a CO transition-tracing molecular gas with a lower density than the CO(6–5) line used in Zhao et al. (2016).

The BPT line diagnostics of the outflow components (Figs. 12 and 13) indicate a Seyfert-like ionisation (except for the crescent narrow component, see below). Their typical velocities and velocity dispersions in a single-component fit (V ≈ 100 km s⁻¹ and σ ≈ 150 km s⁻¹, respectively; Knapen et al. 2019) are comparable to those in other Seyfert galaxies (e.g. Müller-Sánchez et al. 2011; Mingozzi et al. 2019) and LIRGs (Arribas et al. 2014; Cazzoli et al. 2014). These pieces of evidence, together with the mass-loss rate and kinetic power estimates discussed below, strongly point to an AGN-driven origin of the outflow.

We find the bipolar emission of the outflow to be asymmetric, with more emission coming from the blueshifted components (Fig. 10). The asymmetry can be either intrinsic or, most likely, a result of absorption by disc dust. The fact that the redshifted outflow components are the most extincted (Table 7) is consistent with the latter possibility. This is in line with previous studies, where AGN line asymmetries have long been interpreted as a sign of an extincted redshifted outflow (Grandi 1977; Heckman et al. 1981). More recently, Crenshaw et al. (2010) and Bae & Woo (2014) established statistically that face-on Seyfert 2 galaxies more often have blueshifted outflow signatures, which was again linked to extinction. This is represented in Fig. 14 with dashed lines indicating plausible locations for completely obscured components.

The blueshifted broad component is seen both north and south of the nucleus. The southern side is less blueshifted than the northern one (Fig. 5). We also observe that, on its northern side, the blueshifted broad component covers a larger solid angle than the blueshifted narrow component. This is represented in Fig. 14 by a larger aperture for the blueshifted broad component. In this model, the southern side is less blueshifted because the axis of the cone is pointing away from the observer. The redshifted broad component would be the (mostly obscured) receding counterpart of the blueshifted broad component.

Why do the narrow components cover a smaller solid angle than the blueshifted broad component? A possibility is that they correspond to different outflows emitted in separate episodes. Another possibility is that the mechanism powering the broad and the narrow components is different. Indeed, the alignment between the narrow components and the synchrotron emission suggests the possibility of a jet-driven component, whereas the broad components could be powered by another mechanism, such as AGN radiation. A third (and complementary) possibility, is that the two components are tracing gas with different properties. The broad components would correspond to lower-density gas, whereas the narrow ones would trace denser gas where efficient energy dissipation would keep the velocity dispersion low. We indeed find that n_e = 300 ± 100 cm⁻³ for the blueshifted broad component, and n_e = 500 ± 100 cm⁻³ for the blueshifted narrow component. Hence, we can say that there is tentative evidence that the blueshifted broad component has a smaller density than the narrow one.

In Knapen et al. (2019), we found a ‘small kinematically decoupled core $0\stackrel{″}{.}2$ (60 pc) in radius ...[that] indicates a tiny inner disc’. This small bipolar structure seen in single-component fits is caused by the combined effects of the two blueshifted components and the redshifted narrow component 1, which dominate the emission in the two blobs just north and just south of the nucleus (Fig. 10), which are also seen in coronal lines (Sect. 3.8). We find that the outflow components that contribute to these knots have Seyfert-like line ratios (Figs. 12 and 13). This, and the fact that coronal lines have large ionisation potentials (≳100 eV; Oke & Sargent 1968), makes it hard to believe that the feature is a nuclear disc. Instead, we propose that we are observing the innermost parts of the outflow.

We estimated the mass outflow rates and kinetic power for each of the outflow components separately. We find that, although both narrow and broad outflow components have virtually the same flux in Hβ (L(Hβ)_corr ∼ (70−80) × 10³⁹ erg s⁻¹), the broad components carry ∼2/3 of the mass outflow and maybe as much as ∼90% of the kinetic power. If we consider only the blueshifted components, which are better constrained due to the reduced extinction, we find that the luminosities of the narrow and the broad components are comparable (L(Hβ)_corr ≈ 30 × 10³⁹ erg s⁻¹ and L(Hβ)_corr ≈ 70 × 10³⁹ erg s⁻¹, respectively). In this case, the broad component carries 90% of the mass loss and almost all the kinetic power (98%). The relatively modest energy output of the ionised gas outflow, compared to the bolometric luminosity of the AGN (F_kin ≈ 0.15% when accounting for all the components, and F_kin ≈ 0.04% when accounting only for the blueshifted components), makes it unlikely that it has a significant impact at a galaxy-wide scale (see discussion in Villar-Martín et al. 2016). The mass-loss rate (Ṁ = 1.5 ± 0.9 M_⊙ yr⁻¹ for all the components and Ṁ = 0.55 ± 0.33 M_⊙ yr⁻¹ for the blueshifted components) is also low compared to the star formation rate, which is estimated to be Ṁ_SFR = 20.93 ± 0.05 M_⊙ yr⁻¹ (Gruppioni et al. 2016) or Ṁ_SFR = 6.7 M_⊙ yr⁻¹ (Diamond-Stanic & Rieke 2012) for the galaxy as a whole, and Ṁ_SFR = 4.3 M_⊙ yr⁻¹ for the innermost kpc (Diamond-Stanic & Rieke 2012).

The estimated mass-outflow rate and the kinetic power are comparable to those computed in several recent studies of AGN by, for example, Villar-Martín et al. (2016) and Rose et al. (2018). The outflow in NGC 7130 appears to be a regular one in the Ṁ versus L_bol and Ė_kin versus L_bol plots in Fig. 1 of Fiore et al. (2017). We note, however, that just as in Villar-Martín et al. (2016) and Rose et al. (2018), we find the mass-outflow rates to be significantly lower than others derived in the literature (e.g., Liu et al. 2013; Harrison et al. 2014; McElroy et al. 2015). Part of the differences come from the fact that, as discussed in Rose et al. (2018), these authors assumed low electron densities on the order of n_e = 100 cm⁻³, which results in an increase of the derived outflow mass (Eq. (7)). Our estimates are also lower than those in Müller-Sánchez et al. (2011), who established their measurements based on coronal lines. On the other hand, they are larger than in Davies et al. (2020), but this is likely because they used a different indicator to measure n_e, which results in larger estimated densities that lower the mass-loss-rate estimate by up to two dex in some cases.

All of the outflow components have Seyfert line ratios, except the crescent narrow component, whose line ratios are compatible with a LINER excitation (Fig. 13). This is demonstrated by the [N II]/Hα and [S II]/Hα line ratios (low-ionisation line ratios; hereafter LILrs) that are larger in the crescent narrow component than elsewhere in the outflow. Assuming a north–south axis for the outflow, the crescent narrow component would be found at the edge of the southern lobe of the ionisation cone (Fig. 6) and may be physically disconnected from the redshifted narrow components 1 and 2 (Sect. 4.4).

The presence of high LILrs at the edges of ionisation cones has been reported by Mingozzi et al. (2019). They also found a correlation between high LILrs and the velocity dispersion of [O III], although this behaviour is not universal. We do not find that the crescent narrow component has a velocity dispersion that is significantly larger than that of the other narrow components (Fig. 5), that is to say, the crescent narrow component average velocity dispersion is σ = 163 km s⁻¹, whereas it is 145 km s⁻¹, 141 km s⁻¹, and 235 km s⁻¹ for the blueshifted narrow component, redshifted narrow component 1, and redshifted narrow component 2, respectively (Table 8). Mingozzi et al. (2019) also found that high LILrs correlate with a low-ionisation parameter, which agrees with the observed [S III]/[S II] line ratio for the crescent narrow component (Sect. 4.7).

Since the velocity dispersion of the crescent narrow component is not particularly large, and because reasonable shock models by Allen et al. (2008) fail to explain line ratios with log ([N II]/Hα) > 0.3 (see discussion in Mingozzi et al. 2019), shock ionisation seems unlikely to explain the peculiar line ratios. Thus, in the case of NGC 7130, the high LILrs could be caused by the radiation field being different from that in other cone regions (Mingozzi et al. 2019). The light reaching the high LILr regions could have been filtered by clumpy ionised clouds between the SMBH and the crescent narrow component.

6. Summary and conclusions

We present a detailed analysis of the circumnuclear medium of NGC 7130, a composite Seyfert 2 galaxy located at a proper distance of 64.8 Mpc. Our work is a follow-up of Knapen et al. (2019). We have used data obtained with MUSE + GALACSI (laser guide star adaptive optics) in narrow-field mode. We achieve an exquisite angular resolution of $0\stackrel{″}{.}17$ which corresponds to ∼50 pc. Figure 1 and the comparison with HST imaging in Knapen et al. (2019) bear witness to the extraordinary quality of the data.

In addition to the well-known UCNR, we found a nuclear ring with a major axis of 3″, which is only 2.5 times smaller than the major axis of the bar (Fig. 1). This nuclear ring is, relative to the bar, exceptionally large.

We analysed the data with adapted tools based on the GIST pipeline. The stellar velocity map (Fig. 2) displays a butterfly pattern with an amplitude of a few tens of km s⁻¹. We used the Tully-Fisher relation to estimate that NGC 7130 is nearly face-on (Sect. 4.1).

We found that for several regions the line profiles are very complex, which is a clear indication of multiple kinematic components in the same line of sight. We have devised an algorithm to automatically fit any number of components between one and six (Figs. 3 and 4). Such multi-component fits were produced for several relevant spectral lines simultaneously (Sect. 3.5). Examples of the fits can be found in Appendix A. We also produced single-component fits of the very weak [Fe VII] λ6087 coronal line (kinematic maps in Fig. 8).

We identified nine distinct kinematic components for the circumnuclear ionised gas. For each of these components we produced kinematic maps (Fig. 5), [O III] maps (Fig. 10), resolved BPT line diagnostic maps (Figs. 12 and 13), the electron density n_e, the extinction E(B − V), and the [S III]/[S II] line ratio (Table 7). Our interpretation of the circumnuclear medium of NGC 7130, in terms of a disc and an AGN-powered outflow, is illustrated in Fig. 14.

The disc of NGC 7130 displays signs of inward streaming motions through the spiral arms (Fig. 7) bringing material to the inner few-hundred parsecs. We also distinguish a zero-velocity narrow component that we interpret as an outflow associated with intense star formation in the UCNR and the spiral arm. The above two components have BPT line ratios compatible with star-forming, or star-forming + AGN, ionisation.

The torus of the AGN is not in the same plane as the disc. The Seyfert 2 spectrum of the central region (at our resolution of ∼50 pc) indicates that the torus is close to edge-on, and therefore almost perpendicular to the plane of the disc. We identify a north–south orientation in the signs of AGN activity: a radio jet and a biconic ionised gas outflow. The lobe length is at least 3″(∼900 pc).

Outflow components on the northern side of the galaxy are dominated by a blueshifted narrow and a blueshifted broad component. The narrow component is much more collimated than the broad one, as indicated by the solid angle covered by each of them. Also, while the narrow component has some well-defined substructure (Fig. 10), the broad one is much smoother. The blueshifted narrow component is probably at least partially jet-powered as indicated by the alignment between its collimated substructure and the radio synchrotron emission (Fig. 11).

To the south of the nucleus, we observe two redshifted narrow components that occupy the same region in projection. A redshifted broad component is also seen in a small blob near the nucleus. The redshifted components are not as extended or bright as the blueshifted ones. Since they are seen through the host’s disc, this is likely due to dust obscuration.

In Knapen et al. (2019), we suggested that the tiny bipolar structure seen in the single-component velocity maps was caused by a nuclear disc. However, its north-south alignment, Seyfert line ratios, and coronal line emission indicate that this feature corresponds to the innermost region of the outflow.

We have measured the ionised gas mass outflow rate and the kinetic power of the outflow (Sect. 4.8). Accounting for all the (the blueshifted only) components, we find them to be Ṁ = 1.5 ± 0.9 M_⊙ yr⁻¹ (Ṁ = 0.55 ± 0.33 M_⊙ yr⁻¹) and Ė_kin = (3.4 ± 2.5) × 10⁴¹ erg s⁻¹ (Ė_kin = (8.8 ± 5.9) × 10⁴⁰ erg s⁻¹), respectively. The kinetic power is F_kin = 0.15 ± 0.11% (F_kin = 0.040 ± 0.027%) of the bolometric AGN output. These values are comparable to those of other AGN (Villar-Martín et al. 2016; Rose et al. 2018) and are roughly a factor of ten lower than the star formation rate. They are probably too low for the outflow to have a galaxy-wide effect. The broad outflow components are responsible for ∼2/3 (∼90%) of the mass outflow rate and about 90% (98%) of the kinetic power output.

All of the outflow components have Seyfert line ratios, except for the crescent narrow component, which has LINER line ratios. This component is located off-axis (presumably at the border of the ionisation cone). We hypothesise that the most likely scenario for its peculiar line ratios is that it is ionised by light that has been filtered by clouds found between the central engine and the crescent narrow component region (see Mingozzi et al. 2019).

Our study has once again proven the extraordinary quality of MUSE data, the finesse of which is such that it fully reveals the complexity of the multi-component outflow. Although multi-component ionised outflows have already been observed (e.g. McElroy et al. 2015), none have, to our knowledge, required the many kinematic components that we have used in our description. This is partly due to the great care taken at scrutinising the data manually before devising an optimal procedure for the automatic multi-Gaussian line fit. It is also because of the great angular resolution provided by the NFM mode of MUSE, which prevents details from being lost due to beam smearing. We conclude that high S/N, high-angular-resolution data, and careful analysis are required to unearth many of the details of AGN outflows. It is to be expected that further AO-assisted IFU observations will show similarly complex ionised outflows in other active galaxies.

Acknowledgments

We thank the anonymous referee for useful comments. We thank Dr. Marja K. Seidel for her contributions to the Very Large Telescope proposal 60.A-9493(A), Dr. Yinghe Zhao for kindly sharing the reduced 8.4 GHz image used in Zhao et al. (2016), Dr. Lodovico Coccato for his help at processing the MUSE data cube, Adrian Bittner for his prompt responses on GIST, Dr. Jesús Falcón-Barroso for pointing to useful references, and Prof. Tom Oosterloo and Dr. Simón Díaz-García for discussions. J.H.K. acknowledges financial support from the European Union’s Horizon 2020 research and innovation programme under Marie Skłodowska-Curie grant agreement No 721463 to the SUNDIAL ITN network, from the State Research Agency (AEI-MCINN) of the Spanish Ministry of Science and Innovation under the grant ‘The structure and evolution of galaxies and their central regions’ with reference PID2019-105602GB-I00/10.13039/501100011033, and from IAC project P/300724, financed by the Ministry of Science and Innovation, through the State Budget and by the Canary Islands Department of Economy, Knowledge and Employment, through the Regional Budget of the Autonomous Community. C.R.A. acknowledges financial support from the Spanish Ministry of Science, Innovation and Universities (MCIU) under grant with reference RYC-2014-15779, from the European Union’s Horizon 2020 research and innovation programme under Marie Skłodowska-Curie grant agreement No 860744 (BID4BEST), from the State Research Agency (AEI-MCINN) of the Spanish MCIU under grants ‘Feeding and feedback in active galaxies’ with reference PID2019-106027GB-C42, ‘Feeding, feedback and obscuration in active galaxies’ with reference AYA2016-76682-C3-2-P, and from IAC project P/301404, financed by the Ministry of Science and Innovation, through the State Budget and by the Canary Islands Department of Economy, Knowledge and Employment, through the Regional Budget of the Autonomous Community. Part of the results are based on public data released from the MUSE NFM-AO commissioning observations at the VLT Yepun (UT4) telescope under Programme ID 60.A-9100, run K. 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 the NASA/IPAC Infrared Science Archive, which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology.

References

Appendix A: Selected examples of multi-component fits

A.1. Example 1: A bin $0\stackrel{″}{.}1$ south of the nucleus with two redshifted narrow components

The bins in the central regions of NGC 7130 show very complex line profiles that require multiple components to be characterised. This is well illustrated in our Example 1 (Fig. A.1). The [O III] λ5007 line profile is a good place to start a visual analysis, because it is isolated and has a high S/N. The [O III] λ5007 line is asymmetric and shows a peak with two redward shoulders, corresponding to the two redshifted narrow components, and a blueshifted tail, corresponding to the blueshifted broad component.

Fig. A.1. Spectrum of Voronoi bin discussed as Example 1 in Appendix A, and indicated by a number ‘1’ in Fig 6. The black line describes the spectrum and the almost coincident purple line describes the fit. The contribution of the six components in this fit are indicated by curves, colour-coded as in Figs. 5 and 6. The shaded areas were masked away while using pyGandALF. The wavelengths are displayed at the rest frame of z = 0.016151.

The peak of the lines in [O III] is displaced blueward with respect to that in Hα and Hβ. This is because the peak is actually made of the superposition of two narrow components, namely the disc component and the blueshifted narrow component. In [O III], the blueshifted narrow component dominates over the disc component, causing the peak of the line to be blueshifted compared to that of the Balmer lines. Fitting the central peak with a single component instead of two slightly shifted ones causes the fit to miss a large fraction of light either in the peaks in the [O III] doublet or in the Balmer lines (i.e. the observed peaks would have a larger amplitude than in the fits).

The redshifted broad component appears as a low-amplitude red tail in both [O III] λ5007 and [N II] λ6583. While it has a small flux, not including it in the fit means that the kinks and shoulders on the red side of [N II] λ6583 are not well described.

A.2. Example 2: A bin at the periphery of the region with two redshifted narrow components

Just outside the region with two fitted redshifted narrow components, subtle signs of redshifted shoulders are seen under visual inspection but have not been fitted because of the χ² limits that we imposed. This is an effect of both the surface brightness decreasing with radius, and also the increasing domination of the disc component with increasing distance from the centre, which ultimately causes the traces of redshifted components to be hidden within the disc component wings. Allowing for less stringent ratios of χ² limits makes the region where these two components are fitted larger, at a cost of introducing spurious components everywhere in the field of view.

In Fig. A.1, we see that the kinks caused by the redshifted narrow components are observed at a higher surface brightness in [O III] than in [N II] λ6583. As a consequence, in noisy bins, clear traces of one or two of these narrow components are seen in [O III] only. An example of that is the spectrum displayed in our Example 2 (Fig. A.2), where the red side of [N II] λ6583 barely shows any signs of the two redshifted narrow components. A global χ² criterion does, in this case, easily overlook the clear indications for the redshifted components seen in [O III]. Hence, this kind of spectrum was the motivation to introduce a χ² criterion referring to the narrow spectral window around [O III] λ5007 (see Sect. 3.5).

Fig. A.2. Same as Fig. A.1 but for bin indicated by number ‘2’ in Fig. 6. This fit is discussed as Example 2 in Appendix A.

A.3. Example 3: A bin in the region with the crescent narrow component

The bins in this region have a conspicuous redshifted shoulder in their [N II] λ6583 line that is also seen, albeit less clearly, in the [O III] lines. The example displayed in Fig. A.3 is an exceptionally good case, but at the periphery of the crescent, the redshifted shoulder of [N II] λ6583 is also obvious. In these periphery bins, the global χ² ratio criterion is not sensitive enough to capture that crescent narrow component. To solve that without loosening the χ² ratio criteria too much, we included the χ² ratio criteria for the wavelength range in the line complex including [N II] and Hα (Sect. 3.5).

Fig. A.3. Same as Fig. A.1 but for the bin indicated by number ‘3’ in Fig 6 and with three components. This fit is discussed as Example 3 in Appendix A.

A.4. Example 4: A bin $0\stackrel{″}{.}1$ north of the nucleus with a conspicuous blueshifted broad component

The innermost bins on the northern side of the nucleus show a clear blueshifted broad shoulder in the [O III] line profiles (Fig. A.4). Hints of the two redshifted narrow components observed south of the nucleus (Example 1) are also seen here, albeit with a smaller amplitude.

Fig. A.4. Same as Fig. A.1 but for the bin indicated by number ‘4’ in Fig 6. This fit is discussed as Example 4 in Appendix A.

The line profile in the [N II] plus Hα complex is much more complicated to interpret because of the blending of lines. A small bump appears almost midway between Hα and [N II] λ6583. The peak of this bump is found at λ = 6571.5 Å, which would correspond to a Hα line redshifted by roughly 400 km s⁻¹, or to an [N II] λ6583 line blueshifted by about 500 km s⁻¹ (with quite some imprecision because of the low contrast of the line). In [O III], we see that both the blueshifted broad component and the redshifted narrow component 2 are shifted by several hundreds of kilometres per second with respect to the disc, hence both could contribute to the bump between Hα and [N II] λ6583.

In practice, if the bump is an [N II] λ6583 line, it is compatible with the blue wing in [N II] λ6548, but a bit too blueshifted to be explained by the component that makes the blueshifted wing in [O III]. What pyGandALF does in this case is transform the component that is set according to the initial conditions as a blueshifted narrow component into the blueshifted broad component, and it takes the component with initial conditions of a blueshifted broad component and blueshifts it a bit more in order for its velocity to be compatible with that of the bump. Thus, the resulting fit has two blueshifted broad components, but the bluest one is too blue to be compatible with the shape of the [O III] lines, so its amplitude is fitted to have a very low value there. In Fig. 5, the most blueshifted component is plotted as the blueshifted broad component, and the least blueshifted one is plotted as the narrow blueshifted component where it can be distinguished from the genuinely narrow component because of its yellow colour, corresponding to σ ∼ 350 km s⁻¹, in the velocity dispersion map. The patch where this occurs is very small (less than $0\stackrel{″}{.}2$ across), and hence it corresponds to an angularly unresolved region.

An alternative to the solution with two blueshifted broad components proposed by pyGandALF would be the bump to be related to the redshifted narrow component 2. However in Figs. A.1 and A.2, the amplitude of this component in Hα is small compared to that of [O III] λ5007. If this is true everywhere for this component, it would be insufficient to explain the bump. On the other hand, the solution offered by pyGandALF poses the problem of why the second blueshifted broad component is not seen in [O III].

In some bins, such as the one discussed here, the redshifted broad component is replaced by a zero-velocity broad component. Even though its addition adheres to our number of component selection criteria, this component is the only one in our fits that does not seem justified upon visual inspection. Including it might be overfitting the data. However, making the criteria more stringent to eliminate the need for this component also causes components that are clearly identified elsewhere under visual inspection to be ignored, so we decided against it. In the particular case shown here, the zero-velocity broad component has a low amplitude. However, in some cases, it has a significant amplitude in [O III]. In these cases the zero-velocity broad component is fitting all the blueshifted and redshifted components with a single very broad Gaussian.

A.5. Example 5: A bin 1″ north-west of the nucleus with a clear blueshifted narrow component

The bin whose spectrum is shown in Fig. A.5 belongs to a region to the north-west of the nucleus where, except for the [O III] doublet, the disc component and the blueshifted narrow component have similar contributions to the line flux. Thus, the blueshifted narrow component does not appear as a blue wing of the disc component as it does in Examples 1, 2, and 3 (Figs. A.1–A.3), but instead the lines are double-peaked. In [O III], the blueshifted narrow component clearly dominates, as seen by comparing the positions of the peaks with those of the vertical marks below the [O III] label that show the positions of lines at rest at z = 0.016151.

Fig. A.5. Same as Fig. A.1 but for the bin indicated by number ‘5’ in Fig 6 and with three components. This fit is discussed as Example 5 in Appendix A.

In this bin, the difference in velocity between the narrow and blueshifted broad components is smaller than 50 km s⁻¹, maybe denoting a single blueshifted component with a core and broad wings. This kind of relatively simple spectrum (only three obvious components) proved surprisingly challenging to fit. Indeed, moving from one to two components did not improve the χ² substantially. A third component was required to obtain a significantly better fit than with a single component. This is the sort of spectrum that motivated us to consider both the Y_{n + 1, n} and Y_{n + 2, n + 1} chi-squared ratios when deciding whether to add an extra component to an n-component fit (Sect. 3.5).

Our algorithm has fitted the redward wings of the [O III] doublet with a low-amplitude redshifted narrow component, which is probably unnecessary. However, adopting more stringent thresholds to add a fourth component would result in visually identified components not being fitted elsewhere.

A.6. Example 6: A bin 1″ north-west of the nucleus with a subtler blueshifted narrow component

The case shown in Fig. A.6 is similar to Example 5, but with a much subtler blueshifted narrow component, which could be confused with a blueward extension of an asymmetric disc component if [O III] were not considered. Indeed, the peak of [O III] is blueshifted compared to the restframe of the galaxy, which indicates that the blueshifted narrow component is genuine. On this occasion, the blueshifted broad component is clearly more blueshifted than the blueshifted narrow component.

Fig. A.6. Same as Fig. A.1 but for the bin indicated by number ‘6’ in Fig 6 and with three components. This fit is discussed as Example 6 in Appendix A.

A.7. Example 7: A bin $2\stackrel{″}{.}5$ north-east of the nucleus with the zero-velocity narrow component

In some regions, the disc component lines have low-amplitude wings that are not well fitted by a single Gaussian function. In this case, the Hα and the [N II] λ6583 lines show the need for a second, slightly blueshifted component (Fig. A.7). The existence of this component (the zero-velocity narrow component) cannot be confirmed in other lines, because, if present, it is buried in the noise.

Fig. A.7. Same as Fig. A.1 but for the bin indicated by number ‘7’ in Fig 6 and with two components. This fit is discussed as Example 7 in Appendix A.

A.8. Example 8: A bin 3″ south-west of the nucleus with only the disc component

Far from the nucleus, especially in regions to its north-east and its south-west, bins with a single component are observed (Fig. A.8). While this lack of complexity might be partly caused by the low surface brightness of these areas, which prevents detections of low-amplitude components, this cannot fully explain the simplicity of these spectra. Indeed, if we use the area of the bins as a proxy for the surface brightness (the lower the surface brightness, the larger the Voronoi bin), then we find bins with a disc plus a blueshifted broad component with sizes comparable to those of the bins showing the disc component only.

Fig. A.8. Same as Fig. A.1 but for the bin indicated by number ‘8’ in Fig 6 and with one component. This fit is discussed as Example 8 in Appendix A.

All Tables

All Figures

Fig. 1. Left panel: white-light image of the centre of NGC 7130 obtained from integrating the reduced MUSE data cube along the spectral direction. Middle panel: continuum-subtracted Hα image obtained from the same data cube. The red ellipses indicate the outline of the nuclear rings (see Sect. 4.3). Right panel: continuum-subtracted [O III] λ5007 image. The intensities are scaled logarithmically. In this and all subsequent maps, the plus sign indicates the inferred galaxy centre and the insets show an enlarged version of the centremost region of the galaxy. The axes are in arcseconds. North is up and east is left.

In the text

	Fig. 2. Left panel: velocity map of the stellar component. The zero in velocity corresponds to z = 0.016221 (cz = 4863 km s−1; see Sect. 4.2). Right panel: velocity dispersion of the stellar component. In this and subsequent maps, the white contour delineates the field of view of the observations.
In the text

	Fig. 3. Flow diagram illustrating how our code decides the number of components that a bin j requires to have its ionised gas lines fitted.
In the text

	Fig. 4. Flow diagram that illustrates the procedure for rerunning the fit of a spectrum j. This process was done three times after the initial run.
In the text

Fig. 5. Kinematic maps of kinematically defined components. They are obtained from fits made accounting for all the emission lines of interest (Sect. 3.6) simultaneously. Nine sets of panels are shown, one for each of the kinematic components identified in the circumnuclear medium of NGC 7130. For each set, the left-hand panel corresponds to the velocity map, and the right-hand panel corresponds to the velocity dispersion map. The velocity colour bar is set so the spacing between velocity steps for |Vgas| > 200 km s−1 is three times larger than for |Vgas| < 200 km s−1. The zero in velocity corresponds to z = 0.016221 (cz = 4863 km s−1; see Sect. 4.2). The labels identifying components are given colours that are used consistently in all the plots in this article.

In the text

Fig. 6. Outlines of the nine kinematic components shown in Fig. 5 overlaid on the continuum-subtracted Hα image from Fig. 1. The contours are colour-coded according to the labels identifying each component in Fig. 5. To ease the readability of the figure, we omitted holes within components and regions smaller than N = 50 spaxels disconnected from the main body of the component. The purple x symbols and numbers indicate the bins whose spectra are discussed in Appendix A. Those corresponding to bins in the inner $1\stackrel{″}{.}5\times 1\stackrel{″}{.}5$ are indicated only in the inset to avoid clutter.

In the text

	Fig. 7. Velocity (left) and velocity dispersion (right) maps of disc component. The data are the same as in the top-left panels in Fig. 5, but the colour scales are adapted to highlight details.
In the text

	Fig. 8. NGC 7130 velocity map (left panel) and velocity dispersion map (right panel) derived from the [Fe VII] λ6087 line. The zero in velocity corresponds to z = 0.016221 (cz = 4863 km s−1; see Sect. 4.2). The colour scales are as in Fig. 5.
In the text

	Fig. 9. From the upper to the lower panel: position angle, axis ratio, systemic velocity (with respect to z = 0.016151), and circular velocity as functions of radius for the disc component of NGC 7130, as obtained with kinemetry. The dashed horizontal lines show the values averaged over radii larger than 2″.
In the text

	Fig. 10. Surface brightness of [O III] λ5007 emission of the nine kinematically identified components.
In the text

Fig. 11. [O III] λ5007 surface brightness of blueshifted narrow component (same map as in the corresponding panel in Fig. 10). The black line overlay shows the VLA 8.4 GHz continuum data obtained for Thean et al. (2000) as processed by Zhao et al. (2016). The contour levels cover a range from 0.001 to 0.007 mJy beam−1 in steps of 0.001 mJy beam−1. The angular resolution in radio is $0\stackrel{″}{.}60\times 0\stackrel{″}{.}19$ and the position angle of the point spread function ellipse is ∼10° (Zhao et al. 2016).

In the text

Fig. 12. Nine sets of two panels are shown, one for each of the kinematic components. For each set, the left-hand panel corresponds to the [O III] λ5007/Hβ versus [N II] λ6563/Hα BPT diagram used to build the resolved BPT map in the right-hand panel. The colour coding indicates the η parameter defined in Erroz-Ferrer et al. (2019). The lines delimiting the regions in the BPT diagram are from Kewley et al. (2001, solid line) and Kauffmann et al. (2003, dashed line).

In the text

	Fig. 13. Nine sets of two panels are shown, one for each of the kinematic components. For each set, the left-hand panel corresponds to the [O III] λ5007/Hβ versus ([S II] λ6716+[S II] λ6371)/Hα BPT diagram used to build the resolved BPT map in the right-hand panel. The lines delimiting the regions in the BPT diagram are from Kewley et al. (2006).
In the text

Fig. 14. Visual representation of our toy model for the ionised circumnuclear gas in NGC 7130. The colours of the lines follow the general colour scheme used throughout this article. Dashed lines indicate features that we hypothesise are obscured by disc dust. The torus is not drawn to scale. The observer’s point of view is indicated by the eye graphic at the top of the image. North is to the right of the cartoon and south to the left.

In the text

Fig. A.1. Spectrum of Voronoi bin discussed as Example 1 in Appendix A, and indicated by a number ‘1’ in Fig 6. The black line describes the spectrum and the almost coincident purple line describes the fit. The contribution of the six components in this fit are indicated by curves, colour-coded as in Figs. 5 and 6. The shaded areas were masked away while using pyGandALF. The wavelengths are displayed at the rest frame of z = 0.016151.

In the text

	Fig. A.2. Same as Fig. A.1 but for bin indicated by number ‘2’ in Fig. 6. This fit is discussed as Example 2 in Appendix A.
In the text

	Fig. A.3. Same as Fig. A.1 but for the bin indicated by number ‘3’ in Fig 6 and with three components. This fit is discussed as Example 3 in Appendix A.
In the text

	Fig. A.4. Same as Fig. A.1 but for the bin indicated by number ‘4’ in Fig 6. This fit is discussed as Example 4 in Appendix A.
In the text

	Fig. A.5. Same as Fig. A.1 but for the bin indicated by number ‘5’ in Fig 6 and with three components. This fit is discussed as Example 5 in Appendix A.
In the text

	Fig. A.6. Same as Fig. A.1 but for the bin indicated by number ‘6’ in Fig 6 and with three components. This fit is discussed as Example 6 in Appendix A.
In the text

	Fig. A.7. Same as Fig. A.1 but for the bin indicated by number ‘7’ in Fig 6 and with two components. This fit is discussed as Example 7 in Appendix A.
In the text

	Fig. A.8. Same as Fig. A.1 but for the bin indicated by number ‘8’ in Fig 6 and with one component. This fit is discussed as Example 8 in Appendix A.
In the text