# Composite Bulges – IV. Detecting Signatures of Gas Inflows in the IFU data:
The MUSE View of Ionized Gas Kinematics in NGC 1097

Tutku Kolcu1,2, Witold Maciejewski1, Dimitri A. Gadotti2,3, Francesca
Fragkoudi2,4, Peter Erwin5,6, Patricia Sánchez-Blázquez7, Justus Neumann8,
Glenn Van de Ven8,9, Camila de Sá-Freitas2, Steven Longmore1, Victor P.
Debattista10
1 Astrophysics Research Institute, Liverpool John Moores University, IC2
Liverpool Science Park, 146 Brownlow Hill, L3 5RF, UK
2 European Southern Observatory, Karl-Schwartzchild-Str. 2, D-85748 Garching
bei München, Germany
3 Centre for Extragalactic Astronomy, Department of Physics, Durham
University, South Road, Durham DH1 3LE, UK
4 Institute for Computational Cosmology, Department of Physics, Durham
University, DH1 3LE, United Kingdom
5 Max-Planck-Insitut für extraterrestrische Physik, Giessenbachstrasse, 85748
Garching, Germany
6 Universitäts-Sternwarte München, Scheinerstrasse 1, D-81679 München, Germany
7 Departamento de Física Teórica, Universidad Autónoma de Madrid, E-28049
Cantoblanco, Spain
8 Max-Planck Institut für Astronomie, Königstuhl 17, D-69117, Heidelberg,
Germany
9 Department of Astrophysics, University of Vienna, Türkenschanzstrasse 17,
1180 Vienna, Austria
10 Jeremiah Horrocks Institute, University of Central Lancashire, Preston, PR1
2HE, UK
Contact e-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Using VLT/MUSE integral-field spectroscopic data for the barred spiral galaxy
NGC 1097, we explore techniques that can be used to search for extended
coherent shocks that can drive gas inflows in centres of galaxies. Such shocks
should appear as coherent velocity jumps in gas kinematic maps, but this
appearance can be distorted by inaccurate extraction of the velocity values
and dominated by the global rotational flow and local perturbations like
stellar outflows. We include multiple components in the emission-line fits,
which corrects the extracted velocity values and reveals emission associated
with AGN outflows. We show that removal of the global rotational flow by
subtracting the circular velocity of a fitted flat disk can produce artefacts
that obscure signatures of the shocks in the residual velocities if the inner
part of the disk is warped or if gas is moving around the centre on elongated
(non-circular) trajectories. As an alternative, we propose a model-independent
method which examines differences in the LOSVD moments of H$\alpha$ and [N
II]$\lambda$6583\. This new method successfully reveals the presence of
continuous shocks in the regions inward from the nuclear ring of NGC 1097, in
agreement with nuclear spiral models.

###### keywords:

galaxies: general - galaxies: star formation - galaxies: kinematics and
dynamics - shock waves -

##  1 Introduction

Gas inflows play an important role in the evolution of galaxies: they enhance
nuclear star formation (SF) (Ellison et al., 2011; Prieto et al., 2019), which
can cause starbursts in nuclear rings (Heller & Shlosman, 1994; Mazzuca et
al., 2008) and change the morphology of the central regions by creating
substructures (e.g., gaseous and stellar disks, rings) (Shlosman et al., 1989;
Athanassoula, 1992; Martini et al., 2003; Kim et al., 2012; Sormani et al.,
2015). Further, the gas reaching the very centre can trigger active galactic
nuclei (AGN) (Davies et al., 2014; Audibert et al., 2020) and lead to highly
energetic feedback (Prieto et al., 2005; Davies et al., 2014; Audibert et al.,
2019). Outflowing gas from AGN feedback can in turn rapidly suppress SF and
cause quenching (Man et al., 2019; Donnari et al., 2021). Hence, to understand
the structural evolution of the galaxies, it is important to understand the
role of the inflowing gas.

Mechanisms triggering gas inflow on various scales have been a prominent focus
of many observational and theoretical studies (Sanders & Huntley, 1976;
Athanassoula, 1992; Emsellem et al., 2003; Maciejewski, 2004; Fathi et al.,
2006; Lin et al., 2013; Fragkoudi et al., 2016; Prieto et al., 2019).
Simulations have revealed that non-axisymmetric potentials such as large-scale
bars and ovals efficiently transport gas towards the inner kiloparsec
(Athanassoula, 1992; García-Burillo et al., 2005; Lin et al., 2013; Audibert
et al., 2021) by generating torques that put gas clouds on intersecting
trajectories. This leads to shocks, in which the gas loses angular momentum
and funnels inwards, though the efficiency of this mechanism inside the inner
kiloparsec is questioned (Kim et al., 2012). It is also evident from both
observations and simulations that central stellar components such as nuclear
bars and disks affect the kinematics of the gas (Shlosman, 2001; Maciejewski
et al., 2002; Fathi et al., 2005), but we require more understanding of
mechanisms transferring gas in the nuclear regions of the galaxies.

Shocks in gas that form in response to the gravitational torques of bars and
ovals extend over large distances along straight or curved lines, the latter
often assuming a spiral shape. Although dust filaments can indicate gas
compression in shocks, shock condition implies a rapid change in velocity
(Athanassoula, 1992; Maciejewski, 2004), therefore coherent structures in
kinematic maps are best indicators of shocks (Zurita et al., 2004; Fathi et
al., 2005; Storchi-Bergmann et al., 2007). In this series of papers, we aim to
identify gas inflows by searching for coherent kinematic structures in
integral-field spectroscopic observations, which allow us to obtain spectral
information over a two-dimensional spatial field-of-view (FoV) of integral
field units (IFUs).

We will conduct our search using a diverse sample of galaxies from a multi-
wavelength imaging and spectroscopic observing campaign Composite Bulges
Survey (CBS, Erwin et al., 2023). CBS focuses on understanding the formation
and growth mechanisms of bulges, and the correlation with their host galaxy
properties and aims to identify and characterise different forms of nuclear
stellar components, from large-scale bars down to nuclear scale
discs/rings/bars and nuclear star clusters, of nearby galaxies (Erwin et al.,
2021). The sample is a mass- and volume-limited set of 53 disk (S0–Sbc)
galaxies with optical and near-IR (NIR) imaging with the Hubble Space
Telescope (HST), of which 40 have been observed with the Multi Unit
Spectroscopic Explorer (MUSE) instrument of the Very Large Telescope (VLT).
Therefore, benefiting from a diverse galaxy sample, we should be able to
determine whether gas inflows in the centres of galaxies occur predominantly
in coherent kinematic structures. Moreover, we should be able to examine the
connection between different central stellar structures and the presence of
coherent kinematic structures indicating inflow.

Tracing the flow of gas in nuclear regions of galaxies has been a challenging
task with former generations of IFUs because of their instrumental
limitations, such as small fields of view and low spatial resolution.
Moreover, in the case of having a low signal-to-noise ratio (SNR), coherent
kinematic structures were lost after binning (Fathi et al., 2005, 2006;
Storchi-Bergmann et al., 2007). In our study, we predominantly rely on data
from MUSE, which with its large, 60″$\times$60″, FoV, fine spatial resolution
and 0.2″$\times$0.2″ sampling, allows us to capture the nuclear regions of
nearby galaxies in great detail. Moreover, MUSE can obtain sufficient SNR for
nearby targets, minimising the need for binning to improve SNR. For our study
on gaseous kinematic structures, we prefer not to bin the data to prevent
losing spatial information. Hence, in comparison to pioneering studies that
suffered from instrumental limitations, advancements of MUSE maximise the
prospect of finding coherent kinematic structures.

In this pilot study, we introduce the methodology that we developed to
identify coherent kinematic structures in the centres of galaxies, and apply
it to the nuclear regions of a galaxy in the CBS sample, NGC 1097 (D=14.8 Mpc,
Freedman et al., 2001), 1″ $\approx$ 72 pc). NGC 1097 is an SBb galaxy with a
strong bar at position angle $\sim$141°(Prieto et al., 2005) and it contains a
prominent star-forming nuclear ring in the inner $\sim$8″–10″(Bittner et al.,
2020). The galaxy hosts a LINER111LINER: Low Ionization Narrow Emission-line
Regions, can be distinguished from the Seyfert type nuclei based on their [O
III] luminosity and kinematics (Kewley et al., 2006) type, low-luminosity
active galactic nucleus (LLAGN – $L_{bol}\lesssim 10^{42}\,erg\,s^{-1}$, Ho,
2008) with bolometric luminosity $L_{bol}=8.6\times 10^{41}\,erg\,s^{-1}$
(Nemmen et al., 2006). The galaxy is interacting with a companion elliptical
galaxy NGC 1097A located to its northwest. Because of its complex
characteristics, NGC 1097 has been a prominent focus of many studies (Quillen
et al., 1995; Prieto et al., 2005; Fathi et al., 2006; Davies et al., 2009;
van de Ven & Fathi, 2010; Lin et al., 2013; Bowen et al., 2016; Izumi et al.,
2017; Prieto et al., 2019; Legodi et al., 2021).

The flow of the paper is as follows: In Section 2 we present the observations
of NGC 1097, and we introduce the data extraction processes for MUSE,
including our approach of fitting the gas emission with multiple components.
We present our main methodology of finding coherent kinematic structures, and
we apply it to the galaxy NGC 1097 in Sect. 3. We provide the discussion of
our findings in Sect. 4.

##  2 Observations and Procedures for deriving the LOSVD moments

###  2.1 Observations of NGC 1097

Figure 1: MUSE coverage of NGC 1097. left: Large-scale (360″$\times$360″)
image of NGC 1097 taken by VIMOS/VLT. The MUSE coverage of 60″$\times$60″ is
highlighted with the red box. right: Colour composite image within the MUSE
FoV, constructed by using the MUSE data cube. Images in both panels are
centred on the brightest point in the centre. The North is to the top and the
East is to the left.

