ALMA captures feeding and feedback from the active galactic nucleus in NGC613

We report ALMA observations of CO(3-2) emission in the Seyfert galaxy NGC613, at a spatial resolution of 17pc, as part of our NUclei of GAlaxies sample. Our aim is to investigate the morphology and dynamics of the gas inside the central kpc, and to probe nuclear fueling and feedback phenomena. The morphology of CO(3-2) line emission reveals a 2-arm trailing nuclear spiral at $r\lesssim$100pc and a circumnuclear ring at ~350pc radius, that is coincident with the star-forming ring seen in the optical images. The molecular gas in the galaxy disk is in a remarkably regular rotation, however, the kinematics in the nuclear region is very skewed. The nuclear spectrum of CO and dense gas tracers HCN(4-3), HCO+(4-3), and CS(7-6) show broad wings up to $\pm$300km/s, associated with a molecular outflow emanating from the nucleus (r~25pc). We derive a molecular outflow mass $M_{out}$=2x10$^6$M$_\odot$ and a mass outflow rate of $\dot{M}_{out}=$27$\rm M_\odot yr^{-1}$. The molecular outflow energetics exceed the values predicted by AGN feedback models: the kinetic power of the outflow corresponds to 20%L_AGN and the momentum rate is $\dot{M}_{out}v\sim400L_{AGN}/c$. The outflow is mainly boosted by the AGN through entrainment by the radio jet, but given the weak nuclear activity of NGC613, we might be witnessing a fossil outflow, resulted from a strong past AGN that now has already faded. The nuclear trailing spiral observed in CO emission is inside the ILR ring of the bar. We compute the gravitational torques exerted in the gas to estimate the efficiency of the angular momentum exchange. The gravity torques are negative from 25 to 100pc and the gas loses its angular momentum in a rotation period, providing evidence of a highly efficient inflow towards the center. This phenomenon shows that the massive central black hole has a significant dynamical influence on the gas, triggering its fueling.