MUSE observations of NGC 1097 used in CBS are provided by the TIMER project as
a part of observing proposal ID:097.B-0640 (PI: Dimitri Gadotti), carried out
in ESO Period 97 on 31st of July 2016. The TIMER project studies nearby barred
galaxies ($D<40Mpc$) larger than 1 arcmin with stellar masses above
$10^{10}M_{\odot}$, brightness threshold of 15.5 B-mag and inclinations
smaller than $\approx$60°, which display prominent central components (e.g.,
nuclear rings/discs, inner bars). Observations are done in MUSE Wide-Field
Mode (WFM), for a nominal wavelength range of $4700-9350$Å. The FoV of WFM
encloses an area of 60″$\times$60″, has a spatial sampling of 0.2″/pixel
(Bacon et al., 2010). The point spread function full width at half maximum of
the observations for NGC 1097 is 0.8″, which corresponds to a spatial
resolution of $\sim$58 pc. Further details on observations and the data
reduction processes, products including stellar kinematics, ages and
metallicities of the sample can be found in Gadotti et al. (2019, 2020) and
Bittner et al. (2020). In the left panel Fig. 1, we present a large-scale
(360″$\times$360″) image222Image credits: ESO Press Release, 2004, date of
observation: 09-10/12/2004, image ID: eso0438d, media source:
https://www.eso.org/public/images/eso0438d/ of NGC 1097, which reveals the
spiral arms on the north and the south, and the companion elliptical galaxy
NGC 1097A in the northwest. In the right panel Fig. 1, we present the colour
composite image within the inner 60″$\times$60″, constructed from the MUSE
data. In both panels, the stellar nuclear ring in the centre and the
connection of dust lanes to the nuclear ring are very well resolved.

###  2.2 Extraction of the kinematic data

Figure 2: Our implementation of the Data Analysis Pipeline. Arrows follow the
order of the workflow within the pipeline. Voronoi binned data is transferred
to stellar kinematics, extinction (not used), stellar population (not used)
and emission line analysis modules respectively. Note that the emission line
analysis can also be performed on single spaxels. In the rightmost blocks, we
present the routine used in each module to perform the analysis alongside the
derived products.

Spatially resolved spectroscopic data hold a great amount of information,
which can be used to probe the fundamental properties of galaxies. To derive
emission line fluxes and to extract stellar and ionised gas kinematics from
MUSE data, we employ the Data Analysis Pipeline333DAP can be accessed at:
https://gitlab.com/francbelf/ifu-pipeline (DAP) of the PHANGS collaboration
(Emsellem et al., 2022), which is an implementation based on the GIST pipeline
(Bittner et al., 2019). DAP is an exceptionally functional software of modular
structure, and it employs the commonly used penalised pixel-fitting (pPXF)
code (Cappellari & Emsellem, 2004) to perform both stellar and emission line
analysis. To optimize the processing time of DAP, we discarded the stellar
population history module, since it was outside of the scope of this work. A
schematic overview of DAP highlighting the main modules, the routines for data
extraction and the products is presented in Fig. 2. A representative spectrum
with the pPXF fit for both stellar and emission lines is shown in Fig. 5 (top
panel).

####  2.2.1 Deriving stellar line-of-sight velocity distribution moments

We extracted the stellar kinematics in the form of line-of-sight velocity
distribution (LOSVD) moments — mean velocity (V), velocity dispersion
($\sigma$) and Gauss-Hermite higher-order moments ($h_{3}$ & $h_{4}$) (van der
Marel & Franx, 1993). Before deriving the LOSVD moments, the data is re-
sampled on a logarithmic wavelength axis, then it is spatially binned with
adaptive weighted implementation (Diehl & Statler, 2006) of the Voronoi
tessellation (Cappellari & Copin, 2003). For the spectral fitting, we limited
the rest frame wavelength range to 4800-7000Å. To fit the stellar continuum
with pPXF, we used the EMILES single-stellar population (SSP) model library
(Vazdekis et al., 2015). SSP models convolved with the LOSVD are fitted to
each spectrum and the best-fitting parameters are determined by $\chi^{2}$
minimisation in pixel space. We used a threshold SNR of 3 for spaxels included
in Voronoi bins, to prevent contribution from low-quality spaxels whose signal
is below the isophote level. We also used a target SNR of 40 for each Voronoi
bin, as this value is commonly used in literature and provides a balance
between the accuracy of the extracted measurements and the spatial resolution
(van der Marel & Franx, 1993).

When fitting the stellar kinematics, we used 8th order multiplicative Legendre
polynomials within the pPXF routine and included no additive polynomials.
Multiplicative polynomials adjust the inaccuracies in the spectral calibration
and make the fits independent of reddening by dust so that we do not need to
provide a particular reddening curve within the routine (Cappellari, 2017).
Although additive polynomials can minimise mismatches between stellar
templates and absorption lines and correct impurities from the sky subtraction
(Cappellari, 2017), they can also alter the equivalent width of the Balmer
absorption lines and might introduce non-physical correction to the line
fluxes (Emsellem et al., 2022). Therefore, we proceeded to use only
multiplicative polynomials, both to avoid any non-physical flux values and for
the efficiency of the routine since combining both additive and multiplicative
Legendre polynomials is computationally expensive.