Introduction
The energy of active galactic nuclei (AGN) is well understood as being due to gas accretion onto the supermassive black hole (SMBH; Antonucci 1993). Gas inflows into the center of galaxies can fuel the SMBH and the energy input by the AGN can trig-the timescales involved, since both feeding processes rely on a common cold gas supply, but in different periods of time (∼10 5 yr for BH growth and ∼10 7−9 yr for star formation, García-Burillo et al. 2016). These timescales are related to the mechanisms that drive the gas from galactic scales (∼10 kpc) to nuclear scales (a few pc), through removal of angular momentum; large nonaxisymmetric perturbations, such as bars or spirals, represent one of these mechanisms, but so far no unique physical process associated with inward transport of gas in galaxy disks has been identified.
On large scales, cosmological simulations show that mergers and galaxy interactions are able to produce strong nonaxisymmetries (e.g., Hopkins et al. 2006;Di Matteo et al. 2008). On kiloparsec scales, bar instabilities, either internally driven by secular evolution or triggered by companions, can first feed a central starburst and then fuel the BH (García-Burillo et al. 2005). On the other hand, gas inflow is impeded by the inner Lindblad resonance (ILR), where the gas is trapped in a nuclear ring (see Piner et al. 1995;Regan & Teuben 2004). On scales of a few hundred parsecs, the "bars within bars" scenario (e.g., Shlosman et al. 1989), together with m = 1 instabilities and nuclear warps (Schinnerer et al. 2000), takes over as a dynamical mechanism (see e.g., Hunt et al. 2008). Observations of nearby low-luminosity AGN (LLAGN) with the NUGA program have revealed smoking-gun evidence of AGN fueling in one third of galaxies (García-Burillo & Combes 2012). This result suggests that galaxies alternate periods of fueling and starvation, and might be found in a feeding phase at 300 pc scales only one third of the time.
As we approach the center of galaxies, other mechanisms can contribute to the fueling: viscous torques, from dense gas in regions of large shear, or dynamical friction can drive massive clouds to the nucleus (e.g., Combes 2003;Jogee et al. 2006). Simulations suggest that fueling involves a series of dynamical instabilities (m = 2, m = 1) on scales of ∼10 pc, and also predict the formation of a thick gas disk similar to the putative torus invoked to explain obscured AGN (Hopkins & Quataert 2010;Hopkins et al. 2012). These fueling episodes are eventually quenched by either nuclear star formation winds or AGN feedback.
The recent discovery of many massive (a few 10 7 M ) molecular outflows in nearby AGN (e.g., Fischer et al. 2010;Feruglio et al. 2010Feruglio et al. , 2017Alatalo et al. 2011;Sturm et al. 2011;Veilleux et al. 2013;Cicone et al. 2014;Sakamoto et al. 2014;García-Burillo et al. 2014;Dasyra & Combes 2012;Dasyra et al. 2014) has promoted the idea that winds may be major actors in sweeping gas out of galaxies. It has already been established that mass-outflow rates increase with AGN luminosity, supporting the idea of a luminous AGN pushing away the surrounding gas through a fast wind. Observational works (Cicone et al. 2014;Fiore et al. 2017;Fluetsch et al. 2019) have shown that molecular outflow properties are correlated with AGN luminosity, where the outflow kinetic power corresponds to about 5%L AGN and the momentum rate is ∼20L AGN /c, in agreement with theoretical models of AGN feedback (Faucher-Giguère & Quataert 2012;Zubovas & King 2012. Outflows have been traced for a long time in ionized or atomic gas (Rupke et al. 2005a;Riffel & Storchi-Bergmann 2011), making it now possible to compare the different gas phases. Carniani et al. (2015) found that ionized gas only traces a small fraction of the total gas mass, suggesting that the molecular phase dominates the outflow mass. This trend is also found by Fiore et al. (2017), but the ratio between molecular and ionized mass-outflow rates is reduced at the highest AGN bolometric luminosities.
Probing AGN feeding and feedback phenomena through the kinematics and morphology of the gas inside the central kiloparsec has only recently been possible due to the unprecedented spatial resolution and sensitivity of the Atacama Large Millimeter/submillimeter Array (ALMA). Evidence of AGN feeding was found in NGC 1566, where a molecular trailing spiral structure from 50 to 300 pc was detected with ALMA Cycle 0 observations, and according to its negative gravity torques, is found to be contributing to the fueling of the central BH . Again using Cycle 0 observations, a molecular outflow is also seen in the LLAGN in the Seyfert 2 NGC 1433. It is the least massive molecular outflow (∼4 × 10 6 M ) ever detected around galaxy nuclei (Combes et al. 2013). Furthermore, a fast and collimated outflow has been detected in HCN(1-0) and CO(1-0) emission in the nucleus of Arp 220, extending up to 120 pc and reaching velocities up to ±840 km s −1 (Barcos-Muñoz et al. 2018).
In the prototypical Seyfert 2 NGC 1068, a clear molecular ouflow has also been detected, entrained by the AGN radio jets (Krips et al. 2011;García-Burillo et al. 2014). NGC 1068 is also the first case where the molecular torus is resolved, using the continuum and the CO(6-5), HCN(3-2), and HCO + (3-2) emission lines observed with ALMA (Gallimore et al. 2016;García-Burillo et al. 2016;Imanishi et al. 2016). The dynamics of the molecular gas in the NGC 1068 torus revealed strong noncircular motions and enhanced turbulence superposed on a slow rotation pattern of the disk. The AGN is clearly off-center with respect to the torus, implying an m = 1 perturbation . Recently, we reported observations of molecular tori around massive BHs in a sample of seven nearby LLAGN (Seyfert/ Low Ionization Nuclear Emission-line Regions-LINERs), at the unprecedented spatial resolution of 3-10 pc (Combes et al. 2019, hereafter Paper I). The ALMA observations bring a wealth of new information on the decoupled molecular tori, which are found to have radii ranging from 6 to 27 pc, to be unaligned with the orientation of the host galaxy, and to frequently be slightly off-center with regards to the AGN position. The kinematics of the gas inside the sphere of influence (SoI) of the central BH also allowed us to estimate the BH masses (M BH ∼ 10 7−8 M ).
In this paper, we present the combined ALMA cycle 3 and 4 observations in the CO(3-2) line of the Seyfert galaxy NGC 613, with a spatial resolution of 17 pc. These observations were part of the sample presented in Paper I, but the data analyzed here have a higher sensitivity, allowing faint broad wings to be detected that are usually associated to outflows. This object presents a nuclear trailing spiral and we discovered a molecular outflow in its nuclear region. Therefore, NGC 613 is a special case allowing us to study the complexity of fuelling and feedback mechanisms in AGN, and to perform a detailed analysis of the gas flow cycle in AGN. All relevant characteristics of NGC 613 are described below. Observations are detailed in Sect. 2 and results are presented in Sect. 3. The properties of the nuclear molecular outflow discovered in the very center of NGC 613 are discussed in Sect. 4. The interpretation of inflowing gas in terms of gravitational torques is discussed in Sect. 5, and conclusions are drawn in Sect. 6.
NGC 613 is a nearby barred SB(rs)bc galaxy (de Vaucouleurs et al. 1991) at a distance of 17.2 Mpc (1 = 83 pc). It has a largescale bar of r bar ∼ 90 with a position angle (PA) of 127 • and a secondary nuclear bar with PA = 122 • (Jungwiert et al. 1997;Seigar et al. 2018). Judging from NED, NGC 613 is only moderately inclined, with inclination ∼38 • (see also Sect. 3.4). Prominent dust lanes are visible along the large-scale bar and the
NGC 613 shows clear evidence of star formation, shock excitation, and AGN activity (Davies et al. 2017). Radio continuum observations show evidence for a collimated jet from the AGN and a nuclear ring with a moderate inclination of i ∼ 55 • (Hummel et al. 1987;Hummel & Jorsater 1992). The presence of the outflow has already been suggested by the high-velocity dispersion of the [Feii] line along the radio jet (Falcón-Barroso et al. 2014). The ring-like structure in Brγ emission in NGC 613, comprising "hot spots" of current massive star formation (Böker et al. 2008;Falcón-Barroso et al. 2014), indicates an ongoing starforming episode. We collect the main properties of NGC 613 in Table 1.
In Cycle 3, NGC 613 was observed (project ID: #2015.1.00 404.S, PI F. Combes) simultaneously in CO(3-2), HCO + (4-3), and HCN(4-3) for both the compact (TC, baselines 15 to 630 m) and the extended (TE, baselines 15 to 1400 m) configurations. The largest recoverable angular scale corresponding to the shortest baseline is about 12 . The TC configuration was observed in April 2016 with 40 antennas and an integration time, including calibration and overheads, of 20 min, providing a synthesized beam of ∼0 . 38. The TE configuration was observed in August 2016 with 41 antennas, total integration of 40 min and a synthesized beam of ∼0 . 14. The correlator setup, designed to simultaneously observe three lines, provided a velocity range of 1600 km s −1 for each line, but did not center the HCO + (4-3) and HCN(4-3) lines (200 km s −1 on one side and 1400 km s −1 on the other, which is adequate for a nearly face-on galaxy), and 1800 MHz bandwidth in the continuum.
The Cycle 4 observations were carried out in November 2016 and July 2017 (project ID: #2016.1.00296.S, PI F. Combes) at higher spatial resolution (∼6.5 pc) aiming at resolving the molecular torus. The tuning configuration of Band 7 was in the CO(3-2) and HCO + (4-3) lines and the continuum to avoid a restricted velocity range in the expected broader spectral lines towards the nucleus. The correlator setup was selected to center the CO(3-2) and the HCO + lines in the 2 GHz bandwidth. The compact configuration (TM2, baselines 19 to 500 m) was observed with 44 antennas for an integration time of 14 min and a synthesised beam of 0 . 31 and the extended (TM1, baselines 19 to 3100 m) was observed with 43 antennas for 1.2 h and a synthesised beam of ∼0 . 08.
The phase center of the observations was that of the nucleus (Table 1), with a single pointing covering a field of view (FoV) of 18 . The galaxy was observed in dual polarization mode with 1.875 GHz total bandwidth per spectral window, and a channel spacing of 0.488 MHz corresponding to ∼0.8 km s −1 , after Hanning smoothing. The flux calibration was done with radio quasars close to the position of the target in the sky, which are regularly monitored at ALMA, and resulted in 10% accuracy.
The data from Cycles 3 and 4 were calibrated and concatenated with the CASA software (version from 4.5.3 to 4.7.2, (McMullin et al. 2007), and the imaging and cleaning were performed with the GILDAS software (Guilloteau & Lucas 2000). In Paper I, we used only the most extended configurations (TM1+TE), but in this work we have combined all the configurations to improve the sensitivity. The analysis was performed in GILDAS together with python packages (radio-astro-tools, using the Hogbom method and a natural weighting in order to achieve the best sensitivity, resulting in a synthesised beam of 0 . 21 × 0 . 19 for the concatenated data cube. The spectral line maps were obtained after subtraction of the continuum in the uv-plane using the tasks uv_continuum and uv_subtract. The data cubes were then smoothed to 10 km s −1 (11.5 MHz). The total integration time provided an rms of 87 µJy beam −1 in the continuum, and 0.43 mJy beam −1 in the line channel maps per channel of 10 km s −1 (corresponding to ∼1 K, at the obtained spatial resolution). The final maps were corrected for primary beam attenuation. Very little CO(3-2) emission was detected outside the full width at half power (FWHP) primary beam. Due to a lack of very short baselines (<15 m), extended emission was filtered out at scales larger than 12 in each channel map. Since the velocity gradients are high in galaxy nuclei, the line measurements are not significantly affected; indeed the size in each velocity channel is not expected to be extended.

Continuum emission
Previous ALMA band 3 and 7 observations by Miyamoto et al. (2017) detect continuum emission from both the circumnuclear disk (CND) and the star-forming ring (250 < r < 340 pc). At 95 GHz with a ∼0 . 6 resolution, these latter authors found a continuum jet with PA = 20 • , which corresponds to the 4.9 GHz and 14.9 GHz jets (Hummel & Jorsater 1992), close to the minor axis of the ring. In the nucleus, the negative spectral index, α ∼ −0.6, is compatible with synchrotron emission, with a small fraction of free-free emission, while the index α ∼ −0.2 along the starforming ring could be from free-free emission (Miyamoto et al. 2017).
At our resolution of 0 . 2 (∼17 pc), the central continuum is resolved at 350 GHz, with some compact emission along the star-forming ring, as display in Fig. 2. The ∼2.2 mJy peak emission is detected at 25σ significance. We determine the peak of the continuum emission by fitting a circular Gaussian source in the uv-plane using the GILDAS uv_fit task. The fitted results in the central continuum emission for the flux is 2.4 ± 0.1 mJy and for the RA and Dec relative to the phase center are ∆RA = −0 . 587 and ∆Dec = −0 . 03219, with a relative uncertainty of ±0 . 003. These values are listed as references of the new adopted center and AGN position in Table 1. Within the error bar, the AGN position is consistent with the positions derived by Miyamoto et al. (2017) and from X-ray observations (Liu & Bregman 2005).

Molecular gas distribution and morphology
Figure 3 displays the CO(3-2) channel maps, with a velocity range of 450 km s −1 and a velocity resolution of 10 km s −1 . The channels show evidence for a regular velocity field in a ring at a radius ∼3.5 (300 pc), with two winding-arm structures coming from the NW and SE directions. These spiral arms coincide with the beginning of the dust lanes along the bar seen in the HST/F814W image (Fig. 1); they are the contact points between the tangent dust lanes and the ring.
At small radii, there is a spiral structure in the central channels (±100 km s −1 ) that can be more clearly seen in Fig. 4. The two-arm nuclear gas spiral at r 100 pc is trailing toward the center; this is discussed further in Sect. 5. Additionally, the velocity distribution is perturbed for channels in the range ±50-100 km s −1 for r ∼ 150 pc, and the morphology shows evidence  for a filamentary structure connecting the ring and the nuclear spiral.
We constructed the moment maps of the CO(3-2) line, clipping the emission at <5σ rms . The integrated intensity (zeromoment) map in the left panel of Fig. 4 shows that the CO emission follows the ∼300 pc star-forming circumnuclear ring. Miyamoto et al. (2017) mapped the ring in CO(1-0) and CO(3-2) with ALMA at 0 . 7 and 0 . 4, respectively, and found a clumpy ring that is globally regular but has spots of active and efficient star formation. In our maps at higher resolution, we find that the molecular emission in the ring is clumpy and incomplete, and coincides with the same star forming clumps observed in Br γ in the near-infrared (NIR; see Falcón-Barroso et al. 2014, and Sect. 3.6). The CO(3-2) emission peak of 25 Jy km s −1 beam −1 corresponds to the AGN position reported in Table 1. Within the CND, there is a clear trailing two-arm spiral structure. The ring reveals two breaks into two winding spiral arms, at NW and SE. The main morphological features are shown in the sketch of the galaxy in Fig. 5.
We superposed the CO(3-2) contours onto the HST maps in the F814W filter 1 shown in Fig. 1. This reveals a remarkable similarity in morphology: the molecular ring seen in the CO emission coincides with the dusty nuclear ring in the HST image, and the winding arms are the beginning of the characteristic dust lanes along the bar. At the very center ( 100 pc) however, the stellar and molecular morphologies are dissimilar.
In the middle panel of Fig. 4, the intensity-weighted velocity (first-moment) map shows a clear rotation pattern in the galaxy plane, with velocities peaking in the range ∼±200 km s −1 from the systemic velocity (v sys 1481 km s −1 ; see discussion below). The velocity distribution and morphology in the central 200 pc are more perturbed due to the filamentary streams and the 1 The HST image was aligned to the ALMA astrometry, and the peak emission in the HST image was recentered to the AGN position in Table 1 Table 1. The color scale is in power stretch (with a power index of 0.5) ranging between 2 and 38 mJy beam −1 . nuclear spiral. The NW winding arm is mostly blueshifted and the SE redshifted, indicating that rotation and possibly gas pileup are taking place; this is discussed further in Sect. 5.
Close to the AGN, the velocity dispersion is high (σ ∼ 130 km s −1 ), as displayed in the right panel of Fig. 4 (secondmoment map). In the nuclear spiral the velocity dispersion ranges from 70 ∼ 120 km s −1 and the average dispersion along the ring is ∼40 km s −1 , with more elevated values in the clumpy regions. Furthermore, we can distinguish a disturbance in the overdense region in the west part of the ring, with an increased dispersion of 150 km s −1 . This region also corresponds to an enhanced spot observed in [Feii] with SINFONI, suggesting a strongly shocked medium (see also Sect. 3.6).

Carbon monoxide luminosity and H 2 mass
The mean intensity map is plotted in Fig. 4 (left). Since the galaxy is more extended than the primary beam, it is difficult to quantify the missing flux. We compare it to the central spectrum obtained with the 15 m single dish of the Swedish-ESO Submillimeter Telescope (SEST) in CO(1-0) and CO(2-1) over a 43 and 22 FoV, respectively. In Fig. 6, we display the total CO(3-2) spectrum integrated over the 18 FoV. Towards the central position, Bajaja et al. (1995) found a CO(2-1) spectrum peaking at T * A = 200 mK with FWHM = 300 km s −1 , yielding a total integrated flux of 1504 Jy km s −1 , in a beam of 22 . Their beam is very similar to our FoV of 18 . The flux comparison is relevant, since our FoV encompasses the entire nuclear ring, and the emission in this nuclear region corresponds to the strongest surface density at different wavelengths (Comerón et al. 2010;Ho et al. 2011;Li et al. 2011), as already discussed by Combes et al. (2014).
We assume a ratio of r 31 = T 3−2 /T 1−0 of 0.82, typical for Seyfert galaxies (Mao et al. 2010) and a ratio r 21 = T 2−1 /T 1−0 compatible with 1 within the error bars. The latter is derived from the SEST CO(2-1)/CO(1-0) observations in the galaxy center, by convolving the beam of CO(2-1) to 43 , implying a  Table 1). The synthesised beam of 0 . 21 × 0 . 19 is shown in red in the bottom left corner of the second-moment map. higher CO excitation at the center of NGC 613. This is expected for thermalized excitation and a dense molecular medium. In that case, the CO(3-2) flux should be higher than the CO(2-1), as we could presume the flux S ν ∝ ν 2 in the Rayleigh-Jeans approximation for gas at greater temperatures than 25 K and densities greater than 10 4 cm −3 . Using these values, the expected CO(3-2) intensity is ∼2266 Jy km s −1 in a 22 beam. When integrated over the spectral range (FWHM ∼ 250 km s −1 ), the integrated emission in our ALMA FoV of 18 , shown in Fig. 6 is 1307 Jy km s −1 . Therefore, we should expect some missing flux by a factor up to ∼40-50%, taking into account the uncertainties of the r 31 and r 21 ratios. We thus find a molecular  Table 2. mass of 5.6 × 10 8 M in our FoV, assuming thermally excited gas and a Milky-Way-like CO-to-H 2 conversion factor of 2 × 10 20 cm −2 /(K km s −1 ); see for example Bolatto et al. (2013).
The total CO(3-2) emission line profile integrated over the observed map (FoV of 18 ) is shown in Fig. 6. We decomposed the spectrum in three components, C1, C2, and C3, and A33, page 6 of 19 A. Audibert et al.: ALMA captures feeding and feedback from the active galactic nucleus in NGC 613 Notes. Results of the Gaussian fits for the three velocity components (C1, C2 and C3) shown in Fig. 6. (a) Peak flux. Fig. 7. Rotation curve of NGC 613. The orange circles represent the CO kinematics from our ALMA observations and the blue diamonds and stars are the Hα measurements by Burbidge et al. (1964), for a PA of 121 • and 115 • , respectively. The dotted lines are the best fit assuming the gas to be on circular orbits in a plane, v c = Ar/(r 2 + c 2 ) p/2 . the results of the Gaussian fits for each component are displayed in Table 2. The total flux is S CO(3−2) = 1307 ± 121 Jy km s −1 .

CO(3-2) kinematics
In Fig. 7 we show the rotation velocities deduced from the CO kinematics and the Hα rotation curve taken from Burbidge et al. (1964). As already pointed out, the dominant feature in the velocity field of the molecular gas appears to be due to circular rotation in the disk (middle panel of Fig. 4), and is consistent with the Hα kinematics. We find good agreement with the PA and inclination from optical studies, in the range of PA = 111 ∼ 124 • and i = 36 ∼ 47 • (Burbidge et al. 1964;Blackman 1981;de Vaucouleurs et al. 1991), and therefore we adopted the values of PA = 120 • and i = 41 • listed in Table 1.
We assume a simple model proposed by Bertola et al. (1991) for the rotation curve, assuming the gas is on circular orbits in the plane v c = Ar/(r 2 + c 2 ) p/2 , where A, c, and p are parameters of the model, and for p = 1 the velocity curve is asymptotically flat and for p = 3/2 the system has a finite total mass, and therefore we expect 1 ≤ p ≤ 3/2. Figure 7 shows the result of fitting the radially averaged velocities to this model; the dotted lines represent the best fit for circular velocity, v c . We can refine the above fit of the radially averaged velocities by fitting the entire velocity field with the same model by Bertola et al. (1991). The observed radial velocity at a position (R, Ψ) on the plane of the sky can be described as where θ is the inclination of the disk (with θ = 0 for a face-on disk), Ψ 0 is the PA of the line of nodes, v sys is the systemic velocity, and R is the radius. We used the tilted-ring model (Rogstad et al. 1974), which consists in dividing the velocity field into concentric rings in radii ∆r, with each ring being allowed to have an arbitrary v c , i, and PA. For each radius, we can independently fit the parameters of Eq.
(1) to the observed velocity field. We show the results of the fitting of the tilted-ring to the velocity map in Fig. 8. We adopted a ∆r = 0 . 3, which corresponds to the deprojected resolution of our observations in the galaxy plane. In the right panel, we display the residuals after subtracting the Bertola et al. (1991) model from the velocity field. As can be seen, the model represents the observed velocity field relatively well, with no significant amplitudes in the residuals along the galaxy disk, except in the west part of the ring where there is an important contribution from noncircular motions. This region coincides with the contact point between the ring and the SE winding arm, and it is also probably perturbed by shocks, as suggested by the high-velocity dispersions in the molecular gas and an enhancement in the [Feii] emission (Fig. 13). These noncircular motions along the minor axis and on the winding spiral indicate streaming motions associated with inflow, if we assume both trailing perturbations and that the north side is the far side of the galaxy, as the NW arm is blueshifted and the SE arm is redshifted.
An additional method was used to derive the CO kinematics using the "3D-Based Analysis of Rotating Objects from Line Observations" ( 3D BAROLO) software by Di Teodoro & Fraternali (2015). 3D BAROLO performs a 3D tilted-ring modeling of the emission line data cubes to derive the parameters that best describe the kinematics of the data. We ran 3D BAROLO on the CO(3-2) data-cube in order to investigate non circular motions, since the code allows us to infer radial velocities in the fit of the rotation curves. We performed several tests running the code, varying the fixed and free parameters of the disk model, and the results that better reproduce the observed velocity field are found when fixing the PA, the inclination, and the central position to 120 • , 41 • , and the position of the AGN, respectively. The best fit reveals radial components of the order of v rad ∼20 km s −1 in the nuclear region and v rad ∼20-100 km s −1 from the end of the circumnuclear ring at 4 up to r ∼ 7 , corresponding to the streaming motions of the winding arms. Overall, we find that in the case of high-resolution ALMA observations, the results using 3D BAROLO are relatively coherent with the 2D approach described above.
The position-velocity diagrams (PVDs) along the major (PA = 120 • ) and minor axis (210 • ) of NGC 613 are shown in Fig. 9. The best fit obtained with 3D BAROLO including radial velocity components is shown in dashed lines. We notice highly skewed kinematics in the center, with velocities 200 km s −1 that cannot be described only by co-planar circular motions in the galaxy disk. Since the velocities in the very center ( 0.5 ) strongly deviate from the rotation curve pattern, we believe this is the signature of an outflow emanating from the AGN. We discuss the detection and properties of the molecular outflow in Sect. 4.

Dense gas: HCO + , HCN, and CS emission
In ALMA band 3, dense gas is detected in various lines, namely HCN(1-0), HCO + (1-0) and CS(2-1), while SiO(2-1) was marginally detected at the edges of the radio jets, probably indicating the existence of shock regions related to the jets, as reported by Miyamoto et al. (2017). Along with CO(3-2), our  band 7 observations detected the HCO + (4-3), HCN(4-3), and CS(7-6) emission lines. The corresponding maps are displayed in Fig. 10, and the integrated spectra in Fig. 11. The CS(7-6) line is mainly detected in the center of the galaxy, while for HCO + (4-3) and HCN(4-3), we detect stronger emission in the center but also some clumps along the star-forming ring. The dense gas tracers detected in the very center are interpreted as a molecular torus, as discussed in Paper I. The HCN line is about twice brighter than the HCO + line in the nuclear region, which is the typical value expected for an AGN (e.g., Kohno et al. 2005;Krips et al. 2008;Imanishi et al. 2016). The detection of these lines reveals the presence of dense gas, since the critical densities of the HCO + (4−3), HCN(4-3), and CS(7-6) transitions are 2.6 × 10 6 , 1.4 × 10 7 , and 3.4 × 10 6 cm −3 , respectively.
As already mentioned, NGC 613 exhibits different ionization mechanisms. Recently, spatially resolved Baldwin, Phillips, & Terlevich (BPT) diagrams (Baldwin et al. 1981) were able to isolate the contributions from star formation, shock excitation, and AGN activity using optical line ratios, as studied with the Siding Spring Observatory Wide-Field Spectrograph (WiFeS) by Davies et al. (2017) and the European Southern Observatory Multi-unit Spectroscopic Explorer (MUSE) in Gadotti et al. (2019). Likewise, in the submillimeter domain, multi-line observations of higher J transitions are fundamental for the derivation of the chemical conditions of the molecular gas and the heating mechanisms (Viti et al. 2014;Imanishi et al. 2018). One useful extinction-free energy diagnostic tool in the centers of galaxies is the submillimeter-HCN diagram (Fig. 12) proposed by Izumi et al. (2016). The diagram uses the HCN(4-3)/HCO + (4-3) and HCN(4-3)/CS(7-6) ratios to distinguish the dominant energy source exciting the molecular gas in galaxies, whether by an AGN (XDR or X-ray-dominated region) or star formation (PDR or photodissociation region). The authors suggest enhanced integrated intensity ratios in circumnuclear molecular gas around AGN compared to those in starburst galaxies (submillimeter HCN enhancement).
Thanks to our high-resolution ALMA observations, we are able to disentangle the emission coming from the nuclear region within the spiral trailing feature observed in the central ∼100 pc and the contribution from a star-forming clump observed in all three molecular tracers at ∼250 pc from the nucleus. The clump is indicated as a circle in Fig. 10. We measured the line intensity ratios R + HCN/HCO and R HCN/CS in these two regions (CND/AGN, for a central aperture of 0 . 3 and clump; listed in Table 3). , and CS(7-6) in the panels from left to right, with the CO contours overplotted. HCN and HCO + present clumpy emissions along a radius of ∼3 , coinciding with the star-forming ring. All the dense gas tracers are detected in the very center as nuclear compact decoupled disks, or molecular tori as discussed in Paper I. The purple circles indicate the position of a clump detected in the HCN, HCO + , and CS emission and used to calculate the ratios in the submillimeter diagram in Fig. 12. The synthesized beam sizes are shown in red in the bottom left corner of each panel. Fig. 11. HCN(4-3), HCO + (4-3), and CS(7-6) emission line profiles (from left to right) integrated along our FoV, using the mask from the zeromoment map in Fig. 10. The integrated values fitting a Gaussian are listed in Table 2.  Izumi et al. (2016) for the high-resolution observations (spatial resolution <500 pc) using the line intensity ratios R + HCN/HCO and R HCN/CS . The red circles represent the measurements of galaxies hosting an AGN and the blue squares indicate the ratios found in starburst galaxies. We include the line ratios of NGC 613 (diamonds) measured in the CND, referred to here as AGN and in a clump detected ∼250 pc northeast of the central position in all the dense tracers shown in Fig. 10. Notes. Ratios of the HCN, HCO + , and CS lines in the CND region and in the clump shown in Fig. 10.
For a 5σ threshold, the clump is barely detected in CS(7-6) emission, and we present the line ratios for a 3σ clump detection and an upper limit for the R HCN/CS .
The line ratios are plotted on the submillimeter diagram of Izumi et al. (2016), as displayed in Fig. 12. We find that the CND region lies in the AGN-dominated part of the diagram, while the clump in the star-forming ring of NGC 613 is indeed dominated by star formation. Ultimately, we do find that the nuclear region of NGC 613 presents line ratios that indicate excitation conditions typical of XDRs in the vicinity of AGN. However, the line ratios in the center could also be explained by shocks driven by the radio jet entraining the gas, as HCN can become enhanced in  Table 1. shocks and this scenario is also suggested by the [Feii] enhancement discussed in Sect. 3.6.

Comparison to the warm molecular and ionized gas
In order to compare the CO(3-2) morphology with the ionized material and the warm molecular gas, we superposed the CO contours onto the NIR maps of [Feii], Brγ, and H 2 λ2.12µm presented in Falcón- Barroso et al. (2014), shown in Fig. 13. There is a remarkable resemblance between the ionized and warm molecular gas and the CO emission along the star forming ring. The positions and ages of the hot spots in the ring suggest a "pearls on a string" scenario of evolution of star formation as proposed by Böker et al. (2008). In this scenario, star formation only occurs in particular overdense regions and the young clusters move along the ring, following the gas movement, and simultaneously age, resulting in an age gradient along the ring. The expected sequence of star formation was indeed observed in the southern part of the ring: the hottest stars are found near the contact point of the dust lanes, and then fewer hot stars are found along the ring (see Fig. 8 of Böker et al. 2008).
The CO(3-2) emission presents lower intensities in the starforming hot spots at the contact points of the winding arms and the ring, corresponding to the hottest and youngest star-clumps traced by the HeI emission in Fig. 8 of Böker et al. (2008). This anti-correlation can be expected between the SFR and the gas surface density on small scales (<200 pc), since the gas has probably been consumed in the formation of new stars and the Kennicutt-Schimidt law (Schmidt 1959;Kennicutt 1998a) breaks down on scales of giant molecular clouds (Onodera et al. 2010). However, we do observe some dense clumps in HCO + and HCN emission, corresponding to other less hot regions of HeI emission (blue regions in their Fig. 8), indicating that some dense material is still being consumed in the southeast and northern part of the ring.
At the center, the nuclear spiral corresponds to the massive reservoir of the bright warm H 2 and [Feii], while contrasting with the weak emission in Brγ. Indeed, as discussed by Falcón- Barroso et al. (2014), the high [Feii]/Brγ ratio in the center indicates that excitation is dominated by shocks and photoionization in the nucleus, and follows the correlation between the strength of the [Feii] and 6 cm radio emission in Seyfert galaxies (Forbes & Ward 1993).
The   Colina et al. (2015). Apertures in the circumnuclear ring are shown with blue circles and the aperture from the nuclear region is shown by a gray diamond. Contours denote regions for young star formation, supernovae, and compact AGN, while solid lines denote upper limits for young star formation and AGN, both derived from integral field spectroscopy data (Colina et al. 2015). Dashed lines denote upper limits for star formation and AGN derived from slit spectroscopy (Riffel et al. 2013). emission originates in supernova-driven shocks. Indeed, the large ratio in the nucleus is similar to those found in AGN, indicating that the most likely mechanism for the production of [Feii] emission is shock excitation from the radio jets and/or supernova remnants typical of Seyfert galaxies (Rodríguez-Ardila et al. 2004). It is well-known that X-ray emission, which is dominant in Seyferts, can penetrate deeply into atomic gas and create extended partly ionized regions where [Feii] can be formed. Models presented by Alonso-Herrero et al. (1997) show that X-rays are able to explain [Feii]/Brγ ratios of up to ∼20, in agreement with the values observed in the NGC 613 nucleus. Colina et al. (2015) developed a 2D diagnostic diagram using integral field spectrograph data to characterize line-emitting regions. The diagnostic uses the line ratios H 2 (1-0)S(1) λ2.122 µm/Brγ and [Fe ii] λ1.644 µm/Brγ. Using this method, these latter authors found that young star-forming regions, older supernova-dominated regions, and the compact AGN-dominated A33, page 10 of 19 region occupy different areas in the line-ratio space. We show the NIR diagnostic diagram proposed by Colina et al. (2015) in Fig. 14. We use the line emission listed in Tables 1 and 2 of Falcón- Barroso et al. (2014) for different apertures: one in the nucleus, seven along the circumnuclear ring (spots 1 to 7), and one aperture between the ring and the nucleus (spot 8; see Fig. 1 of Falcón-Barroso et al. 2014). Placing the line ratios for the different apertures in the diagrams, we note that all spots in the circumnuclear ring are located in the regions of young star formation or SNe-dominated stellar populations. The spot in the intermediate region has higher line ratios, and is shifted into the compact AGN region. On the other hand, the nucleus of NGC 613 takes its place in the shock ionized LINER regime, where we expect H 2 /Brγ ratios greater than six (Mazzalay et al. 2013;Riffel et al. 2013). As suggested by Falcón- Barroso et al. (2014), the enhancement of the H 2 emission in the nucleus is possibly due to the interaction with the radio jet. The high value of H 2 /Brγ measured in the nucleus of NGC 613 is consistent with a LLAGN Seyfert/LINER composite, which is strongly influenced from shock heating.

The CO(3-2) kinematics of the outflow
The PVDs along the major and minor axes of the galaxy shown in Fig. 9 illustrate the highly skewed kinematics in the central part of NGC 613. The nuclear region contains gas reaching velocities up to ∼±300 km s −1 in projection, which is much higher than the rest of the nuclear disk gas, as we can see in the middle panel of Fig. 4.
Examination of the individual spectra at pixels around the AGN position reveals conspicuous blue-and redshifted wings in all the spectra within a radius of r ∼ 0.3 . These broad line wings are characteristic of gas ejection out of equilibrium; in this case an indication of a molecular outflow. The integrated spectrum extracted in a central circular aperture of 0.28 (∼23 pc) is shown in Fig. 15. We can see a bump in the emission around v ∼ 400 km s −1 , which corresponds to the isotope H 13 CN(4-3).
One might question whether the outflow signature could instead arise from a central mass able to induce such fast rotation in the center. If we assume that these high-velocity components are within the SoI, in a radius equivalent to our nuclear aperture of ∼0.3 = 25 pc, or about r = 38 pc in the galaxy plane, a massive black hole located in the center should have a mass of at least M BH = v 2 R/G for the rotational velocity in the galaxy plane v = 460 km s −1 (corresponding to ±300 km s −1 projected velocities), or M BH = 1.9 × 10 9 M . This value is about two orders of magnitude higher than the values reported in the literature; for example, M BH = 7.4 × 10 6 M derived by Davis et al. (2014) using the spiral pitch angle and M BH = 4 × 10 7 M using the central stellar velocity dispersion by van den Bosch (2016). In Paper I, we derived a BH mass of M BH = 3.7 × 10 7 M within the SoI of 50 pc, and for the whole NUGA sample the values derived tend to follow the pseudo-bulge region in the M BH − σ plane (Ho & Kim 2014).
Furthermore, if these high-velocity features were due to the rotation, they would not be observed along the minor axis of the galaxy; however, the PVDs in Fig. 9 clearly show a gradient of high-velocity emission along the minor axis that cannot be explained as being due to co-planar rotation. Could these features be a signature of inflowing gas? We cannot exclude a priori that the high-velocity components are due to a coplanar inflow, since the blueshifted gas is observed in the north and the We fit a Gaussian that takes into account the main disk (or "core") contribution to the CO emission and subtract this from the observed spectrum (residuals in gray). The regions considered in the computation of the molecular outflow properties in the blue and red wings are shown in color. Bottom panel: zoom view of the blue (−400 to −120 km s −1 ) and red (120 to 300 km s −1 ) wings. The red wing is integrated only up to +300 km s −1 to avoid the contamination from the H 13 CN emission at ∼400 km s −1 . redshifted gas is seen in the south. If the far side of the galaxy is the north side, this gradient could be explained as an inflow. However, the deprojected velocities onto the galaxy disk would be v ∼ 300/sin(i) ∼ 460 km s −1 ; the order of magnitude of the co-planar inflow would thus be too large, and therefore we can exclude this hypothesis.
Additionally, in Fig. 16 we show the contours for the blue and red wing emission, and we see that the contours overlap, the blue component arises in the northern part of the nucleus, and the red component in the south. Due to the small size of the outflow, we can barely resolve each wing contribution. However, there is an indication that they do not follow the rotation pattern of the mean velocity field. The direction of velocities are also opposite to what is found in the molecular torus (cf. Paper I).
We compare the molecular to the radio emission observed with the Karl G. Jansky Very Large Array (VLA) at 4.86 GHz (Hummel & Jorsater 1992) in Fig. 16. The VLA radio jet emission is shown in the right panel with the molecular outflow blue and red wing emission contours overplotted. The molecular outflow emission coincides with the central blob of the radio jet and appears to be aligned with the orientation of the radio jet at PA = 12 • .
In order to derive conservative values for the flux related to the outflow in the broad emission, instead of fitting one Gaussian component for the "core" and one broad component for the outflow (including some low-velocity emission that might not be associated with the outflowing material), we have taken two different approaches. First, we fitted a Gaussian to the nuclear spectrum to take into account the contribution of the main rotating disk (black line in Fig. 15). The residual spectrum is shown by the gray line, after subtracting the core contribution from the CO spectrum. We then integrated the contribution of the blue and  Notes. (a) Results of the Gaussian fit for the main disk contribution of the nuclear spectrum extracted in a r = 0.28 aperture shown in Fig. 15. (b) In this case we assumed that S peak is the maximum flux in the nuclear spectrum within the selected velocity range of the blue and red wings. red regions in the residual spectrum. The results of the fit of the main core and blue and red wings in the residuals are listed in Table 4. The derived molecular mass corresponding to the main disk component inside the radius of 23 pc is 4.8 × 10 7 M , which agrees with the mass of the 14 pc molecular torus of 3.9×10 7 M found in Paper I (using the same conversion factor and r 31 as in Paper I).
The second approach consists in creating moment maps only taking the velocity channels from −400 to −120 km s −1 for the blue component and from +120 km s −1 to 300 km s −1 for the red wing. We try to avoid the contribution from the H 13 CN by limiting the velocity channels up to 300 km s −1 . We display the blue and red wing contour maps in Fig. 16, superimposed over the intensity weighted velocity map of the original CO map, which represents the mean velocity pattern. The integrated fluxes corresponding to velocity intervals of the red and blue maps are also listed in Table 4. The peak temperatures of the blue and red wings represent ∼5% of the peak of the main core component, however they are still detected with a 50σ significance in the case of the wings in the residuals.
The maximal velocities of the red and blue components are up to about ±300 km s −1 in projection (∼460 km s −1 if in the galaxy plane, i = 41 • ). Given their location near the nucleus, we tentatively interpret these high-velocity features as the two sides of an outflow. Globally, these features represent as much as 8% of the total molecular emission in the nuclear ring region.

CO-to-H 2 conversion in the nuclear region of NGC613
From the integrated flux, S CO ∆V(Jy km s −1 ), listed in Table 4, we can derive the molecular mass involved in the outflow using the equation from Solomon & Vanden Bout (2005): where ν rest = 345.796 GHz, and D L is the luminosity distance in Mpc. The molecular mass, including helium, is then derived from M(H 2 ) = α CO L CO r 13 (Tacconi et al. 2013). This implies a molecular mass of M out = 1.9−2.2 × 10 6 M . This mass was obtained using r 31 = 0.82, a luminosity distance of 17.2 Mpc, and the standard Galactic CO-to-H 2 conversion factor (α CO,MW = 4.36 M (K km s −1 pc 2 ) −1 Dame et al. 2001;Bolatto et al. 2013). However, this mass could be an upper limit if the flow is made of more diffuse optically thin gas. The standard α CO = 4.36 for the Milky Way is the recommended value to use in the inner disks of galaxies. However, several observational works (e.g., Israel 2009a,b;Sandstrom et al. 2013) found that in the center of galaxies (R 1 Kpc) the conversion X CO can be up to three to ten times lower than X CO,MW . As pointed out by Bolatto et al. (2013), the recommended value to be applied in galaxy centers is α CO,cen ∼ 1 4 α CO,MW , with a 0.3 dex uncertainty. In our case, the masses involved in the outflow would be four times lower, that is, in the range of M out = 4.8−5.5 × 10 5 M , providing a more conservative estimative of the mass.
The assumption of a smaller α CO has already been discussed in the literature. In the case of the outflow in NGC 1068, the A33, page 12 of 19 1 4 α CO,MW factor was also assumed , in agreement with the Large Velocity Gradient (LVG) approximation analysis of the CO line ratios in the central region of this galaxy (Usero et al. 2004). Another example is the molecular outflow detected in M 51, where the authors assumed α CO = 1 2 α CO,MW (Querejeta et al. 2016;Matsushita et al. 2007). In a study of molecular gas excitation in the jet-driven winds of IC 5063, Dasyra et al. (2016) found that the outflowing molecular gas is partly optically thin, implying a α CO that is one order of magnitude smaller than that of the Galaxy.
Most molecular outflows are detected in ultra-luminous infrared galaxies (ULIRGs; Cicone et al. 2014). Therefore, the CO-to-H 2 conversion factor usually assumed in the literature is α CO = 0.8, approximately five times lower than the Milky Way factor α CO,MW . Since NGC 613 has a rather moderate IR luminosity (L IR = 3 × 10 10 L ), there is no reason a priori to adopt the lower factor applied to ULIRGs. In fact, we would like to highlight that uncertainties in α CO impact the comparison of scaling factors between outflows and host galaxies properties (e.g., Fiore et al. 2017;Fluetsch et al. 2019) by a factor of approximately five.

Mass outflow rate
To estimate the mass outflow rate, along with the observational quantities (outflow mass M out , size R out and velocity v out ), we need to assume a certain geometry. Following Fiore et al. (2017) and Cicone et al. (2014), for a spherical or multi-conical geometry, in which the outflowing clouds are uniformly distributed along the flow, the mass outflow rateṀ out can be calculated as: M out = 3v out (M out /R out ) . (3) If instead we assume a time-averaged thin expelled shell geometry (Rupke et al. 2005b), also adopted in the study of molecular outflows in the local Universe (Veilleux et al. 2017;Fluetsch et al. 2019), we havė which corresponds to the outflow mass averaged over the flow timescale, t flow = R out v out . The difference in the mass loading factor between the two proposed scenarios for the outflow geometry is a factor three times larger in the multi-conical/spherical description.
In the following estimates, we use Eq. (4) to derive more conservative outflow energetics, since the observations cannot constrain the geometry, and we are not able to favor one scenario over the other.
As discussed in Sect. 4.1, the outflow is found in a region of R out = 0 . 28 (∼23 pc) and here we consider the maximum projected velocity of the wings, v out = 300 km s −1 . If the outflow direction is between the observer line of sight and the galaxy plane, even assuming the maxima projected velocities, the de-projected velocities will encompass the adopted value, and therefore v out = 300 km s −1 is a conservative value. The flow timescale is consequently t flow ∼ 10 4 yr, which is comparable to the timescales of the BH growth-burst episodes of nuclear activity, with a duration of 10 4−5 yr (Wada 2004).
For an outflow mass of M out = (1.9−2.2)×10 6 M (assuming the standard α CO,MW ), we find a mass outflow rate ofṀ out = (25 − 29) M yr −1 . If instead, we use the mass derived assuming the typical values for galaxy centers, α CO = 1 4 α CO,MW , we find that the mass load rate isṀ out ∼ 7 M yr −1 .

The nuclear molecular outflow in dense gas tracers
We also detected the presence of broad wings in the nuclear spectrum of dense gas tracers HCN(4-3), HCO + (4−3), and CS(7−6), as indicated in Fig. 17. We also show the CO nuclear spectrum for comparison, and the high-velocity components cover the same velocity width of the CO wings (±300 km s −1 ). We find that the ratio between the peak fluxes for the main core in the dense gas and CO are ∼5, ∼8.5, and ∼36, for HCN(4-3), HCO + (4−3), and CS(7-6), respectively. We also see some indication in Fig. 17 that, at least for HCN(4-3) and CS(7-6), the line ratios of the blue wings tend to be higher than the core.
In order to quantify the line ratio in the core and wings, we show the HCN(4-3)/CO(3-2) along the nuclear spectra in Fig. 18. The ratio is shown up to velocities of +200 km s −1 due to the tuning of Cycle 3 described in Sect. 2. The core component, defined by the disk rotation with velocities up to ±100 km s −1 , has a ratio of about 0.2 and the ratio increases for the high velocities towards the wings, up to values of ∼0.6, suggesting an enhancement of HCN in the outflow. As discussed in Sect. 3.5, from analysis of the dense gas ratios in the submillimeter-HCN diagram  in Fig. 12, the nuclear region of NGC 613 presents excitation conditions typical of XDRs in the vicinity of AGN. However, we find evidence that the HCN in the outflow can be three times higher than the values found in the nuclear CND.
A similar trend was also reported in the molecular outflow of the QSO galaxy Mrk 231. The detection of the outflow in Fig. 18. HCN(4-3)/CO(3-2) ratio in the nuclear spectra. The ratio is shown up to velocities of +200 km s −1 due to the tuning of Cycle 3 described in Sect. 2. We can see that in the core, the ratio is ∼0.2 and increases towards the wings up to values of ∼0.6, indicating that the outflow is entrained mostly in a dense gas (n 10 4 cm −3 ), as discussed in Aalto et al. (2012).
HCN(1-0) by Aalto et al. (2012) covers the same velocity range (±750 km s −1 ) of the CO(1-0) outflow (Feruglio et al. 2010), and they found a high ratio of the HCN/CO∼0.3-1 in the outflow, higher than in the line core. The HCN is enhanced in the line wings by factors of between two and five, and they suggest that the outflow is mostly entrained in dense gas n 10 4 cm −3 , which is consistent with the molecular gas being compressed and fragmented by shocks (Aalto et al. 2012). High-resolution observations of HCN and HCO + in the higher J = 3 → 2 transition exhibit prominent, spatially extended line wings for HCN(3-2) in Mrk 231 (Aalto et al. 2015). In Mrk 231 there were no line wings detected in HCO + (3-2), while in NGC 613, there is some indication of high-velocity gas. The HCN(4−3)/HCO + (4−3) ratios in the core are approximately equal to two, while in the wings can be two to three times higher, similar to what is found in Mrk 231 using the HCN(1-0)/HCO+(1-0) ratios (Lindberg et al. 2016). Aalto et al. (2015) claimed that the elevated HCN abundance in the outflow is possibly caused by high temperatures in the X-ray-irradiated gas regions surrounding AGN (Harada et al. 2013).
Another possibility to explain the HCN enhancement in the outflow of NGC 613, is that the HCN emission stems from shocks potentially originating from the interaction of the ouflowing gas with the radio jet. The fact that the molecular outflow is spatially aligned with the central blob of the radio jet detected by Hummel & Jorsater (1992) (see Fig. 16), a region where there is evidence of shock excitation as discussed in Sect. 3.6, corroborates this scenario. A further detailed analysis of the line ratios for the whole NUGA sample will be discussed in a future paper.

Driving mechanism: AGN or star formation?
The origin of the outflow might be related to star formation, which is concentrated in the nuclear ring region. The SFR can be estimated from the IR luminosity and the calibration from Kennicutt (1998b). From the IRAS fluxes, the IR luminosity is L IR = 3 × 10 9 L (Table 1, Sturm et al. 2002), and the total SFR is equal to 5.3 M yr −1 . Based on the Hα luminosity associated with star formation in the central 3 × 3 kpc measured by Davies et al. (2017), L(Hα) SF = 5.28 × 10 41 erg s −1 , we can also deduce, from Kennicutt's calibration, SFR = 4.1 M yr −1 , which is consistent with the estimate from IR luminosity. If we isolate the contribution only coming from the circumnuclear starforming ring, the SFR estimated is SFR = 2.2 M yr −1 in the ring of r 400 pc (Mazzuca et al. 2008). Evidence of young starforming regions in the "hotspots" along the ring has recently been found by Falcón-Barroso et al. (2014). However, the value reported for NGC 613 in the nuclear region of r ∼ 40 pc in aperture is very low: SFR∼ 0.015 M yr (Falcón-Barroso et al. 2014). This low value could be due to extinction factors or a possible AGN contribution. The same nuclear region has a large reservoir of warm molecular gas (e.g., see Fig. 13), also found in other Seyfert galaxies (Hicks et al. 2013). Falcón-Barroso et al. (2014) suggest a cyclical episode of starburst about ∼10 Myr ago, followed by another episode of nuclear activity.
We estimate a mass rate ofṀ out ∼ 27 M yr −1 for the nuclear molecular outflow. Although this estimate is uncertain by a factor of a few given the unknown projection and the assumptions previously discussed in the text, the SFR in the nuclear region is about three orders of magnitude lower thanṀ out . In general, galactic winds driven by starbursts correspond to mass-outflow rates of the same order as the SFR (e.g., Veilleux et al. 2005). Given the discrepancy between the SFR in the nuclear region and the mass-load rate of the outflow, we conclude that star formation alone is not able to drive the nuclear molecular outflow in NGC 613.
It has already been established that the mass-outflow rate increases with the AGN luminosity, supporting the idea of a luminous AGN pushing away the surrounding gas through a fast wind. Previous observational studies (Cicone et al. 2014;Carniani et al. 2015;Fiore et al. 2017) have shown that the molecular outflow properties are correlated with AGN luminosity, where the outflow kinetic power corresponds to about 5%L AGN and the momentum rate is ∼20L AGN /c, in agreement with theoretical models of AGN feedback (Faucher-Giguère & Quataert 2012; Zubovas & King 2012). For a sample of molecular and ionized outflows, Carniani et al. (2015) found that the ionized gas only traces a small fraction of the total gas mass, suggesting that the molecular phase dominates the outflow mass. This trend is also found by Fiore et al. (2017), but the ratio between molecular and ionized mass-outflow rates is reduced at the highest AGN bolometric luminosities. Nevertheless, these latter authors analyzed different samples of galaxies, and this conclusion could be affected by selection bias.
From XMM-Newton observations of NGC 613, Castangia et al. (2013) reported an X-ray luminosity of log L X (2−10 keV) = 41.3 erg s −1 . Applying a bolometric correction from Marconi et al. (2004) gives an AGN bolometric luminosity of L AGN,X = 1.7 × 10 42 erg s −1 . The bolometric luminosity derived by Davies et al. (2017) using the [Oiii] emission associated only with the AGN contribution, traced as an extended ionization cone aligned with the radio jet, is L AGN,[Oiii] = 4 × 10 42 erg s −1 . If we include the shock and star forming contributions of the total O[iii] emission, we obtain L bol,[Oiii] = 3.75 × 10 43 erg s −1 . The shock contribution most likely arises from the radio jet launched by the AGN, but here we cannot disentangle the contribution from star formation; the latter probably overestimates the bolometric AGN luminosity, while the former probably sets a lower limit.
In a recent study, Fluetsch et al. (2019) identified 45 molecular outflows in the local Universe using previous results from the literature and new detections from ALMA archive. These latter authors propose an even tighter empirical relation between the mass-outflow rate and the SFR, stellar mass, M * , and the bolometric L AGN : log(Ṁ out ) = 1.14 log 0.52 SFR M yr −1 + 0.51 where the SFR is calculated from the total IR luminosity anḋ M out is in M yr −1 . We adopt the bolometric AGN luminosity derived from the O[iii] emission line, L AGN = 4 × 10 42 erg s −1 (Davies et al. 2017), the total SFR inferred from the IRAS fluxes, SFR = 5.3 M yr −1 , and the stellar mass of M * = 4.5 × 10 10 M calculated in Paper I, derived from the S 4 G3.6 µm IR image and the GALFIT decomposition (Salo et al. 2015). Hence, according to Eq. (5), we should expect an mass-outflow rate of 4.8 M yr −1 in NGC 613, which is six times lower than our estimate. The kinetic power of the nuclear outflow can be estimated as P K,out = 0.5v 2Ṁ out . For the α CO,MW assumption, M out ∼ 27 M yr −1 , and we find P K,out = 8 × 10 41 erg s −1 , which corresponds to 20%L AGN . This value exceeds the predictions from AGN feedback models and cosmological simulations that require that a fraction of the radiated luminosity be coupled to the surrounding gas ∼5%L AGN (Di Matteo et al. 2005;Zubovas & King 2012). However, these predictions are based on powerful AGN accreting close to the Eddington limit, and we hypothesize that the coupling efficiency between AGNdriven outflows and L AGN should be outweighed in LLAGN. The Eddington luminosity derived from the BH mass (M BH = 3.7 × 10 7 M , Paper I) is L Edd = 4.6 × 10 45 erg s −1 , leading to a low accretion rate L AGN /L Edd 1 × 10 −3 . If we assume the lower α CO = 1/4α CO,MW , the kinetic power of the outflow is P K,out ∼ 2 × 10 41 erg s −1 , and the coupling is 5%L AGN . In both assumptions, the results indicate that the AGN can power the nuclear outflow in NGC 613, but the former requires a higher coupling efficiency.
The momentum flux of the outflow can be computed from P out =Ṁ out v ∼ 5 × 10 34 dynes. Compared to the momentum provided by the AGN photons, L AGN /c = 1.3 × 10 32 dynes, it is higher by a factor ofṀ out v ∼ 400L AGN /c. In case of energyconserved winds, AGN feedback models predict a momentum boost by factors up to 50 (e.g., Faucher-Giguère & Quataert 2012). Even assuming the lower α O , the value would exceed the predictions (∼100L AGN /c). One possibility to explain the high apparent energetics of the outflow in NGC 613 is that the AGN activity was stronger in the past. As discussed by Fluetsch et al. (2019), these authors found that 10% of the galaxies in their sample present outflows that exceed the theoretical predictions. They suggest these galaxies could have fossil outflows, resulting from a strong past AGN that has now faded. In the case of the LLAGN NGC 1377, a collimated molecular outflow is detected at 150 pc scales , possibly entrained by a faint radio jet, and these latter authors suggest that the nuclear activity of NGC 1377 may also be fading.
Alternatively, it is possible that the outflow is driven by the AGN radio jets. The radio power of the jet can be estimated from the 1.4 GHz luminosity, P 1.4 , using the relation of Bîrzan et al. (2008). If we compute the P 1.4 luminosity using the NRAO VLA Sky Survey (NVSS) flux measurement of 179.6 mJy at ν = 1.4 GHz (Condon et al. 1998), we find a radio power of P jet = 1.2 × 10 43 erg s −1 . Since this value is elevated, it might include emission from the circumnuclear star-forming ring. To avoid this contribution, we used only the flux of the linear component of the radio jet at 4.86 GHz by Hummel & Jorsater (1992). At 4.86 GHz, the flux associated to the jet is S ν = 7.2 mJy, and using a spectral index of α = −0.8 to derive the flux at 1.4 GHz, we find P jet = 5.5 × 10 42 erg s −1 . Hydrodynamical simulations of the interaction of the jet with a clumpy interstellar medium has shown that the jet is able to drive a flow efficiently as soon as the Eddington ratio of the jet P jet /L Edd is larger than 10 −4 (Wagner et al. 2012). This condition is required for the jet-driven outflow velocity to exceed the velocity dispersion of the M − σ relation but it also depends on cloud size in the ISM. In NGC 613, this ratio is about 1.2 × 10 −3 , and the jet power is about one or two orders of magnitude higher than the kinetic luminosity of the outflow. We conclude that the jet is able to drive the molecular outflow, even with low coupling. The molecular outflow in NGC 613 is an intriguing case where a very powerful molecular outflow is detected in a LLAGN. The SFR is very weak in the nuclear region and therefore unable to drive the flow. Radio jets are found to play a role in the driving mechanism to accelerate the molecular gas in other LLAGN. This is the case for the LINER NGC 1266, as suggested by Alatalo et al. (2011), which shows the presence of a weak radio jet (Nyland et al. 2013), and NGC 6764 (Leon et al. 2007) and the Seyfert 2 NGC 1433 (Combes et al. 2013). The properties of the flow require a contribution from the AGN through the entrainment of its radio jets.

Torques and AGN fueling
In Paper I, we reported the ALMA observations of the CO(3-2) line for all the galaxies in the NUGA sample. In three cases, NGC 613, NGC 1566, and NGC 1808, the CO emission revealed a nuclear trailing spiral, as in the large-scale one. Inside the nuclear ring at the ILR of the bar, usually a leading spiral is expected, developing transiently and generating positive torques, which drive the inner gas onto the ring. However, when the gravitational impact of the black hole is significant, the spiral can then be trailing, and the torques negative, and the gas is inflowing to fuel the nucleus (Buta & Combes 1996;Fukuda et al. 1998).
We find filamentary structures in NGC 613, pointing to radial gas transport from the circumnuclear star-forming ring into the core region dominated by the AGN, as displayed in Fig. 4. The ring morphology appears disturbed by a radial outflow of material from the AGN, which is confirmed by the existence of a weak jet in archival radio maps. However, the radio jet does not seem to have any significant effect on the morphology of the large reservoir of molecular gas that has accumulated inside the central 100 pc.
If gas is inflowing from the bar dust lanes into the ring, as expected from gravity torques (García-Burillo et al. 2005), there is also an inflow in the CND, due to the nuclear trailing spiral, as already observed for NGC 1566 . Inside the nuclear spiral structure, there is a very dense and compact (radius ∼14 pc) rotating component, which might be interpreted as the molecular torus (Paper I).
In order to explore the efficiency of feeding in NGC 613, we estimated the gravitational torques exerted by the stellar potential on the molecular gas, following the methodology described by García-Burillo et al. (2005). The gravitational potential is computed in the plane of the galaxy using a red image (F814W) from HST, since 2MASS images have insufficient angular resolution. We did not separate the bulge from the disk contribution since NGC 613 is a late-type (Sbc) galaxy. This means that the bulge was effectively assumed to be flattened. Dark matter can be safely neglected inside the central kiloparsec. The image was rotated and deprojected according to PA = 120 • and i = 41 • , and then Fourier transformed to compute the gravitational potential and forces. A stellar exponential disk thickness of approximately one twelfth of the radial scale-length of the galaxy (h r = 3.8 kpc) was assumed, giving h z = 317 pc. This is the average scale ratio for galaxies of this type (e.g., Barnaby & Thronson 1992;Bizyaev & Mitronova 2009). The potential was obtained assuming a constant mass-to-light ratio of M/L = 0.5 in the I-band over the considered portion of the image of 1 kpc in size. This value is realistic in view of what is found statistically for spiral galaxies (Bell & de Jong 2001). The pixel size of the map is 0.06 = 5 pc, the average value between the ALMA and HST resolutions. The stellar M/L value was fit to reproduce the observed CO rotation curve.
The potential Φ(R, θ) can be decomposed into its different Fourier components: where Φ m (R) and φ m (R) are the amplitude and the phase of the m-mode, respectively. The strength of the m-Fourier component, Q m (R), is defined as Q m (R) = mΦ m /R|F 0 (R)|, that is, by the ratio between tangential and radial forces (Combes & Sanders 1981).
The strength of the total nonaxisymmetric perturbation Q T (R) is defined similarly with the maximum amplitude of the tangential force F max T (R). Their radial distributions and the radial phase variations are displayed in Fig. 19.
The derivatives of the potential yield the forces per unit mass (F x and F y ) at each pixel, and the torques per unit mass t(x, y) are then computed by t(x, y) = x F y −y F x . The sign of the torque is determined relative to the sense of rotation in the plane of the galaxy. The product of the torque and the gas density Σ at each Fig. 20. Top: map of the gravitational torque, (t(x, y)×Σ(x, y), as defined in the text) in the center of NGC 613. The torques change sign as expected in a four-quadrant pattern (or butterfly diagram). The orientation of the quadrants follows the nuclear bar's orientation. In this deprojected picture, the major axis of the galaxy is oriented parallel to the horizontal axis. Bottom: deprojected image of the CO(3-2) emission at the same scale and with the same orientation, for comparison. The color scales are linear, in arbitrary units. pixel allows one to then derive the net effect on the gas at each radius. This quantity t(x, y)×Σ(x, y), is shown in Fig. 20, together with the deprojected CO map.
The torque weighted by the gas density Σ(x, y) is then averaged over azimuth, that is, The quantity t(R) represents the time derivative of the specific angular momentum L of the gas averaged azimuthally (García-Burillo et al. 2005). Normalizing at each radius by the angular momentum and rotation period (T rot ) allows us to estimate the efficiency of the gas flow, as shown in Fig. 21. The torques are negative in the winding arms at r ∼ 500 pc (corresponding to Fig. 21. Radial distribution of the torque, quantified by the fraction of the angular momentum transferred from the gas in one rotation-dL/L, estimated from the CO(3-2) deprojected map. The vertical dashed line at 25 pc radius delimitates the extent of the central gas outflow, and the computation has no meaning here. The torque is negative inside the 100 pc nuclear spiral and in the winding arms and positive in the filaments. the dust lanes) and are indeed contributing to drive the outer gas into the star-forming ring at ∼350 pc. The observations of nuclear rings are more common among barred galaxies (Kormendy & Kennicutt 2004;Peeples & Martini 2006;Jogee et al. 2006), and can be explained by the gas slowing down as it crosses the ILR, consequently weakening the gravitational torques, and the gas piles up in rings (Combes & Gerin 1985;Byrd et al. 1994). We can see in Fig. 21 that the efficiency certainly drops in the inner ILR region. The filaments within the nuclear region, between 100 and 200 pc, show that the gas gains angular momentum in one rotation, which is T rot ∼ 12 Myr.
The nuclear bar strength is moderate, however the BH has a strong influence on the nuclear molecular gas, as shown in Fig. 21, and the fueling efficiency is high in the nuclear spiral. Between 25 and 100 pc, the gas loses its angular momentum in one rotation, which is T rot (100 pc) ∼ 9.5 Myr. Inside 25 pc, the torques are positive and correspond to the region where we detect the molecular outflow. Since the gas in this region is not in quasi-stationary orbits, but is rather ejected from the galaxy plane, the computation cannot be interpreted in terms of the average torque here.
As shown in Fig. 20, the nuclear spiral structure inside the ILR ring of the bar is of a trailing nature and is located inside the negative torque quadrants. This might appear surprising, since in many cases the spiral structure is predicted to be leading in this region, and the torque positive, maintaining the gas in the ILR ring (Buta & Combes 1996). However, in the presence of a sufficiently massive black hole, this behavior can be reversed: the spiral becomes trailing, and the torque negative. Indeed, the indicator of the precession rate of elliptical orbits in the frame of the epicyclic approximation, Ω − κ/2, is significantly modified by a central massive body. Instead of decreasing regularly towards zero at the center, the Ω − κ/2 curve increases steeply as r −3/2 . The gas undergoes collisions, loses energy and spirals progressively towards the center. When the precession rate decreases, the series of elliptical orbits precess more and more slowly, and then lag at smaller radii, forming a leading structure; when the precession rate increases, they form a trailing structure .
One way to detect spiral structure in the nuclear gas is to amplify the dust extinction features in the HST images.
However, we would like to stress that observations of spiral dust lanes interior to the star-formation rings are not very common (Martini et al. 2003;Peeples & Martini 2006). Notwithstanding, a clear case is seen in NGC 2207, where dusty spirals extend from ∼50 to 300 pc, suggesting that the gas continues to sink inside the ILR and might be promoting gas accretion onto the nucleus (Elmegreen et al. 1998). Multiple spiral arms were also detected inside the star-forming ring (r ∼ 700 pc) in the LINER/Sy 1 galaxy NGC 1097 by Prieto et al. (2005). These latter authors compared their results with hydrodynamic models of Maciejewski (2004), concluding that inflows were occurring along the nuclear spirals. Although nuclear star forming rings are common, their interior structure is difficult to distinguish in dust maps, making the molecular gas emission a better tracer. The NUGA maps show clear evidence of trailing spirals in three objects of the sample, highlighting the importance of the highresolution ALMA observations. NGC 613 is therefore an example of a trailing spiral inside the nuclear ring of a bar; for example, see the case described in Buta & Combes (1996). This means that the mass of the black hole should be sufficiently high to have an influence on the gas dynamics on a 100 pc scale. In summary, the gravity torques are negative in the winding spirals, and the gas accumulates in the star forming ring at the inner ILR of the nuclear bar. The filamentary structure gains angular momentum and then the nuclear spiral at 100 pc present a very high efficiency in fuelling the central BH.

Conclusions
Here we present the combined ALMA cycle 3 and 4 observations for the Seyfert/nuclear starburst galaxy NGC 613. The combined observations in CO(3-2) reach a spatial resolution of 0.2 ∼ 17 pc. We study the morphology and the kinematics of the molecular gas in the central 1 kpc and our main findings are summarized below: -The 350 GHz continuum map shows a compact, barely resolved emission peak at the position of the AGN, surrounded by a patchy ring, matching the star-forming ring seen on optical images. -The morphology of CO(3-2) line emission reveals several components: a two-arm nuclear spiral at r 100 pc trailing the gas toward the center, and a circumnuclear ring at ∼350 pc corresponding to the star-forming ring. Also, we find evidence of a filamentary structure connecting the ring and the nuclear spiral. The ring reveals two breaks into two winding spiral arms towards the NW and SE corresponding to the dust lanes in the HST images. -The kinematics of the CO emission show a rather regular rotational velocity field in the inner kiloparsec disk. We applied a tilted-ring model to fit the velocity map and the residuals indeed show the bulk of the molecular gas in circular motion except in the west part of the ring, the perturbations of which are due to the contact point between the ring and the SE winding arm. -There is a remarkable coincidence between the molecular gas and the warm H 2 and ionized gas traced by the [Feii] and Brγ emission lines in the star-forming ring. Line diagnostics in the NIR indicate that the clumps in the ring are in agreement with young star-forming excitation. On the other hand, the nucleus of NGC 613 presents an excitation mechanism typical of Seyferts or LINERs. -We measured the line intensity ratios R + HCN/HCO and R HCN/CS in the nuclear region and in a clump along the ring. We find that the ratio for the nuclear region points to the AGN-dominated part of the HCN-submillimeter diagram, while the ratio for the clump is located in the starburstdominated part. These results indicate that the nuclear region of NGC 613 presents line ratios in agreement with excitation conditions typical of XDRs in the vicinity of AGN.
-In the PV diagrams, we find skewed kinematics in the nuclear region of r ∼ 25 pc. This feature is seen as broad wings (v ± 300 km s −1 ) in the CO nuclear spectrum, and the wings are also present in the dense gas tracers. We identify this feature as a molecular outflow emanating from the nucleus. The molecular outflow is co-spatial with the central blob detected in the radio jet. -We derive the molecular mass associated to the outflow as M out = 2 × 10 6 M . The mass outflow rate isṀ out = 27 M yr −1 . If instead, we use the mass derived assuming the typical values for galaxy centers, α CO = 1 4 α CO,MW , we find that the mass load rate isṀ out ∼ 7 M yr −1 .
-We find a HCN enhancement in the outflow probed by an increase in the HCN(4-3)/CO(3-2) ratio along the wings in the nuclear spectra. While the core has ratios of about ∼0.2, this value increases in the wings up to ∼0.6, indicating that the outflow is entrained mostly in a dense gas (n 10 4 cm −3 ). Another possibility is that the HCN emission stems from shocks potentially originating from the interaction of the ouflowing gas with the radio jet. -The molecular outflow energetics exceed the values predicted by AGN-feedback models. The kinetic power of the nuclear outflow corresponds to P K,out = 20%L AGN and the momentum rate isṀ out v ∼ 400L AGN/c . We speculate that, given its current weak nuclear activity, NGC 613 might be a case of fossil outflows resulting from a strong AGN that has now faded. -The outflow can be entrained by its radio jet. We find that P jet /L Edd ∼ 10 −3 and the jet power is about one or two orders of magnitudes higher than the kinetic luminosity of the outflow. In these conditions, the jet is able to drive the molecular outflow. -The trailing spiral observed in CO emission is inside the ILR ring of the bar. We computed the gravitational potential from the stars within the central kiloparsec from the I-band HST image. Weighting the torques on each pixel by the gas surface density observed in the CO(3-2) line allowed us to estimate the sense of the angular momentum exchange and its efficiency. The gravity torques are negative from 25 to 100 pc. Between 50 and 100 pc, the gas loses its angular momentum in a rotation period, providing evidence of fueling of the AGN. The molecular outflow in NGC 613 is an intriguing case where a very powerful molecular outflow is detected in a LLAGN. The SFR is very weak in the nuclear region, and is therefore not able to drive the flow. The properties of the flow require the contribution of the AGN through the entrainment of its radio jets. On the other hand, there is a clear trailing spiral observed in molecular gas inside the ILR ring of a bar, indicating that the SMBH is influencing the gas dynamics. Instead of maintaining the ILR ring density, the torques are then driving gas towards the nucleus, a first step towards possibly fueling the AGN. NGC 613 is a remarkable example of the complexity of fuelling and feedback mechanisms in AGN, and reinforces the importance of detailed analysis of nearby galaxies with ALMA capabilities to shed light on the gas-flow cycle in AGN.