The stellar LOSVD moments can be seen in Fig. 3. For the systemic velocity of
NGC 1097, we used $1271$ km s-1, taken from the NASA/IPAC Extragalactic
Database (NED444NED can be accessed via https://ned.ipac.caltech.edu/) and
based on the redshift measurement in Allison et al., 2014 ($z$ $\approx$
0.00424). Further in this work, all wavelengths are corrected for the redshift
of the galaxy, equivalent to the rest frame wavelengths of the lines.

Figure 3: Stellar LOSVD moment maps of NGC 1097. Stellar velocity (top left)
velocity dispersion (top right) Gauss-Hermite moments h3 (bottom left) and h4
(bottom right)

To define the centre of the galaxy, we first attempted to fit ellipses to the
flux distribution by using the ellipse routine of IRAF555IRAF: the Image
Reduction and Analysis Facility software (Tody, 2000) software, but we
realised that the flux distribution of the galaxy is asymmetric, so the
parameters of the ellipses were poorly constrained. Hence, we defined the
photometric centre as the brightest spaxel in the nucleus based on the flux
integrated over the entire MUSE spectral range. The kinematic centre derived
in Sect. 3.3 differs only by $\sim$3 spaxels from thus defined photometric
centre. Hereafter, we give the coordinates $\Delta\alpha$, $\Delta\delta$ in
relation to the photometric centre, i.e. the photometric centre is located at
[$\Delta\alpha$,$\Delta\delta$]=[0,0].

####  2.2.2 Deriving ionised gas line-of-sight velocity distribution moments

Table 1: Characteristics of the fitted emission lines. Wavelengths of the emission lines are taken from the National Institute of Standards and Technology, (NIST; https://www.nist.gov/pml/atomic-spectra-database). Line ID |  | Wavelength
---
(Å)
| Amplitude Ratio
---
(fixed)
Group 1 - Hydrogen Balmer Lines
H$\beta$ |  4861.35 |
H$\alpha$ |  6562.79 |
Group 2 - Low Ionisation Lines
[N I]$\lambda$5198 |  5197.90 |
[N I]$\lambda$5200 |  5200.26 |
[N II]$\lambda$5754 |  5754.49 |
[O I]$\lambda$6300 |  6300.30 |
[O I]$\lambda$6364 |  6363.78 |  0.33[O I]$\lambda$6300
[N II]$\lambda$6548 |  6548.05 |  0.34[N II]$\lambda$6583
[N II]$\lambda$6583 |  6583.45 |
[S II]$\lambda$6716 |  6716.44 |
[S II]$\lambda$6731 |  6730.82 |
Group 3 - High Ionisation Lines
[O III]$\lambda$4959 |  4958.91 |  0.35[O III]$\lambda$5007
[O III]$\lambda$5007 |  5006.84 |
He I $\lambda$5876 |  5875.61 |
[S III]$\lambda$6312 |  6312.06 |

Figure 4: Flux and LOSVD moment maps of H$\alpha$ (left) [N II]$\lambda$6583Å
(middle), [O III]$\lambda$5007Å (right). We show the flux on logarithmic scale
in units $10^{-20}$ergs-1cm-2spaxel-1 (top row) velocity (middle row) and the
velocity dispersion (bottom row).

We derived the ionised gas kinematics by fitting initially single Gaussian
components to the emission lines, with the stellar continuum and absorption
lines fitted simultaneously. The fits were performed on individual spaxels,
while the stellar velocity moments were fixed to the values determined for
each Voronoi bin, specified in the prior stellar kinematics fits. Computing
Gauss-Hermite higher-order moments is also possible for the emission line
analysis, and we show the results of the fit including the h3 and h4
coefficients in Appendix A (Fig. A.3). However, we find fitting emission lines
with Gaussians to be a better approach when searching for multiple components
(see Sect. 2.3), whose presence otherwise could be incorporated in the higher-
order moments.

DAP can fit each emission line separately or fit several emission lines
together by assigning each the same velocity and velocity dispersion and by
imposing flux ratios. We left the fluxes free, except for the constraints from
atomic physics (Osterbrock & Ferland, 2006). We gathered the emission lines
into three main groups, with kinematics shared among lines in each group.
Extracting kinematics for each group separately allows us to study and
understand how gas regions of different physical properties move with respect
to each other. The grouping of the emission lines is done as follows: Group 1
consists of Balmer lines —H$\beta$$\lambda$4861, H$\alpha$$\lambda$6562
(indexed as H$\beta$ and H$\alpha$ in figures); Group 2 consists of low-
ionisation lines —[N I]$\lambda\lambda$5197,5200, [N II]$\lambda$5754, [O
I]$\lambda\lambda$6300,64, [N II]$\lambda\lambda$6548,83 , [S
II]$\lambda\lambda$6717,31 and Group 3 consists of high-ionisation lines —He
I$\lambda$5875, [O III]$\lambda\lambda$4959,5007, [S III]$\lambda$6312\. We
summarise the characteristics of the fitted emission lines in each group in
Table 1.

For each group, in Figure 4, we present their kinematics accompanied by the
flux maps of the strongest lines in each group, which are H$\alpha$, [N
II]$\lambda$6583 and [O III]$\lambda$5007 respectively. To avoid contribution
from spaxels dominated by noise, we used a threshold SNR$=$10; spaxels of
lower SNR are discarded from the maps. All three groups show high flux in the
nuclear ring and in the nucleus. The observed velocity field of each group has
a typical spider velocity pattern and is dominated by rotation, but it also
shows strong distortions from circular motion caused by the presence of the
bar. Such distortions in the gas kinematics have been reported in various
nearby galaxy studies (Bosma, 1978; Emsellem et al., 2001; Zurita et al.,
2004; Fathi, 2004; Fathi et al., 2006). The velocity dispersion is the lowest
in the nuclear ring and in the outer regions and significantly increases close
to the nucleus.

###  2.3 Improving the quality of the fits by fitting multiple Gaussians

Figure 5: Representative spectra zoomed on the region of H$\alpha$ and [N
II]$\lambda$6583, and the best fits from pPXF/DAP. The lines represent the
MUSE spectrum (black), the combined best fit of the stellar template and
emission lines (red), the best fit of emission lines (green), and the
difference between the spectrum and the combined best-fit stellar template
with emission lines (orange). top: Spectrum of a spaxel without the excess,
marked in the upper left panel of Fig. 6 with a ’+’. middle: Spectrum of a
spaxel with the excess affecting the single Gaussian fits of H$\alpha$ and [N
II]$\lambda$6583, marked in the upper left panel of Fig. 6 with an ’x’.
bottom: Spectrum of the spaxel shown in the middle panel, with fits that
include the secondary component. Figure 6: Comparison between fits using a
single Gaussian and fits including a secondary [N II]$\lambda$6583 component.
Grey contours in all panels highlight the H$\alpha$ emission. The left and
right columns represent H$\alpha$ and [N II]$\lambda$6583 respectively. top:
the corrected velocity dispersion. ’x’ and ’+’ markers on the left panel
represent the spaxels whose spectra are displayed in Fig. 5. middle: The ratio
of velocity dispersion obtained from single Gaussian fits to the corrected
velocity dispersion. bottom: The velocity differences between the corrected
velocity field and the velocity field obtained with single Gaussians.  Figure
7: left: The mask of spaxels for which the fits including secondary component
were used in further analysis. Spaxels in the grey-shaded region fall outside
the defined conditions and are excluded from the mask. Blue and red colours,
cut respectively at the 5 km s-1 and 4 km s-1 thresholds, show the difference
between velocities extracted by fitting single Gaussians and after including
the secondary component for H$\alpha$ and [N II]$\lambda$6583\. Grey contours
represent the H$\alpha$ emission. The centre of the galaxy is marked with a
cross. right: Flux and kinematics of the secondary component within the mask.
Only spaxels within the mask are shown. Maps respectively show the derived
flux (top left), velocity (top middle), velocity dispersion (top right), the
ratio of the fluxes of the secondary component and the main [N
II]$\lambda$6583 line (bottom left), the difference in velocities of the
secondary component and the main [N II]$\lambda$6583 line (bottom middle), the
ratio of velocity dispersion of the secondary component and the main [N
II]$\lambda$6583 line (bottom right).

In the regions inner to the nuclear ring, the kinematics of all emission line
groups reveal a region of high velocity dispersion, extending from the centre
of the galaxy $\sim$3″ towards the northeast (see Fig. 4), which drove us to
closely examine the spectra in the affected spaxels. In Figure 5, we show the
H$\alpha$ and [N II]$\lambda\lambda$6548,83 region of the spectrum in one of
the affected spaxels, as well as in a spaxel outside of the affected region.
These spaxels are located at the positions marked by an ’x’ and a ’+’,
respectively, in the top left panel of Fig. 6. In the spectrum of the affected
spaxel, as shown in the middle panel of Fig. 5, we notice excess emission in
the blue wing of [N II]$\lambda$6583 line, with a peak at the rest frame
wavelength $\sim$6574.50Å, which hereafter we call the excess. Because of the
excess, the line profile of [N II]$\lambda$6583 becomes strongly non-Gaussian,
and the excess encroaches on the region between [N II]$\lambda$6583 and
H$\alpha$, so the H$\alpha$ fits are also strongly affected by its presence.
As we demonstrate in the top panel of Fig. 5, when the excess is not present
the emission lines are well-fitted, but when the excess is present the
minimisation algorithm of pPXF produces a broad Gaussian, especially for
H$\alpha$, leading to overestimated velocity dispersion values. The presence
of the excess shifts the centre of the fitted Gaussians for the main [N
II]$\lambda$6583 component to slightly lower velocities, and for the H$\alpha$
component to higher velocities. We found no signs of an excess in the red wing
of [N II]$\lambda$6583, or either of the wings of H$\alpha$. Although we see a
weak excess in the blue wing of [N II]$\lambda$6548, its SNR is at the same
level as noise, and it is insignificant compared to the excess in the blue
wing of [N II]$\lambda$6583\. Moreover, we found no signs of an excess in the
wings of the [O III]$\lambda$5007, though its presence might be lost to the
noise. However, in the spaxels of the region with high velocity dispersion,
where the presence of the excess in the blue wing of [N II]$\lambda$6583 is
the most pronounced, [O III]$\lambda$5007 line is well-fitted. This suggests
that the high velocity dispersion of [O III]$\lambda$5007 in those regions
(Fig. 4) is intrinsic.

In order to account for the excess and to improve the line fits, we introduced
one additional emission component to the DAP configuration at wavelength
$\lambda\sim$6574.50Å, which we call the secondary component. Throughout the
rest of this work, we work under the assumption that the secondary component
is a Doppler-shifted component of [N II]$\lambda$6583\. We tested how the
derived kinematics depends on the initial guess of the wavelength of the
secondary component by shifting it by $\pm$3.3Å (roughly equal to the
dispersion broadening); we find no significant changes to the results. In the
bottom panel of Fig. 5, we demonstrate that for the spaxel, whose spectrum and
single Gaussian fits were affected by the presence of the excess, the addition
of the secondary component accounts for the excess and fixes the erroneous
fits. However, one must be careful, as for the spaxels where the excess is not
strong the fits of the secondary component would be poorly constrained, and we
must carefully determine the spaxels truly affected by the presence of the
excess. This can be done by comparing the measurements extracted from fits
with and without the secondary component. As we explained earlier in this
section, the excess emission on one side of the line not only increases the
velocity dispersion of a single Gaussian fit but also alters its velocity.
While high velocity dispersion can be intrinsic, a large change in the derived
velocity caused by including the secondary component indicates strong
asymmetry of the line profile, hence, likely the presence of multiple line
components. Therefore, we search for the spaxels that are most affected by the
presence of the excess by looking for a change in derived velocities, which
should be more reliable than using the high velocity dispersion measurements.

In order to identify which spaxels in and around the region of the elevated
velocity dispersion in the nucleus are affected by the presence of the excess,
we fitted the secondary component in the inner 8″$\times$8″. In Figure 6, we
show the change in velocity (V) and velocity dispersion ($\sigma$) between the
original values obtained from single Gaussian fits and the corrected values
obtained from fits with the secondary component included. For the corrected
measurements, in addition to the respective emission line ID, we use the
notation ”corr”.

The velocity dispersion in the nuclear regions remains high after including
the secondary component (upper panels), with the values for [N
II]$\lambda$6583 and H$\alpha$ velocity dispersion reduced by $\sim$20–30$\%$
and $\sim$30–40$\%$, respectively, from the values obtained from single
Gaussian fits which might indicate the velocity dispersion is intrinsically
high. The persisting high velocity dispersion might also suggest that our
approach is simplistic and that using more than two Gaussian components
possibly can account for the excess better and can reduce the estimated
velocity dispersion more. However, increasing the number of components
introduces more free parameters, and we do not have enough constraints or
resolution to fit more than two components. Moreover, we notice that the
largest reduction in velocity dispersion (middle panels) occurs in spaxels
where adding the secondary component significantly changes the velocity of the
main lines (lower panels).

Therefore, based on the bottom panels of Fig. 6 we defined the following
conditions to create a mask region where single Gaussian fits are not
reliable, and values from the run with the secondary component will be used:

  1. i.

$V_{\rm[N\ II]\lambda 6583}-V_{\rm[N\ II]\lambda 6583,corr}<-4$ km s-1

  2. ii.

$V_{\rm H\alpha}-V_{\rm H\alpha,corr}>5$ km s-1

  3. iii.

$-6\arcsec\leq\Delta\alpha\leq 1.2\arcsec$ and $-3\arcsec\leq\Delta\delta\leq
4\arcsec$

The application of the conditions above is shown in the left panel of Fig. 7.
The spatial cuts are applied to discard the isolated spaxels outside the
region of the elevated velocity dispersion in the nucleus where the velocity
correction is low. The velocity limits are chosen to maximise the area where
the secondary component can be reliably fitted. By replacing within the
defined mask region, the measurements from the single Gaussian fits with those
from fits with the secondary component included, we created a combined data
set of corrected flux, velocity and velocity dispersion. The resultant maps
are presented in Fig. A.2. For both H$\alpha$ and [N II]$\lambda$6583, the
apparent change in flux, velocity and velocity dispersion maps are not
significant. However, the accurate corrected values are crucial for our
further analysis in Sect. 3. Hereafter, we use measurements from the combined
data set. [O III]$\lambda$5007 can still be represented by the maps in Fig. 4
since single Gaussian fits were not affected by the presence of the excess.

We present the flux (F) and kinematic maps of the secondary component (indexed
by [N II]w) within the mask region in the right panel of Fig. 7, both absolute
and relative to the main [N II]$\lambda$6583 component. Note that few spaxels
within the mask show inconsistent values compared to the main gradient of flux
and kinematic measurements. This implies that the secondary component fits in
those spaxels are not reliable and those spaxels are discarded in the analysis
below. The secondary component has high flux at the centre of the galaxy that
extends towards the northeast side of the mask, coinciding with the high
velocity dispersion of the H$\alpha$ and [N II]$\lambda$6583 components. In
addition, the flux ratio of the secondary component to the main component
increases away from the centre. In the velocity map of the secondary
component, we see a gradient consistent with rotation, but after we subtract
the disk rotation, the secondary component appears blueshifted by 300–400 km
s-1 with respect to the main component, and the blueshift decreases away from
the galaxy centre. The velocity dispersions of the secondary and main
components are similar, though there is a slight indication that in spaxels
where the contribution from the secondary component is the strongest, its
velocity dispersion is smaller than that of the main component by up to 20%.

##  3 Analysis of the derived kinematics and fluxes

###  3.1 Searching for coherent kinematic structures in deviations from
circular motion

Coherent structures in velocity maps can be best traced when the rotational
component of motion is removed from the data. It is most commonly done by
subtracting the line-of-sight (LOS) velocity of a simple rotating disk from
the observed LOS velocity. The result is displayed in a spatial map, which is
called the residual velocity map. Coherent structures in the residual velocity
map will indicate large-scale motions of gas that depart from circular
trajectories. In particular, shocks in gas result in velocity discontinuities,
and large-scale shocks are expected to appear as coherent structures in the
residual velocity map.

In Figure 8, we demonstrate what we define as coherent structures in a
residual velocity map. In its upper panel, we present the LOS velocity of a
representative thin flat rotating disk which we constructed using Kinemetry
software (Krajnović et al., 2006), described in detail later in this section.
The disk orientation parameters, position angle (PA) of the line of nodes
(LON) of 137${\degr}$ and inclination (i) of 35${\degr}$, are taken from the
literature (Fathi et al., 2006). In the bottom panel of Fig. 8, we present the
constructed residual velocity map after subtracting the LOS velocity of the
modelled disk from the observed H$\alpha$ velocity. The residual velocity map
reveals two straight continuous structures, one with positive and one with
negative residual velocities, which are roughly parallel to the LON and are
offset from the LON to the top left and bottom right. Such large residual
velocities appearing over continuous stretches in the velocity map are
associated with bar-induced shocks in gas which cause velocity jumps. The
shocked gas loses angular momentum and funnels from the outer regions of the
galaxy towards the inner regions. As shocks compress the gas, which is well
mixed with dust, they also coincide with the loci of the dust lanes. Therefore
the alignment of the kinematic structures with large residual velocities shown
in Fig. 8 and the dust lanes shown in Fig. 1 signifies the presence of large-
scale shocks that extend from the outer regions of the bar towards the central
regions. Although gas inflow in bar-induced shocks is efficient, in Fig. 8 it
can be seen that the straight shocks weaken in the central regions. As the
result, inflow slows down, and the density of the gas increases which can
trigger SF and form stellar substructures like the nuclear rings (Shlosman et
al., 1989; Athanassoula, 1992; Martini et al., 2003; Kim et al., 2012). This
effect is seen in the central $\sim$8″ of NGC 1097, where the stretches of
large residual velocities are connecting to the nuclear ring (shown in Fig.
8).

Thus, although we have a good understanding of the bar driven gas inflows to
the central regions, we lack an understanding of the inflow to the innermost
nuclear regions. In this work, we aim to identify the inflow of gas from the
nuclear ring towards the nucleus of the galaxy by searching for coherent
structures in the residual velocity fields. Since subtracting the LOS velocity
of an incorrectly fitted disk from the observed velocity will obscure the real
features (van der Kruit & Allen, 1978), we need to construct the most accurate
disk model. To do this, we must cautiously determine the galaxy centre and
disk orientation parameters: PA of the LON and i. Below, we explain in detail
the procedures of constructing a flat disk in circular motion.

Figure 8: top: LOS velocity of the rotating disk fitted to observed H$\alpha$
velocity field with PA of the LON $=137{\degr}$ and $i=35{\degr}$. bottom:
Residual velocities after subtracting the LOS velocity of the modelled disk
from the observed H$\alpha$ velocity. Grey contours highlight the H$\alpha$
emission. The nuclear ring is represented by the black dashed circle. The
centre of the galaxy is marked with a cross. Figure 9: Radial distributions of
parameters fitted to the H$\alpha$ velocity field with Kinemetry. The blue-
shaded region represents the nuclear ring and the yellow-shaded region marks
the radial range where PA of the LON and i vary significantly (the Dip).
Results from the run with all parameters in the Kinemetry fit allowed to vary
freely are shown in pink, results for the centre and $A_{0}$ fixed are shown
in blue, and results for the PA of the LON, i, centre and $A_{0}$ fixed are
shown in purple. The left panels display the flattening, q (top) and the PA of
the LON (bottom). The right panels show the coordinates of the kinematic
centre (top), the $A_{0}$ coefficient (middle) and the rotation curve
(bottom).

For a flat thin disk in circular motion, the distribution of LOS velocities
can be expressed by:

$V(R,\psi)=V_{sys}+V_{C}(R)\ sin(i)\ cos(\psi)$ (1)

where $V_{sys}$ and $V_{C}$ are the systemic and circular velocity
respectively, i and $\psi$ are the inclination of the disk and azimuthal angle
measured counterclockwise from the LON in the galaxy plane, and $R$ is the
radius of a circular ring in the galaxy plane. Inclination is directly related
to the flattening parameter (or axial ratio) q ($cos(i)=q$) of the disk
projected on the sky, which is an ellipse with major axis on the LON and
ellipticity $\epsilon=1-q$.

To determine the orientation of the rotating disk model, we employed
Kinemetry, which is software specifically designed for IFU data and is a
generalisation of surface photometry to higher-order moments of the LOSVD
(Krajnović et al., 2006). Kinemetry performs harmonic expansion along best-
fitting ellipses of the optimal major axis’ PA, optimal ellipticity and centre
using Fourier analysis. When applied to the velocity field, it produces the
distribution of velocity along ellipses of semi-major axis a, K(a,$\psi$),
described by a finite number of harmonic terms:

$K(a,\psi)=A_{0}(a)+\sum_{n=1}^{N}A_{n}(a)sin(n\psi)+B_{n}(a)cos(n\psi)$ (2)

where $\psi$ is the eccentric anomaly, equivalent to that in Equation 1. To
determine the best-fitting ellipse, the appropriate harmonic terms are
minimised. Under the assumption of emission from a thin disk structure and
fixed centre, Kinemetry reduces to the tilted-rings method, which breaks the
disk into individual rings with different i, PA of the LON and rotational
velocity (Rogstad et al., 1974).

Figure 10:  Residual velocity map of the inner 30″$\times$30″. The dashed
circle represents the nuclear ring. The centre of the galaxy is marked with a
cross. In each panel, we subtract the LOS velocity of the disk modelled with a
different combination of PA of the LON and i. top: PA=137${\degr}$ and
i=35${\degr}$ ($q=0.82$) as in Fig. 8. The grey contours represent the
H$\alpha$ velocity field. middle: i is fixed at $35{\degr}$ and PA of the LON
is varied. bottom: PA of the LON is fixed at $137{\degr}$ and i is varied.

Since we want to fit the rotating disk to the gas motions, we use the
H$\alpha$ velocity field as the input for Kinemetry. To ensure the thin flat
rotating disk is accurately fitted we performed three Kinemetry fits, starting
from one with all orientation parameters allowed to vary and gradually fixing
the global orientation parameters in the consecutive fits. In Figure 9, we
summarise the estimated parameter distributions from each fit. In determining
the global values for the disk orientation parameters, we discard parameters
estimated within the inner 2″ as the estimations within this radius are
unreliable because of the insufficient data that ellipses are being fitted to,
and because of the impact of the AGN.

We first fitted the velocity field by allowing the orientation parameters of
each ellipse, including the centre, to vary freely. As shown in the top right
panel of Fig. 9, in the regions inward from the nuclear ring, the horizontal
($\Delta\alpha$) and vertical coordinates ($\Delta\delta$) of the kinematic
centre show only small variations ($\sim$0.7″ for $\Delta\alpha$ and
$\sim$0.5″ for $\Delta\delta$) from the photometric centre. This indicates
that in the regions inner to the nuclear ring, the estimated kinematic centre
agrees with the photometric centre within $\sim$3 spaxels. Also, as shown in
the middle right panel of Fig. 9, the $A_{0}$ values for ellipses fitted
interior to the nuclear ring depart from the systemic velocity of the galaxy
by less than 25 km s-1. Therefore, we reduced the number of free parameters
within Kinemetry by fixing the centre of the galaxy to the coordinates of the
photometric centre, and the $A_{0}$ to the velocity at the photometric centre
($-$3.55 km s-1), and performed a second fit. In this second fit, we intended
to estimate a representative PA of the LON and i which can accurately describe
the ellipses in nuclear regions, and we aimed to build our final circular disk
model by fixing the PA of the LON and i to these estimated representative
values. However, as shown in the left panels of Fig. 9, regardless of keeping
the centre free or fixed, in the regions inward from the nuclear ring, the
distributions of the PA of the LON and i strongly fluctuate within a radial
region marked in yellow, which we refer to as the Dip. Within the Dip, the PA
of the LON from the second fit varies between 125°–140°, and q varies between
0.65–0.90, which translates to an i range of $\sim$25°– 48°. While the impact
of bar-induced noncircular motions can justify the fluctuation of these
parameter values outside the nuclear ring, inward from the nuclear ring one
would expect motion close to circular and therefore well fitted by single
values for the PA of the LON and i. However, this expectation might be
simplistic, as NGC 1097 is proposed to be hosting a nuclear bar (Shaw et al.,
1993; Buta & Crocker, 1993) which would impact the gas flow in the regions
inward from the nuclear ring.

As we indicate above, in the third Kinemetry fit we want to fix the PA of the
LON and i to fit a flat disk, but it is challenging to deduce reliable
parameter values from such fluctuating distributions. Moreover, in Figure 8,
which presents residual velocities for a thin disk fit with parameters taken
from the literature (PA of the LON $=137{\degr}$ and $i=35{\degr}$, Fathi et
al., 2006), in the region inner to the nuclear ring we notice a periodic
azimuthal variation resembling an $m$ = 3 mode. Our initial explanation for it
was that i used in the disk model was higher than the true value, which would
give artefacts in the residuals in the form of $m$ = 3 harmonic terms (van der
Kruit & Allen, 1978; Franx et al., 1994; Schoenmakers et al., 1997). Artefacts
can be removed by choosing an i value in the disk model that is closer to the
actual i of the galaxy, but variations in Kinemetry parameters are preventing
us from settling on a reliable i.

Therefore, to find the i value for which the residuals of the Kinemetry fits
are free from the $m$ = 3 azimuthal variations, we investigated the effect of
changing the i and PA of the LON on residual velocities by running Kinemetry
first with a fixed PA of the LON at 137${\degr}$ and different q (0.7, 0.9),
and then with a fixed q at 0.82 ($i\approx$35°) and different PA of the LON
(125${\degr}$, 140${\degr}$). We constructed the corresponding residual
velocity maps and examined the residuals, as we show in Fig. 10, to deduce
which parameter set is best for the final disk model. When the disk is
modelled with the PA of the LON $125{\degr}$ and i $35{\degr}$, we see the $m$
= 1 harmonic terms in the residual velocity map as an effect of wrongly fitted
PA of the LON, with one side of the map having mainly positive velocities and
the other side having mainly negative velocities (van der Kruit & Allen,
1978). However, these artefacts can be largely eliminated by adopting the
optimal PA of the LON. We also see $m$ = 3 harmonic terms in all panels. They
are strongest in the bottom right panel where we notice both positive and
negative residual velocity regions extending radially from the centre beyond
the nuclear ring. We must emphasise that no matter what inclination is used in
the disk model, $m$ = 3 azimuthal variation in the residual velocity is only
minimised and not eliminated. The i within the range of 25${\degr}$–
35${\degr}$ minimises its amplitude, yet we do not see a clear minimum.

The range of acceptable parameter values agrees well with the parameters in
the literature: the kinematic position angle, PA=135${\degr}$, derived in
Storchi-Bergmann et al. (1996) and Emsellem et al. (2001); inclination,
i=34${\degr}$, derived in Storchi-Bergmann et al. (2003) for the nuclear
accretion disk, and i=35${\degr}$ derived in Fathi et al. (2006) from
exponential thin disk model. Therefore we conclude that the disk parameters,
PA=137${\degr}$ and $i$=35$\degr$, taken from Fathi et al. (2006) and applied
in Fig. 8, are optimal. Although the parameters, PA=122${\degr}$ and
$i$=48${{\degr}}$ adopted by Lang et al. (2020) differ from our optimal
parameters, we show in Fig. 10 that each of these parameters would give
artefacts in the residual velocity.

Due to the smoother nature of the stellar velocity field compared to the gas
velocity field, it is generally expected to obtain more continuous parameter
distributions when fitting a disk to the stellar velocity field. However, the
large velocity dispersion inward from the nuclear ring of NGC 1097 (Fig. 3)
reduces the LOS velocity values over an elliptical region inconsistent with a
projected disk in an acceptable inclination range, and therefore is likely to
result in an incorrect fit, which we indeed observed when using Kinemetry.
Therefore, parameters estimated by fitting the stellar velocity field would
not be reliable inward from the nuclear ring. On the other hand, the PA of the
LON and i that we derive with Kinemetry for the stellar velocity field outside
the nuclear ring is consistent with our adopted values, although the PA of the
LON diverges at larger radii, likely because of the influence of the bar.
Furthermore, disk orientation parameters cannot be reliably obtained from
large-scale photometry due to the presence of companion and asymmetric spiral
arms.

Since, even for the optimal parameters, a flat disk in circular motion does
not provide a good fit to the data, because periodic azimuthal variations of
residual velocities persist in those regions, we revisited the LOS velocity
fields of each emission line group to check if there are any signs of these
features directly in the observed velocity fields. In the top panel of Fig.
10, we superimpose H$\alpha$ velocity contours on the residual velocity map to
show that in the regions inward from the nuclear ring, the H$\alpha$ velocity
field shows wiggles at loci coinciding with two most prominent negative
residual velocity regions that we interpreted above as part of the $m$ = 3
harmonic terms. This indicates that the features which we initially
interpreted as artefacts of the fit arise from real features present in the
data. Hence, in Sect. 3.2, we focus on examining the physical processes in the
regions inner to the nuclear ring, which may give rise to some of the
kinematic structures observed in the velocity map.

Figure 11:  left: The BPT diagram of the inner 30″$\times$30″ region. The
SNR$=$10 and AON$=$7 thresholds are applied. Each point represents one spaxel.
Dashed, solid and dotted lines represent the Kauffmann et al. (2003), Kewley
et al. (2006) and Cid Fernandes et al. (2011) fits respectively which
establish the dominating ionisation mechanisms. The colour coding is based on
the distance from the Kewley et al. (2006) line. middle & right: Maps of the
inner 30″$\times$30″, with colour coding of the spaxels taken from the left
panel. Contours represent H$\alpha$ flux (middle) and LOS H$\alpha$ velocity
(right).

###  3.2 Dominating ionisation mechanisms in the centre of NGC 1097

To determine the dominating source of ionisation from optical emission lines,
such as AGN or star-forming regions, one can use the Baldwin, Phillips &
Terlevich (BPT) diagram (Baldwin et al., 1981). In the left panel of Fig. 11,
we show the BPT diagram for the inner 30″$\times$30″ region of NGC 1097. Each
point on the diagram corresponds to a spaxel in the data, with its position
determined by the fluxes measured for that spaxel. We mask spaxels with poor
flux measurements by applying an amplitude-over-noise (AON) ratio (Sarzi et
al., 2006) threshold of 7, which excludes spaxels with AON$\leq$7 from the BPT
diagram, choosing a smaller threshold would give noisy results already
occurring outside the nuclear ring in Fig. 11. This allows us to focus on the
emission that is coming from the nuclear ring and the region inner to it. To
distinguish AGN from star-forming regions, we used the Kewley et al., 2001
(hereafter Kewley+01) relation which indicates what line ratios can be
produced by stellar photoionization, as points above the Kewley+01 line likely
originate from the AGN activity or shocks. However, the Kewley+01 relation
does not take into account the composite emission which has a contribution
from SF, hot old low-mass stars and AGN. Therefore, to include the composite
emission, we used the modified Kewley+01 line from Kauffmann et al., 2003
(hereafter Kauffmann+03) so that the region between Kaufmann+03 and Kewley+01
lines is the location of the composite sources, and points below the
Kaufmann+03 line indicate emission purely coming from SF. We also separated
the AGN region into two classes, LINER and Seyfert region, based on the
relation given by Cid Fernandes et al. (2011).

Each point in the BPT diagram represents a spaxel and is colour coded
according to its distance from the Kewley+01 line. This colour coding is then
used in the spatial map to highlight which sources of ionisation dominate in
what regions (Fig. 11 central and right panels). We also included H$\alpha$
flux and velocity contours on the spatial maps to combine the BPT diagram
results with kinematics and light distribution.

In the nuclear ring, the emission is strongly dominated by SF, while the
regions inward from the nuclear ring predominantly show LINER emission.
However, three locations in the regions inward from the nuclear ring ($\sim$5″
east and west and $\sim$2″ northwest from the centre) show composite emission
due to contribution from SF and have high H$\alpha$ flux (see also Fig. 4).
The region of composite emission to the west has the highest contribution from
SF (is closest to the Kaufmann+03 line in the BPT diagram) and it appears to
be part of a structure with elevated SF that extends from the nuclear ring
directly north of it towards a smaller region of the composite emission to the
south. The H$\alpha$ velocity field throughout the extent of this structure is
strongly disturbed (Fig. 11 right panel), which can be caused by stellar
outflows associated with elevated SF. Moreover, this structure also coincides
with the strong negative velocity residual (Fig. 10), suggesting that the
negative velocity residuals, which we initially interpreted as artefacts
caused by adopting the wrong inclination of the fitted disk, arise from real
features present in the data. The region of composite emission to the east
shows a weaker SF component than the region to the west. This region is also
associated with a disturbed H$\alpha$ velocity field, and unlike the region to
the west which is extended tangentially in the north-south direction, it is
curved towards the nucleus of the galaxy. Lastly, the region of the composite
emission to the northwest has the weakest contribution from SF, and the
H$\alpha$ velocity field in that region is only mildly disturbed. Compared to
the other two regions, this region appears to be more isolated yet it still
clearly shows in the velocity residuals in Fig. 10.

###  3.3 Searching for coherent kinematic differences between distinct
emission line groups

Figure 12:  left: Line flux ratio of [N II]$\lambda$6583 and H$\alpha$.
middle: The velocity difference map, which is constructed by subtracting the
LOS velocity of H$\alpha$ from the LOS velocity of [N II]$\lambda$6583\. Black
contour marks the region where the velocity difference between [N
II]$\lambda$6583 and H$\alpha$ was more than 30 km s-1for single Gaussian
fits. right: Velocity dispersion ratio of [N II]$\lambda$6583 and H$\alpha$.
In all panels, grey contours highlight the H$\alpha$ emission.

We showed in Sect. 3.1 that velocity residuals after subtracting a flat disk
in circular motion do not provide conclusive information about coherent
kinematic structures inner to the nuclear ring. Here, we explore what
kinematic information can be extracted from the data without relying on
kinematic or mass models. In Fig. 12, we compare the LOSVD moments of the
Balmer and low ionisation lines by constructing the maps of [N
II]$\lambda$6583/H$\alpha$ line flux ratio (left panel), velocity difference
(middle panel) and velocity dispersion ratio (right panel).

[N II]$\lambda$6583 emission, with higher ionisation potential, may
predominantly originate from post-shock gas, while H$\alpha$ emission, mostly
generated by stellar photoionization, is less affected by the shock.
Therefore, in the regions where shocks are present, the observed [N
II]$\lambda$6583/H$\alpha$ line flux ratios should be elevated. As shown in
the left panel of Fig. 12, [N II]$\lambda$6583/H$\alpha$ is the lowest in the
nuclear ring where SF dominates the emission, and it is the highest in the
nucleus where the AGN dominates. There are no other regions inward from the
nuclear ring with elevated [N II]$\lambda$6583 emission. Instead, we notice
regions of elevated H$\alpha$ flux associated with enhanced SF (see Sect.
3.2), which might be obscuring the possible enhancement in [N II]$\lambda$6583
in the same region if SF is triggered by shocks. Nevertheless, if [N
II]$\lambda$6583 around shocks originates mostly from the post-shock gas,
large-scale shocks can manifest themselves by velocity differences between [N
II]$\lambda$6583 and H$\alpha$ that are coherent over large distances, as
explained below. The velocity of each tracer evaluated in a spaxel is the
first moment of the LOSVD arising from integrating the emission along the LOS
and over the spaxel area. If this includes regions of gas with different
physical conditions which move with respect to each other, the integrated
LOSVD recorded in a spaxel should vary between the tracers, as relative
contributions to line emission from regions with different physical conditions
differ between emission lines. Consequently, the value of its first moment may
differ for each tracer.

A good example of two neighbouring gas regions that differ in physical
conditions and move with respect to one another are regions upstream and
downstream from a shock in gas. Commonly temperature and pressure are higher
in the post-shock gas than in gas upstream from the shock, and jump conditions
in the shock impose a velocity difference between these two regions. Hence,
for a spaxel which records emission from pre-shock and post-shock regions, the
derived [N II]$\lambda$6583 and H$\alpha$ velocities will differ, as the
former will be more representative of the post-shock region. Then, if large-
scale shocks are present, they will appear as continuous structures in the
velocity difference map, in which the velocity difference takes consistently
positive (or negative) values. In the velocity difference map shown in the
middle panel of Fig. 12, we notice a coherent blue feature of negative
velocity difference spiralling out from the nucleus towards the southeast side
of the nuclear ring, hereafter called the blue spiral arm. If it is a
signature of shock in gas, it may be a spiral shock (Maciejewski, 2004)
extending inwards from the large-scale straight shocks in the bar of NGC 1097
(Fathi et al., 2006). The largest velocity difference in the blue spiral
occurs at the location where the blue spiral curves towards the nucleus and is
cospatial with the H$\alpha$ enhancement to the east and with enhanced SF (see
Fig. 11). This may indicate that the enhanced SF is associated with the inner
part of the spiral shock.

Furthermore, shock conditions give rise to random motions, which increase the
velocity dispersion of [N II]$\lambda$6583 in the post-shock regions. This can
possibly be seen in the southeast region of the blue spiral arm, where the
velocity dispersion of [N II]$\lambda$6583 is higher than H$\alpha$ but in the
region where the blue spiral curves towards the nucleus, the velocity
dispersion of [N II]$\lambda$6583 is only slightly above $\sim$1\. In the
region $\sim$5″west from the nucleus we notice the largest velocity dispersion
ratio, which is cospatial with the largest velocity differences to the west,
and the elevated [N II]$\lambda$6583/H$\alpha$ line flux ratios in the same
region, as well as the enhanced SF (Sect. 3.2), but this feature is less
spatially coherent than the blue spiral arm.

Although the observed velocity differences between emission lines could result
from obscuration, H$\alpha$ and [N II]$\lambda$6583 emission wavelengths
differ only by 20Å, so differential extinction is negligible. While incorrect
line fitting can also cause the difference in the derived velocities between
tracers, we minimise that possibility by addressing the non-Gaussianity of the
emission lines by including the secondary component of [N II]$\lambda$6583
line in spaxels where the excess emission was present, as we explain in Sect.
2.3. In fact, fitting emission lines with single Gaussian components resulted
in the largest velocity difference between H$\alpha$ and [N II]$\lambda$6583
in spaxels where the excess in the spectra was present. In the middle panel of
Figure 12, we mark with a black contour the region where the velocity
difference ranged from $-$30 to $-$50 km s-1for single Gaussian fits.
Including the secondary component reduced this difference to $\sim$10 km s-1.
The spectra throughout the blue spiral arm, which is outside the contour
region, are well-fitted with single Gaussians, which indicates that the
velocity difference between emission lines in those regions is not caused by
the incorrect fitting.

##  4 Discussion

###  4.1 Elliptical flow pattern in regions inward from the nuclear ring

Figure 13: left: H$\alpha$ (top) and stellar (bottom) LOS velocity maps. The
black line represents the corresponding zero velocity line fit. Grey contours
represent the H$\alpha$ emission. The centre of the galaxy is marked with a
cross. right: Distributions of the PA of the LON obtained from Kinemetry and
PA of the zero velocity line. For the stars the distribution of the PA of the
LON is plotted in small red circles and the PA of the zero velocity line in
red triangles. For gas, the PA of the LON is plotted in small blue circles and
the PA of the zero velocity line in blue triangles. The blue-shaded region
represents the nuclear ring, and the yellow-shaded region represents the Dip
defined in Sect. 3.1.

We based our studies of residual velocities (Sect. 3.1) on the assumption that
inward from the nuclear ring the flow of unperturbed gas is circular. However,
we showed that in this region Kinemetry cannot accurately constrain the
parameters of the unperturbed circular motion in a thin flat disk. The
possible reason for this is that the unperturbed flow is either not flat or
not circular there. In order to obtain meaningful residuals, one needs to
first correctly identify the underlying flow.

If the flow is intrinsically circular, the PA of the LON and PA of the zero
velocity line are perpendicular, but this is not the case if the global flow
is intrinsically not circular (e.g. elliptical; see the appendix of Mazzalay
et al., 2014). Here, we test whether in the regions inward from the nuclear
ring the underlying flow is intrinsically noncircular, by determining the
angle between the LON and the zero velocity line, where the distribution of
the PA of the LON is inferred from Kinemetry.

In order to derive the distribution of the PA of the zero velocity line, for
each column of spaxels in the velocity field (corresponding to constant Right
Ascension), we identified the spaxel with the velocity closest to the systemic
velocity of the galaxy. We then connected these spaxels to produce the fit of
the zero velocity line. Further, we determined the PA of each spaxel on the
zero velocity line which gives us the distribution of the PA of the zero
velocity line with radius. We applied these steps on both stellar and
H$\alpha$ velocity fields. We expect that the angle between the LON and the
zero velocity line will be less affected by non-circular motion for stars than
for gas because high velocity dispersion in stars hides the underlying non-
circularity.

In the left panels of Fig. 13, we show the H$\alpha$ and stellar LOS velocity
maps and the corresponding zero velocity line fitted on each map. In the right
panel, we present the individual distributions of the PA of the LON and the PA
of the zero velocity line for both velocity fields. In the figure, we show the
PA of the zero velocity line values after subtracting 90°, which should be
equal to the PA of the LON for circular motion. The distributions of the PA of
the LON and PA of the zero velocity line for stars roughly align with each
other over the whole radial range, differing slightly inside the Dip (see
Sect. 3.1) by at most $\sim$5°. The values of the PA of the LON and PA of the
zero velocity line for H$\alpha$ velocity field roughly converge in the
nuclear ring, but within the Dip, the two differ by $\sim$15°–20°. This
implies that in the regions inward from the nuclear ring, the underlying
orbits of the gas flow are intrinsically noncircular, which prevents Kinemetry
from accurately constraining the disk orientation.

Thus, assuming an intrinsically circular flow of gas may compromise the search
for coherent kinematic structures in the velocity residuals, and an oval or
elliptical flow can be a better representation of the gas flow. Such
assumption of an oval flow is allowed by the DiskFit software (Spekkens &
Sellwood, 2007), which is specifically developed for bar-like and oval
distortions and allows fitting non-axisymmetric flow patterns to two-
dimensional velocity fields. Although exploring this approach is plausible, it
is outside the scope of this work.

###  4.2 Evidence from the MUSE data for nuclear spiral shocks

Studying the near-infrared images from the adaptive optics assisted NACO
camera/spectrograph of the VLT, Prieto et al. (2005) found dust extinction in
the form of a complex three-arm filamentary structure in the regions inner to
the nuclear ring. The spiral arms of this structure are reaching to the centre
from the north, south and west. Moreover, in good agreement with Prieto et al.
(2005), a similar three-arm network is observed in the images from the
Advanced Camera for Surveys (ACS) of the HST (Fathi et al., 2006). Hereafter
we call this structure the nuclear dust filaments. In Figure 14, we show the
NACO image of Prieto et al. (2005). Inside the oval area, the image displays
the residuals after subtracting the light profile of a simple ellipse model.
In the figure, we include a dashed line to represent the proposed nuclear bar
of the galaxy (Shaw et al., 1993; Buta & Crocker, 1993), with a major axis PA
of 28°(Quillen et al., 1995). We notice that the northern and southern arms of
the nuclear dust filaments are roughly parallel to this nuclear bar, which is
consistent with the expected geometry of dust lanes in a bar when the dust
lanes are associated with gas inflow in shocks. The arm to the west of the
nucleus is not readily associated with the nuclear bar.

In Figure 14, in blue contours, we outline the blue spiral arm extending from
the inner southeast part of the nuclear ring towards the nucleus which is
revealed in the velocity difference map (Fig. 12, middle panel). The region of
the blue spiral where it curves towards the nucleus is cospatial with the
region of the composite emission directly to the east of the nucleus, which is
associated with enhanced SF, and with the disturbed H$\alpha$ velocity field
(Fig. 11). If the velocity difference is caused by the shock, as we argued in
Sect. 3.3, then the blue spiral arm can be a nuclear spiral shock in which gas
inflow to the nucleus occurs. For bisymmetric mass distribution, such shock
should have a symmetric red component on the other side of the nucleus
(Canzian, 1993; Schoenmakers et al., 1997) in the velocity difference map. In
Fig. 12, middle panel, we notice a red structure at $\sim$5″ west of the
nucleus extending from the nuclear ring towards the south, outlined in pink in
Fig. 14. This structure is not curved towards the nucleus and does not have as
defined spiral continuity as the blue spiral arm, but it has the highest
velocity differences, it is cospatial with the region of the composite
emission to the west, which has the highest contribution from SF, and it
coincides with the strongest distortions of the H$\alpha$ velocity field.
Hereafter, we call this structure the red spiral arm. We hypothesise that in
the red spiral arm the gas density is higher compared to the region to the
east, which leads to higher SF there, whose consequences, particularly
outflows, dominate the shock signature there, while the SF is less dominant in
the blue spiral.

Figure 14: NACO J-band image from Prieto et al. (2005) with the light profile
of a simple ellipse model subtracted inside the oval. The white dashed line
represents the proposed nuclear bar with PA=28${\degr}$ and length $\approx$
20″ (Quillen et al., 1995). Blue and pink contours represent the blue and red
spiral arms in the velocity difference map (Fig. 12). Grey contours represent
the H$\alpha$ emission. Figure 15: Residual velocities after subtracting the
LOS velocity of the disk model from the observed [N II]$\lambda$6583 velocity
field. The black dotted lines represent the kinematic spiral network of Fathi
et al. (2006). The black dashed circle represents the nuclear ring. The centre
of the galaxy is marked with a cross.

The blue and red spiral arms are shifted $\sim$2-3″ downstream from the
northern and southern arms of the nuclear dust filaments. Hence, the
combination of the structures observed in kinematics and morphology might be
the imprint of the gas flowing to the nucleus under forcing from the nuclear
bar. However, it is unclear whether the kinematic arms are directly associated
with the dust filaments.

Hydrodynamical and kinematic modelling (Davies et al., 2009; Fathi et al.,
2013) support a net molecular gas inflow to the nucleus, channelled via the
arms of the nuclear dust filaments of Prieto et al. (2005). Using data from
the Gemini Multi-Object Spectrograph (GMOS) on Gemini South Telescope, Fathi
et al. (2006) studied the [N II]$\lambda$6583 kinematics within the inner
15″$\times$7″ region of NGC 1097, and revealed a three-arm spiral network in
residual velocities centred on the galaxy nucleus with velocity amplitudes
reaching $\sim$20$\%$ of the rotational velocity, consistent with the nuclear
spiral models of Maciejewski (2004). Moreover, Davies et al. (2009) studied
the H2 kinematics within the inner 4″$\times$4″ region and could resolve a
two-arm kinematic spiral structure in their residual velocities. Benefiting
from the large FoV of MUSE, we can trace these observed kinematic structures
in our residual velocity map beyond their spatial coverage. We do not find any
correlation between the structures found in the residual velocities by Davies
et al. (2009) and our results, but the FoV of SINFONI is much smaller than the
region to which we fitted the disk – Note that the SINFONI data have higher
spatial resolution than the MUSE data, which might possibly make it harder to
find a match as well.

Since Fathi et al. (2006) used [N II]$\lambda$6583 data in their analysis, to
compare our findings, we fitted our disk model with fixed PA of the
LON=137${\degr}$ and i=35${\degr}$ (see Sect. 3.1) to [N II]$\lambda$6583
velocity field. The resulting residual velocity map is presented in Fig. 15,
which also includes the proposed three-arm kinematic spiral of Fathi et al.
(2006). We do not find a coherent spiral network extending from the galaxy
nucleus towards outer radii. Arguably, the proposed arm of the kinematic
spiral to the south appears connected with a coherent structure in our
residual velocity map, which might be extending from the nuclear ring towards
the nuclear regions, but instead of curving towards the nucleus, it continues
to the northeast. As already seen in Fathi et al. (2006) this arm runs along
the arm of the nuclear dust filaments to the west. Hence, although in the
region inward from the nuclear ring a three-arm spiral network is present in
morphology (Prieto et al., 2005; Fathi et al., 2006), we do not find evidence
of a spiral network in the residual velocities of the ionised gas. As we
explained in Sects. 3.1, 3.2 and 4.1, the data inward from the nuclear ring do
not allow an accurate fit of a flat disk in circular motion, therefore
residual velocities have to be interpreted with caution. On the other hand, in
this work we find indications of nuclear spiral shocks in velocity difference
maps.

###  4.3 Multi-phase AGN outflows in the nuclear regions

As confirmed by a number of studies, the gas funnelling into the nucleus of
NGC 1097 feeds its LINER core (Davies et al., 2009; van de Ven & Fathi, 2010;
Fathi et al., 2013; Izumi et al., 2013; Martín et al., 2015; Izumi et al.,
2017; Prieto et al., 2019). Some of those studies also attempt to compute the
mass inflow rate ($\dot{M}$) from the nuclear ring of the galaxy into its
central black hole. For example, van de Ven & Fathi (2010) estimated the
ionised $\dot{M}$ at $\sim$70pc from the galaxy nucleus as
$\sim$0.011$M_{\odot}yr^{-1}$, Fathi et al. (2013) estimated a net $\dot{M}$
as $\sim$0.2$M_{\odot}yr^{-1}$ at $\sim$40pc and $\sim$0.6$M_{\odot}yr^{-1}$
at $\sim$100pc from the galaxy nucleus by using combined molecular and ionised
gas observations, and Davies et al. (2009) estimated a net $\dot{M}$ of
0.06$M_{\odot}yr^{-1}$ to the central few tens of parsecs regions. Magnetic
fields in the region inner to the nuclear ring that can affect gas inflow
rates have been observed (Lopez-Rodriguez et al., 2021; Hu et al., 2022).

It is known that the inflow of gas into the central black hole can feed the
AGN and trigger powerful outflows (García-Burillo et al., 2005; Davies et al.,
2014). Even though this phenomenon is more commonly seen in Seyfert nuclei
(García-Burillo et al., 2014; Alonso-Herrero et al., 2019; Audibert et al.,
2020), outflows can still be powered in less luminous cores (Nemmen et al.,
2006; Santoro et al., 2020), as supported by recent multi-wavelength studies
of nearby LINER galaxies (Márquez et al., 2017; Hermosa Muñoz et al., 2022).
In this section, we discuss the imprint of ionised gas outflows in the nucleus
of NGC 1097, which may be triggered due to the inflow of gas into the central
black hole.

As we discussed in Sect. 2.3, we noticed a high velocity dispersion region
extending from the centre of the galaxy towards $\sim$3″ northeast (see Fig.
4). One reason for the high velocity dispersion values is the excess to the
blue wing of [N II]$\lambda$6583 line, which is affecting the line profiles
and the Gaussian fits (see Fig. 5). However, even after modelling the excess
as a separate secondary emission component, the velocity dispersion of the
main component remained high (see Fig. 6). Contribution from the secondary
component to the total [N II]$\lambda$6583 flux appears to increase with
radius and its velocity is high and blueshifted compared to the main [N
II]$\lambda$6583 line (Fig. 7).

X-ray/Chandra observations adopted by Prieto et al. (2019), highlight the
X-ray-rich nucleus of the galaxy along with the energetic extension blasting
from the nucleus towards the northeast, northwest and southwest. The northeast
X-ray extent coincides with the high velocity dispersion region. A similar
correlation between ionised H$\alpha$ emission and soft X-ray observations has
been reported by previous studies on nearby LINER galaxies (Márquez et al.,
2017; Hermosa Muñoz et al., 2022). Moreover, high velocity dispersion of
various molecular gas tracers is also observed in the regions inward from the
nuclear ring (Izumi et al., 2017; Leroy et al., 2021).

We hypothesise that the strong disturbances in the kinematics of the multiple
gas phases and ionised gas disturbances coinciding with the loci of X-ray
emission imply the presence of AGN outflows, triggered by the gas inflows into
the black hole. In particular, emission from the secondary component could be
coming from the jet-like outflow, whose interaction with the disk increases
the velocity dispersion in the latter.

##  5 Summary & Conclusions

In this pilot study, we presented methods that can be used in the search for
gas inflow in centres of galaxies, and we applied them to a detailed study of
the kinematics of ionised gas in the inner 60″$\times$60″ region of NGC 1097.
Using data from MUSE/VLT, we derived kinematic maps for emission lines
representing three distinct groups: Balmer, low ionisation and high ionisation
lines.

In the regions inward from the nuclear ring, we improved the quality of line
fitting by including a secondary component to the [N II]$\lambda$6583 emission
line in order to account for its asymmetric line profile. The morphology and
kinematics of gas giving rise to the secondary component, when compared with
previous multi-wavelength studies, indicate that this emission could be
related to the outflows from the AGN.

Shocks in gas, especially those of coherent span over large distances, like
straight shocks in bars and spiral shocks, have long been associated with
inflow, therefore we argue that gas inflow can be revealed by coherent
kinematic structures that large-scale shocks impress in kinematic maps. These
structures are best highlighted when the revolving motion around the centre is
removed, but we showed that subtracting the circular velocity of a flat disk
may not reveal coherent kinematic structures when gas is not moving on
circular orbits and its revolution around the centre is better approximated by
motion along an oval. In addition, local perturbations, like outflows from SF
regions can distort the appearance of coherent kinematic structures induced by
large-scale shocks in gas.

Instead of subtracting the circular velocity of a modelled disk, the
contribution from the revolving motion around the centre can be minimised in a
model-independent way by subtracting velocity maps of two different emission
lines, which should highlight differences in kinematics between two tracers.

The map of velocity difference between [N II]$\lambda$6583 and H$\alpha$
constructed for the central kiloparsec of NGC 1097 reveals a coherent
structure in the form of a spiral, extending from the inner southeast part of
the nuclear ring to the north and curving towards the nucleus (the ”blue
spiral arm”), whose southeastern region shows higher velocity dispersion of [N
II]$\lambda$6583 compared to H$\alpha$. We associate the coherent velocity
difference in this structure with shock conditions and argue for the blue
spiral arm being a large-scale shock in the gas. In the velocity difference
map, the blue spiral arm appears to have a symmetric companion on the other
side of the nucleus (the ”red spiral arm”). Although less coherent than the
blue spiral arm, it coincides with the region of high SF and strongly
disturbed H$\alpha$ velocity field, which indicates that if the red spiral arm
is also a shock in gas, it is strongly affected by SF activity. The
bisymmetric geometry of the blue and red spiral arms is consistent with
predictions from nuclear spiral models. The kinematic arms that we found in
NGC 1097 run along the southern and northern spiral arms of the nuclear dust
filaments, located $\sim$3″ downstream from them. Together, they appear
aligned well with the nuclear bar, which suggests that the nuclear bar can
play a role in gas inflow in the centre of NGC 1097.

The methodology developed in this pilot study provides new means of finding
coherent kinematic structures in gas associated with shocks and inflows; by
applying it to the diverse galaxy sample of nearby galaxies contained in CBS,
in the following paper we intend to determine how common gas inflows are in
coherent kinematic structures.

## Acknowledgements

TK acknowledges the joint studentship support from Liverpool John Moores
University the Faculty of Engineering and Technology, and the Science
Technology Facilities Council. TK acknowledges the studentship support from
European Southern Observatory. TK thanks Dr Almudena Prieto for sharing the
processed NACO J-band image of NGC 1097, Professor Philip A. James for useful
comments on the manuscript, Professor Eric Emsellem for suggestions and
discussions on data analysis, Dr Francesco Belfiore for valuable help on
general enquiries with the Data Analysis Pipeline, and Dr Adrian Bittner for
help on the GIST pipeline. TK especially thanks Alonso Luna Ruiz Fernández for
his continuous support and encouragement.
JN acknowledges funding from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme (grant
agreement No. 694343).

##  6 Data Availability

MUSE data underlying this work is accessed via the data archive of the
European Southern Observatory, under proposal ID:0097.B-0640 (PI: Dimitri
Gadotti) provided by the TIMER collaboration, at
http://archive.eso.org/wdb/wdb/adp/phase3$_$main/form.

## References

  * Allison et al. (2014) Allison J. R., Sadler E. M., Meekin A. M., 2014, MNRAS, 440, 696
  * Alonso-Herrero et al. (2019) Alonso-Herrero A., et al., 2019, A&A, 628, A65
  * Athanassoula (1992) Athanassoula E., 1992, MNRAS, 259, 345
  * Audibert et al. (2019) Audibert A., et al., 2019, A&A, 632, A33
  * Audibert et al. (2020) Audibert A., et al., 2020, arXiv e-prints, p. arXiv:2011.09133
  * Audibert et al. (2021) Audibert A., Combes F., García-Burillo S., Dasyra K., 2021, IAU Symposium, 359, 307
  * Bacon et al. (2010) Bacon R., et al., 2010, in McLean I. S., Ramsay S. K., Takami H., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III. p. 773508, doi:10.1117/12.856027
  * Baldwin et al. (1981) Baldwin J. A., Phillips M. M., Terlevich R., 1981, PASP, 93, 5
  * Bittner et al. (2019) Bittner A., et al., 2019, A&A, 628, A117
  * Bittner et al. (2020) Bittner A., et al., 2020, A&A, 643, A65
  * Bosma (1978) Bosma A., 1978, PhD thesis, University of Groningen, Netherlands
  * Bowen et al. (2016) Bowen D. V., Chelouche D., Jenkins E. B., Tripp T. M., Pettini M., York D. G., Frye B. L., 2016, ApJ, 826, 50
  * Buta & Crocker (1993) Buta R., Crocker D. A., 1993, AJ, 106, 939
  * Canzian (1993) Canzian B., 1993, ApJ, 414, 487
  * Cappellari (2017) Cappellari M., 2017, MNRAS, 466, 798
  * Cappellari & Copin (2003) Cappellari M., Copin Y., 2003, MNRAS, 342, 345
  * Cappellari & Emsellem (2004) Cappellari M., Emsellem E., 2004, PASP, 116, 138
  * Cid Fernandes et al. (2011) Cid Fernandes R., Stasińska G., Mateus A., Vale Asari N., 2011, MNRAS, 413, 1687
  * Davies et al. (2009) Davies R. I., Maciejewski W., Hicks E. K. S., Tacconi L. J., Genzel R., Engel H., 2009, ApJ, 702, 114
  * Davies et al. (2014) Davies R. I., et al., 2014, ApJ, 792, 101
  * Diehl & Statler (2006) Diehl S., Statler T. S., 2006, MNRAS, 368, 497
  * Donnari et al. (2021) Donnari M., et al., 2021, MNRAS, 500, 4004
  * ESO Press Release (2004) ESO Press Release 2004, The President and the Galaxy, ESO Press Release, 12/2004
  * Ellison et al. (2011) Ellison S. L., Nair P., Patton D. R., Scudder J. M., Mendel J. T., Simard L., 2011, MNRAS, 416, 2182
  * Emsellem et al. (2001) Emsellem E., Greusard D., Combes F., Friedli D., Leon S., Pécontal E., Wozniak H., 2001, A&A, 368, 52
  * Emsellem et al. (2003) Emsellem E., Goudfrooij P., Ferruit P., 2003, MNRAS, 345, 1297
  * Emsellem et al. (2022) Emsellem E., et al., 2022, A&A, 659, A191
  * Erwin et al. (2021) Erwin P., et al., 2021, MNRAS, 502, 2446
  * Erwin et al. (2023) Erwin P., et al., 2023, in prep
  * Fathi (2004) Fathi K., 2004, PhD thesis, University of Groningen
  * Fathi et al. (2005) Fathi K., van de Ven G., Peletier R. F., Emsellem E., Falcón-Barroso J., Cappellari M., de Zeeuw T., 2005, MNRAS, 364, 773
  * Fathi et al. (2006) Fathi K., Storchi-Bergmann T., Riffel R. A., Winge C., Axon D. J., Robinson A., Capetti A., Marconi A., 2006, ApJ, 641, L25
  * Fathi et al. (2013) Fathi K., et al., 2013, ApJ, 770, L27
  * Fragkoudi et al. (2016) Fragkoudi F., Athanassoula E., Bosma A., 2016, MNRAS, 462, L41
  * Franx et al. (1994) Franx M., van Gorkom J. H., de Zeeuw T., 1994, ApJ, 436, 642
  * Freedman et al. (2001) Freedman W. L., et al., 2001, The Astrophysical Journal, 553, 47
  * Gadotti et al. (2019) Gadotti D. A., et al., 2019, MNRAS, 482, 506
  * Gadotti et al. (2020) Gadotti D. A., et al., 2020, A&A, 643, A14
  * García-Burillo et al. (2005) García-Burillo S., Combes F., Schinnerer E., Boone F., Hunt L. K., 2005, A&A, 441, 1011
  * García-Burillo et al. (2014) García-Burillo S., et al., 2014, A&A, 567, A125
  * Heller & Shlosman (1994) Heller C. H., Shlosman I., 1994, ApJ, 424, 84
  * Hermosa Muñoz et al. (2022) Hermosa Muñoz L., Márquez I., Cazzoli S., Masegosa J., Agís-González B., 2022, A&A, 660, A133
  * Ho (2008) Ho L. C., 2008, ARA&A, 46, 475
  * Hu et al. (2022) Hu Y., Lazarian A., Beck R., Xu S., 2022, arXiv e-prints, p. arXiv:2206.05423
  * Izumi et al. (2013) Izumi T., et al., 2013, PASJ, 65, 100
  * Izumi et al. (2017) Izumi T., et al., 2017, ApJ, 845, L5
  * Kauffmann et al. (2003) Kauffmann G., et al., 2003, MNRAS, 346, 1055
  * Kewley et al. (2001) Kewley L. J., Dopita M. A., Sutherland R. S., Heisler C. A., Trevena J., 2001, ApJ, 556, 121
  * Kewley et al. (2006) Kewley L. J., Groves B., Kauffmann G., Heckman T., 2006, MNRAS, 372, 961
  * Kim et al. (2012) Kim W.-T., Seo W.-Y., Stone J. M., Yoon D., Teuben P. J., 2012, ApJ, 747, 60
  * Krajnović et al. (2006) Krajnović D., Cappellari M., de Zeeuw P. T., Copin Y., 2006, MNRAS, 366, 787
  * Lang et al. (2020) Lang P., et al., 2020, ApJ, 897, 122
  * Legodi et al. (2021) Legodi L. S., Taylor A. R., Stil J. M., 2021, MNRAS, 500, 576
  * Leroy et al. (2021) Leroy A. K., et al., 2021, The Astrophysical Journal Supplement Series, 257, 43
  * Lin et al. (2013) Lin L.-H., Wang H.-H., Hsieh P.-Y., Taam R. E., Yang C.-C., Yen D. C. C., 2013, ApJ, 771, 8
  * Lopez-Rodriguez et al. (2021) Lopez-Rodriguez E., et al., 2021, ApJ, 923, 150
  * Maciejewski (2004) Maciejewski W., 2004, MNRAS, 354, 892
  * Maciejewski et al. (2002) Maciejewski W., Teuben P. J., Sparke L. S., Stone J. M., 2002, MNRAS, 329, 502
  * Man et al. (2019) Man Z.-y., Peng Y.-j., Kong X., Guo K.-x., Zhang C.-p., Dou J., 2019, MNRAS, 488, 89
  * Márquez et al. (2017) Márquez I., Masegosa J., González-Martin O., Hernández-Garcia L., Pović M., Netzer H., Cazzoli S., del Olmo A., 2017, Frontiers in Astronomy and Space Sciences, 4, 34
  * Martín et al. (2015) Martín S., et al., 2015, A&A, 573, A116
  * Martini et al. (2003) Martini P., Regan M. W., Mulchaey J. S., Pogge R. W., 2003, ApJ, 589, 774
  * Mazzalay et al. (2014) Mazzalay X., et al., 2014, MNRAS, 438, 2036
  * Mazzuca et al. (2008) Mazzuca L. M., Knapen J. H., Veilleux S., Regan M. W., 2008, ApJS, 174, 337
  * Nemmen et al. (2006) Nemmen R. S., Storchi-Bergmann T., Yuan F., Eracleous M., Terashima Y., Wilson A. S., 2006, ApJ, 643, 652
  * Osterbrock & Ferland (2006) Osterbrock D. E., Ferland G. J., 2006, Astrophysics of gaseous nebulae and active galactic nuclei. University Science Books
  * Prieto et al. (2005) Prieto M. A., Maciejewski W., Reunanen J., 2005, AJ, 130, 1472
  * Prieto et al. (2019) Prieto M. A., Fernandez-Ontiveros J. A., Bruzual G., Burkert A., Schartmann M., Charlot S., 2019, MNRAS, 485, 3264
  * Quillen et al. (1995) Quillen A. C., Frogel J. A., Kuchinski L. E., Terndrup D. M., 1995, AJ, 110, 156
  * Rogstad et al. (1974) Rogstad D. H., Lockhart I. A., Wright M. C. H., 1974, ApJ, 193, 309
  * Sanders & Huntley (1976) Sanders R. H., Huntley J. M., 1976, ApJ, 209, 53
  * Santoro et al. (2020) Santoro F., Tadhunter C., Baron D., Morganti R., Holt J., 2020, A&A, 644, A54
  * Sarzi et al. (2006) Sarzi M., et al., 2006, MNRAS, 366, 1151
  * Schoenmakers et al. (1997) Schoenmakers R. H. M., Franx M., de Zeeuw P. T., 1997, MNRAS, 292, 349
  * Shaw et al. (1993) Shaw M. A., Combes F., Axon D. J., Wright G. S., 1993, A&A, 273, 31
  * Shlosman (2001) Shlosman I., 2001, in Knapen J. H., Beckman J. E., Shlosman I., Mahoney T. J., eds, Astronomical Society of the Pacific Conference Series Vol. 249, The Central Kiloparsec of Starbursts and AGN: The La Palma Connection. p. 55
  * Shlosman et al. (1989) Shlosman I., Frank J., Begelman M. C., 1989, Nature, 338, 45
  * Sormani et al. (2015) Sormani M. C., Binney J., Magorrian J., 2015, MNRAS, 449, 2421
  * Spekkens & Sellwood (2007) Spekkens K., Sellwood J. A., 2007, ApJ, 664, 204
  * Storchi-Bergmann et al. (1996) Storchi-Bergmann T., Rodriguez-Ardila A., Schmitt H. R., Wilson A. S., Baldwin J. A., 1996, ApJ, 472, 83
  * Storchi-Bergmann et al. (2003) Storchi-Bergmann T., et al., 2003, ApJ, 598, 956
  * Storchi-Bergmann et al. (2007) Storchi-Bergmann T., Dors Oli L. J., Riffel R. A., Fathi K., Axon D. J., Robinson A., Marconi A., Östlin G., 2007, ApJ, 670, 959
  * Tody (2000) Tody D., 2000, in Murdin P., ed., , Encyclopedia of Astronomy and Astrophysics. p. 2923, doi:10.1888/0333750888/2923
  * Vazdekis et al. (2015) Vazdekis A., et al., 2015, MNRAS, 449, 1177
  * Zurita et al. (2004) Zurita A., Relaño M., Beckman J. E., Knapen J. H., 2004, A&A, 413, 73
  * van de Ven & Fathi (2010) van de Ven G., Fathi K., 2010, ApJ, 723, 767
  * van der Kruit & Allen (1978) van der Kruit P. C., Allen R. J., 1978, ARA&A, 16, 103
  * van der Marel & Franx (1993) van der Marel R. P., Franx M., 1993, ApJ, 407, 525

##  Appendix A Supplementary Maps

Figure A.1: Flux maps of the weaker lines in each emission line group outlined
in Table 1. Units are in $10^{-20}$ergs-1cm-2spaxel-1 units.

Figure A.2: Corrected flux and LOSVD moments after combining the results from
single Gaussian fits with the fits including the secondary component (as
explained in Sect. 2.3). For H$\alpha$ line (left), and for [N
II]$\lambda$6583 line (right). We display the line flux on a logarithmic scale
in units $10^{-20}$ergs-1cm-2spaxel-1 (top row), velocity (middle row) and
velocity dispersion (bottom row). Figure A.3: Flux and LOSVD moments for fits
including higher order moments (h3 and h4). Columns show H$\alpha$ (left), [N
II]$\lambda$6583 (middle) and [O III]$\lambda$5007 (right). The units of the
flux measurements are $10^{-20}$ergs-1cm-2spaxel-1.