© ESO 2019

1. Introduction

Supermassive black holes (SMBH) are thought to be ubiquitous in galaxies with bulges and may be key to the formation and evolution of galaxies (Kormendy & Ho 2013). It has been argued that the correlation between the host bulge and central black hole mass (Ferrarese & Merritt 2000; Gebhardt et al. 2000; Tremaine et al. 2002; Ferrarese & Ford 2005; Gültekin et al. 2009; Kormendy & Ho 2013) implies a direct causal relationship between the accretion of material by the black hole, the host galaxy star formation, and AGN-driven feedback, but direct observational evidence of the mechanisms responsible has remained elusive (Heckman & Best 2014). Over the past decade there have been a growing number of facilities providing 3D spectroscopic imaging observations which support the study of gaseous and stellar kinematics in active and inactive galaxies at radio (e.g., Morganti et al. 2009; Nesvadba et al. 2010), infrared wavelengths (e.g., Storchi-Bergmann et al. 2010; Riffel et al. 2013; Diniz et al. 2015) and optical (e.g., Dumas et al. 2007; Storchi-Bergmann et al. 2007; Dicaire et al. 2008; Westoby et al. 2012; Schnorr-Müller et al. 2014a; Lena et al. 2015; Roche et al. 2016). The combination of enhanced sensitivity at the unprecedented spatial and spectral resolution provided by ALMA has opened a new window on molecular gas dynamics to study the central kiloparsec of local galaxies where the dynamical and AGN-activity timescales become comparable, and where nuclear fueling, AGN feedback, and host galaxy quenching can be probed directly.

NGC 1566, a nearly face-on barred spiral galaxy (morphological type SAB) is the dominant (de Vaucouleurs 1973) and brightest member of the Dorado group (Bajaja et al. 1995; Agüero et al. 2004; Kilborn et al. 2005), and one of the nearest and brightest Seyfert galaxies. Although it has many features of a Seyfert 1, several studies (Alloin et al. 1985; Bottema 1992; Ehle et al. 1996; Kawamuro et al. 2013) have indicated this nature as uncertain. NGC 1566 has an intermediate-strength bar (projected radius 33″ or ∼1.5 kpc and P.A. ∼0°; Agüero et al. 2004), and two strongly contrasted spiral arms. Both the assumption of trailing spiral arms and the more marked dust obscuration on the NW side (dust in the disk obscuring light from the bulge) seen in Hubble Space Telescope (HST) imaging (Malkan et al. 1998), point to the NW side as the near side and the SE as the far side of the disk.

Despite its proximity, the distance of NGC 1566 is controversial. Several studies using the Tully Fisher relation (TFR) have claimed distances between 18 Mpc (EDD¹) and around 6 Mpc (Sorce et al. 2014; Tully et al. 2013). Even though the H I spectra used in these studies have high signal-to-noise ratios, they clearly show a double-peaked structure, which could lead to significant underestimations of the rotation velocity and thus the TFR-based distance. Consequently, in this work, we use the mean distance of 10 Mpc from NED², in agreement with the distance used by Combes et al. (2014; hereafter C14). At this adopted distance, the linear scale in our images is 0. pc arcsec⁻¹.

The systemic velocity of NGC 1566 is 1504 km s⁻¹ from H I observations (NED), but there is a wide range in the optical spectroscopy-based recession velocity values found by different authors. C14 found a systemic velocity of 1516 km s⁻¹ for CO J:3-2; offset ∼12 km s⁻¹ from the H I- derived value. We find that the galaxy’s integrated CO J:2-1 profile is centered on a systemic velocity of 1485 km s⁻¹ (see below). As we discuss in this work, the nuclear CO profiles are highly perturbed and non-axisymmetric about the nucleus. Thus, the molecular gas derived systemic velocities do not necessarily trace the true systemic velocity of the nucleus of NGC 1566.

The position angle of the major axis of NGC 1566 is ∼45° (HyperLEDA; Makarov et al. 2014)³, and the inclination of the disk was found to be i = 35° (C14, Agüero et al. 2004). In this work we use this major axis position angle and a galaxy disk inclination of i = 33°, as derived from our CO data (Sect. 2).

H I studies of the local group of NGC 1566 (a sub-part of the Dorado group) show that NGC 1566 is interacting with its smaller companions (Kilborn et al. 2005), and this finding is reinforced by the strong correlation found between galaxies with prominent barred structures and companions in the Dorado group (Kendall et al. 2011). Given its strong and symmetric spiral arms, its active nucleus (AGN), and its proximity, NGC 1566 has been the subject of great interest within the community, and has extensive studies of its spiral arm formation (Korchagin et al. 2000; Ma 2001; Erwin 2004; Kendall et al. 2011), its gas kinematics (Pence et al. 1990; Bottema 1992; Bajaja et al. 1995; Agüero et al. 2004; Dicaire et al. 2008; Mezcua et al. 2015) and the feeding and feedback of the SMBH in its center (Elvis et al. 1989; Schmitt & Kinney 1996; Combes et al. 2014; Smajić et al. 2015; Davies et al. 2016; da Silva et al. 2017).

An early kinematic study of NGC 1566 in H I and Hα (Pence et al. 1990) showed that the most significant Hα kinematic feature (after subtraction of regular rotation) was a spiral arc located 26″ from the nucleus towards the main spiral arm on the southeast (far) side of the galaxy. This spiral arc has a redshifted velocity of 60 km s⁻¹, i.e., gas moving away from the nucleus under the assumption of motion in the disk of the galaxy. Under this assumption, Pence et al. (1990) estimated outflow velocities, most plausibly driven by the AGN, of 130 km s⁻¹ in the plane of the galaxy, i.e., an equatorial outflow. Schmitt & Kinney (1996) also supported the presence of an outflow when analyzing the morphology of the nuclear [O III] emission; they observed a total extension of ∼0$0\stackrel{″}{.}7$, mainly to the SE, which they interpreted as the base of a conical NLR originating in the nucleus and oriented perpendicular to the plane of the disk, i.e., a polar outflow. An extension in the nuclear [O III] emission to the SW was also found by da Silva et al. (2017) in integral field unit (IFU) imaging. They interpreted the morphology and blueshifted kinematics of the [O III] line as being consistent with an outflow of ∼500–800 km s⁻¹ driven by the AGN perpendicular to the plane of the disk. Agüero et al. (2004) found a H II deficiency in the inner regions of NGC 1566 (see also Pence et al. 1990), reinforcing the evidence of outflows to the SE, and posited that the blueshifted knot found ∼8″ from the nucleus on the far side of the disk signaled the presence of inflows along the galaxy minor axis. Using optical integral field spectroscopy, Davies et al. (2016) found a high dispersion (∼100–200 km s⁻¹) in the Hα line over a region ∼200 pc to the SW of the nucleus, which they interpret as most likely due to an outflowing gas illuminated by the radiation field of the AGN. They find that the latter is sufficiently high to drive outflows in this galaxy. In X-rays, Elvis et al. (1989) found extended X-ray emission centered at a position ∼10″ from the nucleus along PA = 308° and at 30″ from the nucleus on the (roughly) opposite side (PA = 130°). Pence et al. (1990) compared their posited outflow model with these extended X-ray emission regions and found that they share the same center. Radio imaging with Australia Telescope Compact Array (ATCA) at 3.5 cm (1.3  ×  0.75 arcsec synthesized beam; Morganti et al. 1999) detected the nucleus in continuum, with a potential extension in PA ∼10°, and a weak radio blob 3 arcsec to the N (PA ∼10°). The nucleus is detected by the Parkes Tidbinbilla Interferometer (PTI) at 13 cm (5 mJy; Roy et al. 1994), i.e., it hosts a compact radio source. The previous reports of outflows in NGC 1566 are consistent with a picture of a nuclear outflow driven by the AGN in which the compact base detected in blueshifted [O III] is primarily from a polar ionization cone tilted towards the observer and close to face-on and a more extended (out to 1 kpc) equatorial outflow component detected in a Hα arc in the disk of the galaxy. The reason that only the blueshifted inner (< 1″) ionization cone has been detected towards the far side of the galaxy disk most likely lies in the dusty features seen on the opposite side (NW) of the nucleus (see Fig. 8 of da Silva et al. (2017) and our structure map in Fig. 3).

C14 have presented CO J:3-2 observations of NGC 1566 using ALMA in Cycle 0: their relatively low spectral resolution (∼10.2 km s⁻¹ per channel) and relatively sparse uv coverage limited the interpretation of the molecular gas kinematics in the nuclear region. Their kinematic study of the CO emission showed a relatively regular rotational velocity field with redshifted streaming motions on the far side of the minor axis, and blueshifted streaming motions on the near side, both centered on, and within a few arcsec of the nucleus. However, they argued that the small velocity amplitudes (total width ≲ 100 km s⁻¹) of these perturbations, and the fact that they were only seen in the central 1″, makes an outflow scenario improbable, and instead attributed the non-circular rotational motions to other factors, e.g., streaming and bars. More recently Smajić et al. (2015) extended the same study by adding SINFONI observations of near-IR molecular emission lines, and came to similar conclusions, i.e., reinforcing the idea that the velocity perturbations are more easily explained by streaming motions along the minor axis as a consequence of the central bar, rather than outflows.

In this work, we reanalyze the nuclear molecular and ionized gas kinematics using new ALMA and Gemini-GMOS/IFU data. We present new ALMA observations of NGC 1566 in the CO J:2-1 emission line which covers the inner 12″ (∼600 pc) at 1.3 km s⁻¹ channel spacing, i.e., a ∼2.6 km s⁻¹spectral resolution. These new CO J:2-1 observations are more sensitive allowing us to create datacubes at the intrinsic channel spacing of the observations and have a higher image fidelity (due to the improved uv coverage from the ∼32 antenna array) as compared to the previously published Cycle 0 CO J:3-2 observations. We compare the distribution and kinematics of molecular gas with that of ionized gas (specifically the [N II] 6583 Å emission line) and stellar absorption lines observed with Gemini-GMOS/IFU at optical wavelengths. We argue that the kinematics can be best explained by a quenched spherical outflow in ionized gas, a decelerating outflow of molecular gas in the plane of the inner (∼300 pc) disk (a scenario we favor over only bar-perturbed kinematics and streaming), and discuss molecular gas streaming inflows to the nucleus.

This work is structured as follows: In Sect. 2 we present the observations and data processing. In Sect. 3 we present our results, including the morphology and kinematics of the ionized and molecular gas and stars, a comparison with our outflow, bar-perturbation and streaming models, and a discussion of the results. Finally, in Sect. 4, we present our summary and conclusions.

2. Observations, data processing, and software

We observed NGC 1566 with ALMA and Gemini-GMOS/IFU in order to obtain a comprehensive picture of the morphology and kinematics of the molecular gas, ionized gas, and stars. Images of NGC 1566, illustrating the FOVs of our observations are shown in Fig. 1.

Fig. 1. Wide field images of NGC 1566. Left panel: (ESO-DSS image from the UK Schmidt Telescope) shows the full galaxy, while the right panel (ID:13364, PI:Calzetti. HST image taken with the F555W filter) shows the inner morphology, and highlights the inner spirals arms and the ∼1′ bar in PA ∼ 0. In both panels the fields of view of the datasets used in this work are shown in yellow (GMOS/IFU), green (inner ALMA FOV: 12″ × 12″), and red (full ALMA FOV: 27″ × 27″) squares.

We observed NGC 1566 with ALMA as part of a survey of five nearby Seyfert galaxies during Cycle 2: project-ID 2012.1.00474.S (PI: Nagar) originally approved for Cycle 1 but carried over to Cycle 2. The observations of NGC 1566 were taken on June 29, 2014, using the ALMA Band 6 receivers on thirty-two 12-m antennas. Four spectral windows (SPWs) were used; two in the lower sideband (LSB) and two in the upper sideband (USB). Three of the SPWs were configured to cover the following lines at relatively high channel spacing (∼1.3 km s⁻¹): CO J:2-1 (ν_obs = 229.401922 GHz), 13CH3OH (ν_obs = 241.548041 GHz) and CS(J:5-4) (ν_obs = 243.728532 GHz). A fourth SPW was used in “continuum” mode to best detect any nuclear continuum emission. The SPWs were thus centered on 229.415 GHz, 227.060 GHz, 241.554 GHz and 243.735 GHz, with bandwidths of 1.875 GHz, 2.0 GHz, 1.875 GHz and 1.875 GHz, respectively, and spectral channel spacing of 1.27 km s⁻¹, 20.53 km s⁻¹, 1.22 km s⁻¹, and 1.21 km s⁻¹, respectively. At these frequencies, the full-width half maximum of the 12 m primary beam is about 26′. Antenna baselines ranged from 17 m to 650 m, resulting in a typical synthesized beam of $0\stackrel{″}{.}6$ × $0\stackrel{″}{.}5$ with a position angle (PA) of 25.3°.

Observations were carried out in two continuous observation blocks, totalling 124 min. The nearby radiogalaxy J0519-4546 (PICTOR A) was used as a phase, bandpass and flux-calibrator. Data were calibrated and imaged using CASA 4.2.1 (McMullin et al. 2007). The CO J:2-1 emission line was strongly detected over a velocity range of ∼ ± 200 km s⁻¹, and we were able to map the CO line at the observed channel spacing of 1.3 km s⁻¹. Thus, our effective spectral resolution (2.6 km s⁻¹) is higher than the internal dispersion of a typical GMC. At this spectral resolution, our highest spatial resolution maps (made with Brigg’s weighting with the robust parameter set to −2) have a synthesized beam of $0\stackrel{″}{.}52$ × $0\stackrel{″}{.}35$ (beam PA = 13°). The r.m.s. noise per channel in line free channels is ∼1 mJy beam⁻¹, and rises by up to a factor 2 in channels with significant line emission. Equivalent “natural weighted” maps (Brigg’s weighting with robust = 2) have a resolution of $0\stackrel{″}{.}6$ × $0\stackrel{″}{.}5$ (beam PA = 15.6°) and an r.m.s. noise per channel of 1.2 mJy beam⁻¹ in line free channels, rising by up to a factor of 4 in channels with significant line emission. The task immoment of CASA was used to create moment (integrated flux, velocity, and dispersion and skewness) maps from the above data cubes.

Gemini-GMOS observations of NGC 1566 were obtained on the night of 27th of September 2011 with GMOS in IFU mode and using the R400_G5325 grating in combination with the r_G0326 filter (program ID: GS-2011B-Q-23; P.I. Nagar). This grating yielded an intrinsic spectral resolution (FWHM) of 123 km s⁻¹, which was sampled on the CCD at ∼30 km s⁻¹ per pixel near the [N II] line. The total spectral coverage was from 5620 Å to 6970 Å. The observations consisted of two adjacent IFU fields covering 7″ × 5″ each, resulting in a total spatial coverage of 7 × 10 arcsec. Six exposures of 350 s were obtained for each field, each slightly shifted in wavelength and position in order to correct for detector defects and fill in CCD chip gaps. The data was processed using specific tasks developed for GMOS data in the gemini.gmos IRAF⁴ package.

We use four software packages for obtaining velocities, velocity fields, and related parameters from the datacubes or moment images. Ionized gas kinematics were obtained by fitting Gauss–Hermite polynomials and double Gaussians to the [N II] 6583 Å emission line using a modified version of the profit⁵ routine (Riffel 2010). The Gauss–Hermite polynomial fits were used to obtain total flux (moment 0), velocity (moment 1), and velocity dispersion (moment 2) maps over the full FOV. The nuclear stellar velocity and velocity dispersion was determined by using the Penalized Pixel Fitting (pPXF)⁶ code (Cappellari & Emsellem 2004), on the integrated (over our full FOV) spectrum of the galaxy, and using templates based on simple stellar populations (SSPs) from Bruzual & Charlot (2003).

We used a modified version of the Kinemetry⁷ package (Krajnović et al. 2006) to constrain the major axis and inclination of NGC 1566 via fits to the CO J:2-1 velocity field, and to determine the best fit circular velocity field via fits to the moment 1 (velocity) maps of both CO J:2-1 and [N II]. This modified version uses an improved global optimization thereby yielding results that are less affected by the starting values and are more robust to missing pixels in the map (discussed in detail in Ramakrishnan et al. (in prep.). Some parameters such as the position angle and inclination of the galaxy can either be fixed or obtained on the fly (see Krajnović et al. 2006, for a more detailed description of the software and its features). On the first Kinemetry run we allowed both PA and inclination to vary with radius. The PA and inclination were then fixed to their median values and a second run of Kinemetry was used to obtain the circular rotation map and the coefficients of each circular velocity and perturbation term. Briefly, Kinemetry fits concentric elliptical rings to the velocity fields under the assumption that it is possible to define the latter such that data extracted along each ellipse can be described by a simple cosine law. Therefore, along each ellipse fitted to our velocity map, the program constructs a Fourier series as a function of azimuthal angle. When using Kinemetry we used six odd terms, i.e., cos(nθ) and sin(nθ) with n = 1, 3, 5.

3. Results

The molecular (CO; ALMA) and ionized gas emission lines (Gemini-GMOS/IFU) are detected at high signal-to-noise out to the edge of the observed FOV. The [N II] emission line is detected in every pixel of the GMOS FOV at signal-to-noise ratios of 3–250 in the moment 0 maps. The CO line is detected in well defined structures which cover a fraction of the FOV: here the signal-to-noise ratio in moment 0 maps ranges between 8 and 35.

3.1. Observed moment maps: ALMA and Gemini-GMOS/IFU

Our ALMA 230 GHz continuum map shows only an unresolved nucleus and a few other weakly detected components. We do not present or discuss these 230 GHz continuum maps further since the sub-mm continuum morphology of the nuclear region can be better appreciated in the 345 GHz continuum maps of C14 due to the dust emission being brighter at this frequency. The principal use of our 230 GHz continuum map is thus to set the position of the nucleus in the CO maps. The extensive dust lanes in the nuclear region of NGC 1566 could cause a small systematic offset between the nucleus and the location of the optical continuum emission peak. Since this systematic offset is most likely to be significantly less than $0\stackrel{″}{.}5$ (see the structure map in the bottom panels of Fig. 3) we here assume that the nucleus is coincident with the stellar continuum peak in the Gemini-GMOS datacube.

The moment 0 (integrated flux), moment 1 (velocity) and moment 2 (velocity dispersion) maps of the CO J:2-1 line in NGC 1566 are shown in Fig. 2. The molecular gas in the nucleus of NGC 1566 has a clearly defined disk-like structure in the inner 3″ (3. pc), even though this region is deficient in both atomic gas and H II regions (Pence et al. 1990; Agüero et al. 2004; Smajić et al. 2015). Within this nuclear disk, the CO J:2-1 traces a two-arm spiral structure in the inner $1\stackrel{″}{.}7$ (2. pc); this spiral structure is also seen in near-infrared and optical images (Smajić et al. 2015) and in previous CO J:3-2 maps (C14). This inner molecular spiral (in which the arms almost close into a ring) sprouts two more extended but fainter CO J:2-1 spiral arms which extend out of the inner disk until roughly 3. pc ($2\stackrel{″}{.}8$). These more extended spiral arms coincide with the dust lanes seen in HST images (Fig. 3; see also C14). The CO velocity map shows velocities ranging over ± ∼140 km s⁻¹. While the disk in its inner 3″ shows a predominantly rotational “spider” velocity diagram, the velocities are asymmetrical, pointing to a warped inner disk or the presence of non-circular velocities. Despite the common association of outer H I disks to warped disks, there exists some evidence for the latter at parsec scales (Greenhill et al. 2003). Warp scenarios at nuclear scales for molecular gas has been explored as, e.g., in Schinnerer et al. (2000) reporting that molecular gas could be warped or bar perturbed in the nuclear region, although without analyzing deeply the origins of thereof. The trailing pattern of spiral arms in the CO velocity maps agrees with that at larger scales in the right panel of Fig. 1, allowing us to assume that the near and far sides of the galaxy disk is to the NW and SE, respectively. The velocity dispersion map reveals a typical dispersion of ∼30 km s⁻¹ in the inner spiral arms (see also Fig. 6), with a high (∼100 km s⁻¹) velocity dispersion region ∼3″ to the SW of the nucleus along the major axis. We note that this region does not correspond to the star-forming region which is clearly detected in the optical observations of Smajić et al. (2015). The nucleus of NGC 1566 shows a velocity dispersion of ∼60 km s⁻¹ (see also Fig. 6).

Fig. 2. Moment maps of the CO J:2-1 emission in NGC 1566. Left to right panels: 0th (integrated flux), 1st (velocity) and 2nd (velocity dispersion) moments. Top row: moment maps derived from the highest spatial resolution maps (over a 12″ × 12″ FOV) to best emphasize the nuclear features. The moment 0 (left) panel has a r.m.s. of 0.17 Jy beam−1 km s−1, and pixels with flux density lower than 0.51 Jy beam−1 km s−1 (3σ) in the moment 0 image were “masked” in all panels of the row by setting them to a value which results in a white color in the panel. Bottom row panels: moments derived from lower spatial resolution (but higher signal-to-noise) maps and show a larger 27″ × 27″ FOV to emphasize the larger scale spiral arms. The moment 0 (left panel) has a r.m.s. of 0.16 Jy beam−1 km s−1, and pixels with flux density lower than 0.48 Jy beam−1 km s−1 (3σ) in the moment 0 images were “masked” in all panels.

The moment maps of the [N II] line in NGC 1566, obtained from the GEMINI-GMOS/IFU data, are shown in the top row of Fig. 3. Overall, these are roughly similar to those of the CO line. The [N II] moment 0 image clearly shows the bright [N II] region to the SW, corresponding to the optically-emitting star-forming region seen in Smajić et al. (2015) and in agreement with a blue region seen to the SW of the nucleus in the CO velocity residual map. The velocity map shows kinematics consistent with rotation with velocities similar to those seen in the CO maps. Once more there are non-symmetrical velocity patterns closest to the nucleus: note especially the excess of blueshifts seen on the far side of the galaxy disk ∼2″ from the nucleus. The map of the [N II] velocity dispersion is more difficult to interpret. The [N II] dispersion is in general higher than that seen for CO, and the inner spiral structure is not as clearly discernible as a higher dispersion region. The star-forming region to the SW has a dispersion of ∼50 km s⁻¹ in [N II] (less than that in CO). The nucleus shows a velocity dispersion (∼120 km s⁻¹) significantly higher than that seen in CO and two regions ∼2″ from the nucleus in the NE and S directions also show relatively high (150–180 km s⁻¹) velocity dispersions. These will be interpreted below in conjunction with the results of the two component fits to the [N II] emission line.

Fig. 3. Top row: moment maps of the [N II] 6583 Å emission line in NGC 1566. From left to right, the panels show the 0th (integrated flux), 1st (velocity) and 2nd (velocity dispersion) moments. The nuclear position (marked with a cross) was derived from the peak flux of the continuum in Gemini-GMOS datacube, and the solid lines indicate the major axis (thick line) and the minor axis (thin line). In the moment 0 map (left panel), the brown contours show the integrated flux (moment 0) of the CO J:2-1 line. Bottom row: from left to right the structure map, and maps of the stellar velocity and stellar velocity dispersion. All maps have N to the top and E to the left (see compass). The structure map was created from a HST F606W filter on which unsharp masking was used to emphasize dust features. It is shown at the same size and orientation as the other panels.

The bottom row of Fig. 3 shows a structure map and the first two moments of the stellar velocity field. The structure map was created by running the IDL routine “unsharp_mask.pro” on a HST image taken through the F606W filter, in order to emphasize sharp changes in the image. While the highest contrast dust arc is seen on the far side of the galaxy as expected, several strong dust features are also visible on the near side of the galaxy. The stellar velocity map, derived from running pPXF on a Voronoi binned datacube (to achieve a minimum signal-to-noise of 25 in the continuum near the [O III] line in each spectrum), shows a clear rotation pattern. The stellar rotation velocities are significantly lower than those seen in the molecular gas. Since the map is relatively noisy even after Voronoi binning we did not attempt to fit a PA and inclination to this velocity field using Kinemetry. Visually, the PA appears consistent with the values we derive from our CO J:2-1 map; this is corroborated by our best fit rotation model to the stellar velocity field (see next section).

3.2. Modeling observed velocities: rotation

We used Kinemetry to analyze the CO J:2-1 velocity field (see Sect. 2), both to constrain the PA and inclination of the CO disk, and to constrain the relative contributions of circular rotation and perturbations. We assumed m = 2 modes and thus use six (odd) Fourier decomposition terms. In the first run both the PA and inclination were allowed to vary with radius, and in the second run we fixed both to their median values from the first run. The Fourier decomposition coefficients of the best fit Kinemetry model are shown in Fig. 5. Here the cosθ term represents the pure circular (rotation) velocity and the other terms are perturbations. In the innermost ∼1″, the sinθ (radial) term is positive and dominates the pure rotation term (below we argue that this is best explained by a nuclear outflow) while the cosθ (circular rotation) term dominates between ∼1 and 4″. Beyond $3\stackrel{″}{.}4$, the CO velocity field is sparse and the results of Kinemetry are thus less reliable. Nevertheless we note that the coefficient of sinθ remains stable at ∼20 km s⁻¹at radii beyond $1\stackrel{″}{.}4$, and the cos3θ coefficient is significantly negative between 4″ and 5″ which could signify an asymmetry about the minor axis. The other terms show relatively small amplitudes, and given that we sample a very small range of radii (significantly less than the bar co-rotation radius) we are unable to reliably interpret their variations. To emphasize the changing reliability of these results with radius we plot two dashed vertical lines. To the left of these is the region with a 100% of data coverage in the ellipse fitting (between 0 and $2\stackrel{″}{.}2$); between these there is a linear decrease from 100% to a 50% and to the right of these the data coverage of the ellipse fitting is lower than 50%, decreasing to 30% at $3\stackrel{″}{.}4$ and beyond. It is thus clear that beyond 3 to 4″ the values of the coefficients are relatively unreliable due to the sparseness of the velocity field, and we thus do not attempt to interpret, e.g., the fact that the s1 term remains positive and almost constant.

Fig. 5. Amplitudes of the Fourier components obtained from the Kinemetry analysis of the CO velocity map, as a function of distance from the nucleus. The solid red line represents the coefficient of the cosθ term (pure circular rotation in a disk) and the other solid lines show the additional components (associated to perturbations) in the Fourier decomposition, following the colors specified in the inset. Only odd Fourier components were used. The vertical lines delineate radii at which we have abrupt changes in the fraction of pixels at a given radius which have values in the velocity map. At radii smaller than the left vertical line, this fraction is 1. Between the two vertical lines the fraction drops linearly from 1 to 0.5, and beyond 4″ the fraction is relatively steady at ∼0.3. The horizontal dashed line delineates zero velocity.

We fit the observed stellar velocity field with a “Bertola” model rotation curve (Eq. (2); Bertola et al. 1991). This model uses six parameters: the maximum amplitude of the rotation curve A; the radius at which this maximum amplitude is achieved c; a p factor which drives the slope of the rotation curve at larger radii (p = 1 gives a flat rotation curve at large radii and p >  1 gives a decreasing rotation curve at large radii, emulating a finite total mass in the disk), the position angle of the major axis; the inclination of the galaxy, and the systemic velocity. Of the six parameters of the model we fixed the c parameter (by visual inspection of the radius at which the rotation amplitude reached its maximum) and the inclination (the same value obtained running Kinemetry on the CO velocity field). The best fit parameters obtained were A = 200 km s⁻¹, p = 1.5, PA = 45° (same as derived by Kinemetry on the CO velocity field), and a systemic velocity offset of −5 km s⁻¹ with respect to that used in this work (1485 km s⁻¹). This rotation model is shown with a dash-dotted line in Fig. 4, and is later (Sect. 3.6) used in the bar perturbation analysis.

Fig. 4. Comparison of the (projected) rotation models and observed velocities extracted along the major axis, in the central kpc of NGC 1566. The solid body model is shown with a dotted line, our ModC2014 model (see text) is shown with a dashed line, the circular (“gascirc”) model obtained by the Kinemetry fit to the CO J:2-1 velocity field is shown with the solid line, and the Bertola model fit to the stellar velocity field is shown with a dashed-dot line. CO J:2-1 and stellar velocities extracted along the major axis are shown with crosses and stars, respectively. Blue is used for the NE (approaching) side of the galaxy, and red for the SW (receding) side of the galaxy. The zero velocity on the y-axis corresponds to 1485 km s−1; at the nucleus, the velocity of the CO J:2-1 line is offset from this by −12.8 km s−1; a consequence of the asymmetric double-peaked profile of the nucleus, most likely caused by the effects of nuclear outflows and/or bar related perturbations.

The asymmetries in the observed CO and [N II] velocity fields are best appreciated once ordered rotational motions are subtracted out. We use three rotation models for our analysis of the observed gas velocity maps: all models use a major axis P.A. of 45° and an inclination of i = 33°, derived from running Kinemetry on our CO J:2-1 velocity field. These three models differ in the parametric form of the circular rotation velocity with radius: (a) solid body rotation with parameters:

$$\begin{array}{c}\hfill {\mathrm{Vc}}_{\mathrm{radial}}=({\mathrm{S}}_{\mathrm{rot}}\ast r)\ast \cos (\varphi )\ast \sin (i).\end{array}$$(1)

where S_rot is equal to 312. km s⁻¹pc⁻¹, r and ϕ are the polar coordinates in the velocity map, and i is the galaxy inclination; (b) An empirical axisymmetric rotation model (hereafter, ModC2014) based on the nuclear rotation curve derived by C14, but with a gradual decrease in circular velocity beyond 2″, roughly following the results from the Kinemetry fit to the CO J:2-1 velocity field. Recall that C14 derived the nuclear rotation curve by using their observed CO J:3-2 velocity field - specifically by minimizing the residual (observed − model) velocities - at small radii, and literature Hα velocities at larger radii (the black line in their Fig. 9). Since the nuclear CO kinematics are highly perturbed and the velocity field is relatively sparse (especially beyond $2\stackrel{″}{.}5$), it is not clear that a Kinemetry fit or a minimization of residuals will produce a reliable circular rotation model. In fact the circular velocity model obtained by applying Kinemetry to our CO velocity field, and that from fitting a Bertola model to the stellar velocity field, are both significantly different from the C14 rotation model (Fig. 4): the most significant difference is a decrease in the rotation velocities beyond 2″. Since a major function of the rotation model in the following sections is to emphasize asymmetries in the observed velocity field, we create a smooth rotation model which follows C14 (and the solid body model above) in the inner 2″ and decreases (to reflect the Kinemetry-derived and Bertola models) at larger distances (dashed lines in Fig. 4). This ModC2014 model is used to produce residual (observed − model) velocity maps which better emphasize deviations from circular rotation as compared to using the C14 rotation model; and (c) the gas circular rotation model (specifically the variable “gascirc”) obtained by running Kinemetry on the CO J:2-1 velocity map.

A direct comparison of the rotation models considered by us in the inner kpc of NGC 1566 is shown in Fig. 4. The CO J:2-1 velocities extracted along the major axis (crosses) show several differences from the solid body and ModC2014 models. First, the velocities are not axisymmetric with the blueshifted velocities (to the NE) higher than the redshifted ones (SW) between 1–2.5 arcsec from the nucleus. Second, both blue- and redshifted sides show wiggles with higher velocities in the inner $0\stackrel{″}{.}5$ and relatively low velocities at distances ≥$2\stackrel{″}{.}5$ arcsec from the nucleus. The circular gas velocity fit of Kinemetry to the CO J:2-1 velocity field (solid line) is by definition an axisymmetric fit; it well follows the solid body rotation model until $2\stackrel{″}{.}4$ (except for an excess of velocities in the inner $0\stackrel{″}{.}5$ (below we argue that this is due to a nuclear outflow) after which it shows rotation velocities slightly lower than ModC2014.

On the other hand, the stellar velocities along the major axis (stars in Fig. 4) presents two stages: the blue- and redshifted sides are quite similar until ∼2″: in the inner arcsec it is completely consistent with our ModC2014 model. Beyond 2″, on the NE (blueshifted) side the stellar velocities flatten at ∼60 km s⁻¹ while on the SW (redshifted) side the stellar velocities reach a similar peak but then decrease to almost zero velocity at $2\stackrel{″}{.}8$, after which they gradually increase again to ∼40 km s⁻¹.

For the three models – solid body rotation, ModC2014, and Kinemetry-derived – the residual (observed minus model) CO J:2-1 velocity maps reveal similar asymmetries in the inner $2\stackrel{″}{.}6$. The differences in the residual maps are seen at greater radii: the Kinemetry model undersubtracts the observed velocities while the ModC2014 and solid body rotation models oversubtract the observed velocities. The Kinemetry model results in the best and most symmetrical residual velocity map for CO J:2-1, but does not work well for the GMOS/IFU data. For GMOS-derived ionized gas velocity fields, the inner region is well subtracted but a radii greater than 2″ the velocities are not well subtracted as a consequence of the very low model velocities; i.e., rotation velocities in ionized and stellar lines do not decrease at radii greater $2\stackrel{″}{.}6$ as in the case of CO J:2-1. Given that we are interested primarily in identifying deviations from axis-symmetric rotation in the innermost region, rather than accurately predicting the true rotation curve, unless otherwise stated, we consistently use the ModC2014 model for all (ALMA and GMOS) kinematical analysis in this work.

The residual (after subtraction of the rotation model) velocity field of the CO J:2-1 emission line is shown in Figs. 6 and 7. The departures from pure rotation are now clearer, especially in the inner 3″. The largest deviations are (a) blue- and redshifted clumps ∼1″ to the NW and SE of the nucleus; (b) red spiral arms ∼4″ to the N and SE with the latter less redshifted, and (c) a diffuse clump some redshifted ∼3″ to the SW of the nucleus. The blue (residual velocity ∼50 km s⁻¹) clump ∼$1\stackrel{″}{.}5$ to the SW of the nucleus along the major axis in our CO residual velocity map (Fig. 7) marks the location of the star-forming region noted by Smajić et al. (2015): recall that this region is discernible in our [N II] moment maps, both as a high-flux region in the [N II] moment 0 map, and a relatively low velocity dispersion region in the [N II] moment 2 map. Further, this is also the region found to have a velocity gradient in its [O III] emission line da Silva et al. (2017). Disturbances produced by the star-forming region could explain why the CO velocity here does not agree with the expectation from pure rotational.

Fig. 6. Top panel: maps of the CO J:2-1 velocity dispersion (left panel) and velocity residuals after subtraction of the ModC2014 rotation model (right panel) for the inner FOV. Contours in both panels show the moment 0 (integrated flux) map of the CO J:2-1 emission, ranging from 0.4 mJy beam−1 km s−1 to 5.4 mJy beam−1 km s−1. Bottom left panel: moment 3 (skewness) map of the CO J:2-1 emission for the larger FOV, shown in color following the color bar: blue colors represent spectra with excess emission towards the blue side of the weighted mean velocity. The CO Moment 0 (total flux) map is overlaid with black contours. Bottom right panel: as in the top right panel but for the larger FOV. Pixels with velocity less than −90 km s−1 in the right panels are shown in white. Contours in bottom panels show the moment 0 (integrated flux) map of the CO J:2-1 emission, ranging from 0.4 mJy beam−1 km s−1 to 8.4 mJy beam−1 km s−1.

Fig. 7. Nuclear velocity field of NGC 1566 in the CO J:2-1 (top row) and [N II] 6583 Å (bottom row) emission lines, in contours and color following the color bar above each panel. The left column shows the expected radial velocity field from the ModC2014 model (see text), the middle column shows the observed radial velocity, and the right column shows the residual velocities (observed − model). A cyan dashed squares were drawing inside bottom panels for showing the FOV of top panels. In these figures we use a systemic velocity of 1485 km s−1and the nuclear position (marked with a cross) is determined by the 230 GHz continuum emission peak. The major (PA = 45°) and minor axes are shown in solid lines. One arcsecond corresponds to 0. pc.

3.3. Modeling observed velocities: molecular outflows?

The CO J:2-1 residual velocity map in Fig. 7 shows excess blueshifts to the NW (near side) and redshifts to the SE (far side) of the galaxy, around 1″ from the nucleus along the minor axis. This is the expected signature of outflows within the plane of the galaxy disk. This feature (but at lower spatial resolution, and at significantly lower spectral resolution and image fidelity) was noted by C14 and Smajić et al. (2015), but they argued that the low velocities seen in their CO J:3-2 maps made an outflow scenario unlikely. Smajić et al. (2015) have shown that residuals in the nuclear kinematics of the H₂ are consistent with outflows along the minor axis, but that these residuals are also easily explained by deviations introduced by the density waves of the nuclear spiral; non-circular orbits, e.g., a closed elliptical orbit with axes not parallel to one of the symmetry axes (minor or major) can produce residual velocities (Smajić et al. 2015).

We nevertheless argue for the presence of a nuclear outflow (which may of course co-exist with other bar-, warp, or spiral-related perturbations) based on the following reasons: (a) the presence of a nuclear outflow in NGC 1566 is not unexpected as previous studies have claimed kinematic and morphological evidence for the presence of outflows in the NLR (Schmitt & Kinney 1996; Davies et al. 2016; da Silva et al. 2017), which are most likely to intersect the disk given the observed geometries, and in the larger scale disk (Pence et al. 1990). Note also the evidence of a blueshifted velocity in [O III]5007 Å near the star-forming region ∼1$1\stackrel{″}{.}5$ SW of the nucleus discussed above which da Silva et al. (2017), see their Fig. 20) interpreted as a consequence of contamination from an AGN outflow; we also detect this blue residual in our CO residual (observed − rotation model) velocity map (Fig. 7); (b) the unresolved nuclear aperture shows a double peaked profile with Full Width at Zero Intensity (FWZI) ∼200 km s⁻¹ (Figs. 13 and 15), higher than that seen in the lower fidelity maps of C14. If these velocities are attributed to an outflow, the fact that opposite velocities are observed on each side of the nucleus implies that the outflow axis is not aligned with our line of sight. Large angles to the line of sight are unlikely as this would imply extremely high true outflow velocities. An outflow in the plane of the disk would imply an outflow with velocities up to 180 km s⁻¹. Conversely, attributing these observed velocities along the minor axis to other perturbations in the plane of the disk, requires radial velocity perturbations of ∼80–100 km s⁻¹ in a nuclear region where the intrinsic (undisturbed) rotation velocities are expected to be ≲40 km s⁻¹. As we will show in Sect. 3.6, our modeling of the bar-related perturbations does not reliably produce both the morphology and the large perturbations seen in the observed velocity field; (c) the pv diagram along the minor axis (bottom right panel of Fig. 9) not only shows the high-velocity components (±60–90 km s⁻¹) in the nuclear aperture but also lower brightness emission which connects these high-velocity components to the zero velocity components seen at r∼$1\stackrel{″}{.}8$ on both sides of the nucleus. To the N (negative offsets in the bottom right panel of Fig. 9) the decrease in velocity in the inner arcsec is clearly seen, and to the S (positive offsets) the decrease is more clearly seen in the r ∼ 1 − 2″ range. (d) the pv diagrams show velocity deviations which are consistent with radial outflows in the plane of the disk over several PAs (and not just the minor axis) and over apertures at distances of several synthesized beams from the nucleus (see Sect. 3.5); (e) as discussed in the next section, the kinematics of the emission line gas in the optical (from GMOS-IFU) are consistent with a nuclear spherical (or bipolar) outflow, which makes the interpretation of a related molecular outflow less surprising. In summary, we support the presence of a nuclear outflow, without ruling out the presence of additional bar- or spiral-related perturbations (Sect. 3.6). Other scenarios, e.g., a warped disk or non-coplanar disks (e.g., Wong et al. 2004; García-Burillo et al. 2014) cannot be constrained by us due to the limited resolution, the sparse velocity field, and the lack of a reliable circular rotation model for the galaxy.

This posited molecular outflow is most likely primarily in the disk of the galaxy for the following reasons: (a) outflows outside the disk are often related to nuclear jets (e.g., Morganti et al. 2013; Sakamoto et al. 2014), but there is no clear evidence for radio-traced jets and outflows in the nucleus of NGC 1566. The potential radio extension in PA 10° (Sect. 1) is not aligned with our posited outflow axis. Further, apart from the Seyfert 1 classification there are no data (e.g., maser disks) to constrain the orientation of the central engine; (b) outflows are also often seen perpendicular to the plane of the disk, especially in the case of starburst driven winds (eg., Veilleux & Rupke 2002; Veilleux et al. 2005; Leroy et al. 2015). However, for such a polar outflow the blueshifted (redshifted) emission would be seen towards the far (near) side of the galaxy disk, the opposite of that seen in NGC 1566; (c) a large (e.g., > 30°) opening angle for the outflow, often the case in radiation-pressure-, jet-, and starburst-driven outflows will produce a high observed velocity dispersion due to the varying projection angles of the outflowing gas to the line of sight. In NGC 1566 we see a relatively low CO J:2-1 velocity dispersion (FWHM of ∼30 km s⁻¹) in the posited nuclear outflow components (see Figs. 6 and 15). A spherical outflow scenario can be rejected as this would produce a high velocity dispersion centered on zero velocity, under the assumption of optically-thin emission in CO J:2-1. It is important to note that in Sect. 3.4, we argue for the presence of spherical outflow in the ionized gas, which does not contradict our claim that the molecular outflow is in the disk; (d) the posited nuclear outflow has a limited extension, and an apparently decreasing velocity, both of which argue for deceleration of the molecular gas in the high-density medium of the disk; (e) the higher molecular gas density in the disk of the galaxy will make this component more easily detectable in short integrations, as compared to more diffuse molecular gas outside the galaxy disk. In summary, while the molecular outflow could have a larger opening angle (and indeed be isotropic) we appear to be preferentially detecting this component within the galaxy disk.

Is the posited outflow AGN or starburst-driven? As mentioned above, there is no clear evidence of a radio jet, and so any AGN-driven outflow would most likely be due to radiation-pressure. We note that the highest velocities in the outflow are detected at the position of the nucleus and not towards the star-forming knot $1\stackrel{″}{.}5$ to the SW. Several authors have presented photometric and kinematic data which argue against the likelihood of a starburst-driven outflow in NGC 1566: Davies et al. (2016) used diagnostic diagrams of [O III]/Hβ to demonstrate that there is no significant contribution from star-forming regions in the nucleus and it is the radiation pressure from the AGN which dominates in the inner scales. Smajić et al. (2015) present similar diagnostic diagrams using Molecular Hydrogen (H₂) and report an AGN domination for the nuclear region inside 1″ and a relatively modest SFR (∼8 × 10⁻³ M_⊙ yr⁻¹) in the inner 3″, implying a relatively low star-forming efficiency given the observed gas reservoir. Both the galaxy-wide SFR (∼4.32 M_⊙ yr⁻¹; Gruppioni et al. 2016) and the SFR surface density (∼0.033 M_⊙ yr⁻¹ kpc⁻²; Hollyhead et al. 2016) are relatively low, and thus star formation is not expected to drive a nuclear outflow (e.g., Cicone et al. 2014).

As seen in the pv diagrams of Fig. 9, a model which sums the ModC2014 model and our empirically derived outflow model (black solid lines in the pv diagrams; see Sect. 3.5) provides a much better fit (as compared to a pure rotational model) to the inner 2″ in the pv diagrams at all PAs. It should be noted that our synthesized beam of ≲0.5 arcsec well resolves the central 4 arcsec of the galaxy (e.g., Fig. 7). Nevertheless there are several specific features which cannot be fit only with the model of radial outflows (in the galaxy disk) plus rotation, e.g., the apparent morphological double structure of each outflow lobe (Fig. 7) and the pv diagram in PA = 75° (middle right panel of Fig. 9) where observed velocities 1.5–2 arcsec from the nucleus to the E are not well fit by the model. These are discussed in Sect. 3.5.

In the residual velocity map of CO (top right panel of Fig. 7) the blue- and redshifted lobes 1″ from the nucleus have a double-peaked morphology. The largest velocity deviations are along PAs of 100° and 140°, i.e., straddling the minor axis. How can this be explained? Do streaming motions into the nucleus along PA ∼120° create the valley between the two peaks? Further, why is it that almost the entire near side of the galaxy (NE) has a blue residual, while almost the entire far side of the galaxy (SW) has a red residual in the top right panel of Fig. 7? This is not due to the use of an incorrect major axis PA (varying the major axis PA does not change these features). Effectively, the NE and SW sides of the galaxy are not axisymmetric in their rotation (see, e.g., Fig. 4).

3.4. Modeling observed velocities: ionized gas outflows?

The [N II] line velocity map was (obtained from the Gauss–Hermite fit version of profit.pro after subtracting the broad Hα line emission). The [N II] velocity residual map constructed using the same rotation model used for the CO J:2-1 line at first glance appears morphologically similar to its CO J:2-1 counterpart (bottom panels in Fig. 7). However the P.A. of the major axis for the [N II] rotation (especially to the SW) appears to be ∼15° smaller than that used in our model. Further, the [N II] residual velocity map clearly shows blue residual velocities at ∼1–2″ on both sides of the nucleus (with the largest velocity deviation to the SE), in contrast to the CO J:2-1 residual velocity map which clearly shows opposite colors on each side of the nucleus (interpreted above as the sign of a nuclear molecular outflow in the disk).

In the [N II] velocity residual map, the blue region on the near side of the galaxy (NW of the nucleus) is in rough agreement with the equivalent blue region in the CO J:2-1 residual map, and thus would be consistent with the outflow scenario posited for the CO J:2-1 data. On the other hand, the blue region in the [N II] velocity residual map on the far side of the galaxy (SE of the nucleus) is located along the minor axis at ∼$1\stackrel{″}{.}8$ from the nucleus: this does not match the feature seen in the CO J:2-1 residual map which is closer to the nucleus (∼1″) and redshifted. It should be noted that this blue SE feature in the [N II] residual velocity map is roughly cospatial with a region of high dust extinction (structure map in Fig. 3). A red region closer to the nucleus matches with the red region to the SE seen in the CO residual map, however, this consists in redshifted radial velocities around 10 km s⁻¹ which is too low to be considered significant as a part of an outflowing gas. There is another small region to the SW at ∼2″ which is redshifted and around 50 km s⁻¹ but in spite of having a similar dynamic range, it does not match with the redshifted outflow region to the SE in the ALMA residual map (∼$1\stackrel{″}{.}5$), being placed too close to the star-forming region.

The [N II] emission in the inner ∼2″ is double-peaked, and thus cannot be well fit with a single Gauss–Hermite polynomial. We, therefore, used the two Gaussian fit version of profit.pro to search for potential independent velocity components in the [N II] emission line. Meaningful two-component fits were obtained in part of the inner 3″ radii nuclear region and the resulting velocity maps are shown in Fig. 8. These maps include only pixels for which a double component fit produced a meaningful result. For regions where the two-component fit was not possible, the single component fit remains valid. It should be noted that in the latter case the single component fit shows primarily the equivalent of the narrow component but in a few regions the velocity dispersion of the single component fit is similar to that of the broad component of the two-component fit. The two-component fit is mainly obtained near the nucleus and on the near side of the galaxy. On the far side (and a small nuclear region to the W on the near side), the regions in which a two component fit is not possible coincide well with the regions of high dust extinction. The “narrow component” shows velocity dispersions ranging between 60 and 90 km s⁻¹ (bottom left panel of the figure) and the velocity field of this component mainly follows that expected from our rotation model. The second component, which we refer to as the “broad component”, shows velocity dispersions of ∼140–160 km s⁻¹ (bottom right panel of the figure) and does not appear to participate in regular rotation.

Fig. 8. Velocity and dispersion maps of the [N II] emission line obtained from the two-Gaussian fit version of profit.pro. Top panels: from left to right, the velocity and velocity dispersion maps of the narrow and broad components, respectively. Bottom panels: velocity dispersion maps of the narrow (left panel) and broad (right panel) components.

The strong blue residual region to the SE in the single component [N II] residual velocity map is not fit with a double component. The weaker blue residual region to the NW in the single component [N II] residual velocity map is now seen to be blue in its broad component only; in its narrow component this region follows regular rotation. It should be noted that both blue knots do not correspond to regions of high dust columns in the structure map.

Based on the velocity field of the broad component of the two component fit to [N II], the velocity field of the high dispersion (≳120 km s⁻¹) areas in the single component fit to [N II] and the distribution of the nuclear dust, we postulate the presence of an expanding sphere of ionized gas, i.e., a spherical ionized outflow, for the reasons given below. With the presence of dust (dominantly in the plane of the galaxy disk) we would preferably see emission from the hemisphere in front of the galaxy and moving towards us, i.e., blueshifted radial velocities. In the absence of dust one would expect a high dispersion and a median velocity close to systemic. It should be noted that given the relatively low inclination of NGC 1566 (33°), dust in the inner 2″ (∼100 pc) of the disk produce an almost equal extinction of light from the bulge for both the “near” and “far” side of the galaxy disk. This is clear in the structure map where dark dust lanes are seen on both the near and far sides in the inner ∼2″. Only at larger radii are the structures of the dust lanes more prominent on the near side of the galaxy disk. In the case of NGC 1566 we do not obtain a two-component fit in areas where the structure map implies marked dust lanes (non-intuitively these are on the far side of the galaxy disk in the innermost arcsec) and find large blueshifts in the broad component of the two-component fit on the near side of the galaxy disk in areas where the structure map shows less marked dust lanes. Further, the maps derived from the single component fit to the [N II] line (Fig. 3) show two regions of high dispersion, about 2″ from the nucleus to the SE (the blue knot referred to above) and to the NE. These two regions effectively correspond to the broad component, and are also blue in their velocity. Thus we effectively see blue velocities in the broad component in almost every compass direction implying a spherical outflow in the inner ∼2″ which is visible to us primarily in areas of lower dust extinction. Similar kinematic signatures, and thus interpretations, were observed and used in previous IFU studies, e.g., in the nearby Seyfert galaxy NGC 2110 (Schnorr-Müller et al. 2014b), and in some radio-quiet quasars (Liu et al. 2013). Nevertheless, despite all these signs, we are not neglecting a potential presence of a bipolar outflow. Namely, it might be reasonable to think that what we see as a spherical outflow might be instead ionized gas ejected from the nucleus in opposite directions but we notice just a part of them as a consequence of a poor resolution in the GMOS data (∼0.5–$0\stackrel{″}{.}6$).

We note that our postulation of this spherical outflow of ionized gas does not contradict our postulation of cold molecular gas outflows in the galaxy disk. The cold molecular gas outflows are expected to be preferentially detected in the plane of the disk for two reasons: a low-density CO outflow would be optically thin and its profile would thus be centered on zero velocity and the significantly higher abundance of molecular gas in the disk as opposed to above the disk, makes the disk molecular gas much easier to detect.

3.5. Observed position-velocity (pv) diagrams: ALMA

In this section we present position-velocity (pv) diagrams for the CO J:2-1 line. These pv diagrams were extracted from the ALMA datacubes using a slitwidth of $0\stackrel{″}{.}2$; they are thus limited in spatial resolution by the intrinsic resolution of the images (∼0$0\stackrel{″}{.}4$). Fig. 9 shows the pv diagrams along the major axis, minor axis, the large-scale bar, and other relevant PAs including PA = 115° (20° from the minor axis) for a direct comparison with the pv diagram of CO and Hα in Fig. 3 of Agüero et al. (2004). In their figure, the “Blueshifted Knot” which they interpret as inflow motion is clearly seen in our CO J:2-1 pv diagram.

Fig. 9. Position-velocity diagrams of the CO J:2-1 emission in NGC 1566 along several PAs are shown in color following the color bar above each panel. The PA of the “slit” over which the pv diagram was extracted is indicated above each panel, as are the PAs corresponding to the major and minor axis of the galaxy, and of the large-scale bar. Negative offsets on the x-axis correspond to the PA listed above the panel, i.e., positive offsets are along the 180° plus the listed PA. The black cross indicates the position of the 230 GHz continuum peak (presumed to be the galaxy center) and the systemic velocity of the CO J:2-1 line (1485 km s−1). To guide the eye horizontal and vertical lines delineate ±2 arcsec from the nucleus and ±50 km s−1 from the systemic velocity. The dashed white and purple lines are the solid body rotation model and the ModC2014 rotation model, respectively. The solid black line shows the expectation of adding our outflow model to the ModC2014 model (see Sect. 3.3). The pv diagrams, were created from a “hybrid” datacube: the inner ∼12″ × 12″ square region centered on the nucleus is taken from a high-resolution (Briggs weighting, Robust = −2) map with spatial resolution $0\stackrel{″}{.}52$ × $0\stackrel{″}{.}35$ and an r.m.s. noise of 0.1 mJy beam−1 per channel (up to 2 mJy beam−1 per channel in channels with strong signal), while the rest of the cube is from a higher signal-to-noise (Briggs weighting, Robust = 2) map with spatial resolution $0\stackrel{″}{.}6$ × $0\stackrel{″}{.}5$ and an r.m.s. noise of 1.2 mJy beam−1 per channel (up to 4 mJy beam−1 per channel in channels with strong signal).

In each pv diagram we have overlaid the predictions of solid body rotation (white dashed line), the ModC2014 model (purple dashed line), and the sum of the ModC2014 model with our outflow model (solid black line). The solid body rotation and ModC2014 are essentially the same over the inner ±3 arcsec, after which the latter flattens in velocity. Our outflow model was derived as follows: we used the pv diagram along the minor axis (bottom right panel of Fig. 9) to measure the radial velocity of the brightest CO emission at a given distance from the nucleus on both the NW (positive offsets in the pv diagram) and SE sides. The absolute values of these velocities as a function of distance from the nucleus were then interpolated and deprojected (assuming that the outflow is in the disk) to construct a function of outflow velocity vs. position. As seen in the bottom right panel of Fig. 9 the redshifted “outflow” velocities provided more constraints closest to the nucleus and the blueshifted velocities provided better constraints at slightly larger distances. The final outflow model starts with outflow velocities in the disk of 157 km s⁻¹ at the (unresolved) nucleus and decreases monotonically to zero velocity 2″ from the nucleus.

First concentrating on the pv diagram along the minor axis in Fig. 9, we see that the NW side of the minor axis shows a clear deceleration in velocities when going from 1″ to 2″ from the nucleus, while the SE side the bright emission at 1.5–2″ (the inner spiral arm) shows a higher velocity dispersion which does not clearly vary with distance. However, this SE side shows a clearer decrease in velocities between 0 and 1″. A similar scenario has been reported in NGC 1068 for both hot and cold molecular gas (Barbosa et al. 2014; García-Burillo et al. 2014, respectively); outflowing nuclear molecular gas, with outflow velocity decelerating from 200 km s⁻¹ to 0, accumulates in an off-centered ring 100 pc from the nucleus (Barbosa et al. 2014).

For the pv diagram along PA = 0, which is aligned with the large-scale bar in NGC 1566, we see (1) gas consistent with outflows, which is well fit by our outflow model; (2) gas which is almost in rotation 1–2″ from the nucleus on either side, but showing a steeper rotation curve which reaches zero velocity at a position offset from the center.

For the pv diagram along the major axis (PA 45°), several velocity components can be seen. These include gas in rotation, and some contamination from the disk outflow component, since at the nucleus, the “slit” (limited by the spatial resolution of our ALMA observations) expectedly picks up the gas outflowing along the minor axis and other angles. Moreover, gas in the inner spiral to the NE is preferentially redshifted and gas in the inner spiral to the SW is preferentially blueshifted, i.e., both spirals show ∼40 km s⁻¹ (in projection) deviations towards values of zero velocity: the most obvious interpretation of this is that gas originally in circular rotation is slowed down on hitting the ends of the nuclear molecular gas “bar”. This loss of momentum could potentially result in inflows. Similar velocity offsets are also seen at slit PAs offset 10° from the major axis (e.g., PA = 55°; Fig. 9), and is very dramatic on the SW side in the pv diagram with slit PA = 75°, at the point where an outer spiral pattern breaks off from the inner spiral pattern.

Molecular gas in the inner spirals always show a high velocity dispersion (around ∼80 km s⁻¹). The inner spiral to the NW (about 1 arcsec from the nucleus) always show velocities which are bluer than that expected from rotation or rotation+outflow. This is clearly seen in all pv diagrams which intersect this arm (e.g., PAs 45°, 55°, 75°). The opposite inner spiral (that to the SE) shows the opposite, i.e., velocities redder than expected from rotation and rotation+outflow (e.g., most obvious in the PA 105° and 115° pv diagrams) but this velocity offset is not as well defined as in the case of the NW arm. Given that the NW arm is mainly on the near side of the galaxy disk and the E arm is on the far side, this is what would be expected from a streaming outflows along the spiral rather than streaming inflows! We speculate that these structures are absorbing the momentum of the nuclear outflow and thus heating up and expanding.

Another region which consistently shows large differences from the rotation+outflow model is the double cavity (between the nuclear bar-like structure and the inner spiral arms) to the NW and SE, for offsets of 1–2 arcsec from the nucleus. When the slit passes through these cavities the pv diagrams (especially those at PA = 0°, and PA = −15°) show a characteristic pattern which can be explained by slower than rotation velocity between 1 and 2 arcsec (but increasing in the correct sense), and a high dispersion of velocities (all higher than those expected from rotation) at 2 arcsec.

While the predictions of our rotation plus outflow model are in general consistent with the position velocity data at different PAs especially in the inner 2″, at PAs close to 0° (see the pv diagrams at PA = 0°,−15° and 15°) one can see significant differences between the data and models in the inner $1\stackrel{″}{.}5$. Here we clearly see a component of gas which follows a steep velocity gradient decreasing to zero velocity at a distance of 1$1\stackrel{″}{.}5$ from the nucleus on both sides. One potential explanation for this anomalous rotation is an inner counter rotating gas disk with major axis in PA ∼ 0°, fed by gas inflowing along the large-scale bar. This possibility is motivated by observational evidence that bars are an efficient pathway for transporting gas from galactic scales to nuclear scales in both active and inactive barred galaxies (Sakamoto et al. 1999; Crenshaw et al. 2003; Regan & Teuben 2004; Sheth et al. 2005). Alternatively, these features are a consequence of perturbations due to the bar, as discussed in the next section.

3.6. Modeling observed velocities: bar perturbations

Velocity perturbations due to bar(s) are believed to play an important role in fueling the SMBH and in triggering nuclear star formation. It is well known that NGC 1566 hosts an intermediate strength nuclear bar with radius ∼1.7 kpc in PA ≈ 0° (Hackwell & Schweizer 1983; Pence et al. 1990; Mulchaey et al. 1997; Agüero et al. 2004; Dicaire et al. 2008; Comerón et al. 2010; Kendall et al. 2011, C14), which could be largely responsible for the velocity perturbations seen in the molecular gas. C14 used torque maps to show that the asymmetries in the velocity field of the nuclear molecular gas are predominantly produced by the bar. They also briefly explore estimates for the bar pattern speed. However, they did not make a detailed kinematic analysis of the bar-produced perturbations.

To analyze the effect of bar-produced perturbations we use both Diskfit (Spekkens & Sellwood 2007) and our own Fourier component decomposition software (Finlez et al. 2018) based on the the linear perturbation analysis described in Wong et al. (2004), Fathi (2004).

The Diskfit⁸ package can be used to fit both the image and the velocity field of a galaxy. In imaging mode, an input image is fit with one or more of a bulge, disk, and bar, resulting in estimates of the relative flux and morphology (ellipticity, brightness profile, and PA) of each component. In velocity mode, Diskfit models asymmetric rotation-dominated velocity fields using a combination of tangential and radial perturbations to a fitted circular velocity model. We fit our CO J:2-1 velocity field using Diskfit considering only m = 2 potential perturbation (i.e., bars) modes and using the galaxy nuclear position, galaxy PA, galaxy inclination and bar PA as fixed values (mm continuum peak position, 45°, 33°, and 0°, respectively). The best fit model obtained by Diskfit, and the velocity residuals (observed − Diskfit model) are shown in the top panels of Fig. 10. The best fit model from Diskfit differs from our toy rotation model (Fig. 7, top left and Sect. 3.2) in that the apparent rotation axis moves to a slightly smaller PA in the inner 4″, the inner 1″ shows twisted isophotes, and there is a resonance at ∼4″, which mainly falls in a region where we do not have observed velocities due to low signal-to-noise.

Fig. 10. Top panel: Diskfit model, with perturbations from a bar at PA 0°, fit to the CO velocity field (left panel) and the resulting residual CO velocity field (observed − model) (right panel). Bottom panel: illustrative bar-perturbed velocity fields resulting from our linearized epicyclic perturbation models (see text) and using λ = 0.2 and Ω = 120km s−1 kpc. The left (right) panel shows the results when setting the intrinsic rotation curve as that from C14 (our Bertola model fit to the stellar velocity field).

The residual velocity map obtained after subtracting the Diskfit (Fig. 10 top right panel) shows smaller deviations from systemic as compared to the velocity residual made from our rotation-only model (top right panel of Fig. 7), especially to the SE of the nucleus, and in general in the inner arcsecs. However, the Diskfit model still does not attain the highest velocities seen in the inner arcsecs. We note that Diskfit only allows us to change basic photometric parameters of the galaxy, e.g., disk PA and inclination and bar PA, and the input observed velocity field. Since all these are relatively well defined for NGC 1566, we are unable to further fine-tune the results of Diskfit.

To better illustrate the differences between the Diskfit model and the observed velocity field, we plot the Diskfit model (orange lines) on the observed pv diagrams at several relevant PAs (Fig. 12). We immediately note that the best-fit Diskfit model was derived from the velocity field (intensity weighted average velocity at each spatial pixel) rather than the full datacube, so that comparing the model directly to the pv diagram is not really fair. Instead it is more correct to compare the model (orange lines) with the velocities from the moment 1 map (intensity weighted velocity; black dashed lines in the figure). While the Diskfit model slightly overpredicts the velocities seen along the major axis (top right panel), and the pattern of the velocities seen along the minor axis (bottom middle panel), it fails to predict (by a factor ∼2) the high peak velocities seen along the minor axis (PA = 135°), or along PAs 0° and −15°. That is, the bar perturbations are unable to explain the ∼90 km s⁻¹ radial velocities seen in the inner 1″ to the NW and SE of the nucleus.

While Diskfit models the observed velocity field with a base rotation model plus perturbations in radial and tangential velocities (one component each in the case of m = 2 modes), it does not use (or at least does not provide details to the user) a physically-based model with, e.g., a mass-based rotation curve or a fixed bar pattern speed. We thus additionally model the observed CO velocity field with linearized epicyclic perturbations produced by a bar (for details see, e.g., Wong et al. 2004; Fathi 2004) applied to a physically derived rotation curve (an exponential disk whose mass is constrained by near-IR photometry) for a given bar pattern speed (Ω), damping factor (associated to a frictional force; λ), and bar PA and ellipticity (Finlez et al. 2018). Our code is based on the algorithms proposed in Franx et al. (1994), Wong et al. (2004) and Fathi (2004). It is important to note that we specifically use the “m = 2” potential (relevant for bars) which introduces changes in the 1st and 3rd harmonic coefficients (Schoenmakers et al. 1997), and that this perturbation analysis is valid only for “weak” bars, i.e., when the bar potential does not dominate the disk potential).

We used Diskfit to decompose an IRAC 3.6 μ image of NGC 1566, obtained from NED, into bulge, disk, and bar components. The galaxy and bar PAs were fixed and other parameters allowed to vary. We further assume a constant mass to light (M/L) ratio for all three components. We find that the disk contains 58% of the total mass (light) of NGC 1566 with a bar to disk mass ratio of 0.7. Given the total mass of NGC 1566 derived by Sheth et al. (2010), the disk mass is 2.2 × 10¹⁰ M_⊙. Alternatively, the 3.6 μ disk luminosity with a M/L_3.6 μ ratio of 0.47 (McGaugh & Schombert 2014) implies a disk mass of 4.4 × 10⁹ M_⊙. An exponential disk with these total masses was then used to derive a first-guess intrinsic (i.e., before bar perturbations) axisymmetric rotation curve (details in Finlez et al. 2018). The disk mass was then slightly adjusted (to 6 × 10⁹ M_⊙) in order to better fit (by eye) the model rotation curve of the CO in the inner disk (i.e., ModC2014) or ∼2 × 10⁹M_⊙ to agree with our best-fit Bertola model (Bertola et al. 1991, Eq. (2)) to the stellar velocity field. It should be noted that these masses are an order of magnitude lower than that predicted by Korchagin et al. (2000) (1.7 × 10¹⁰ M_⊙).

We then ran our linear perturbation code, in m = 2 mode, using as inputs the intrinsic axisymmetric rotation curve(s) derived above, the PAs of the galaxy disk and bar, and the bar ellipticity (ϵ = 0.42), the latter derived from our Diskfit decomposition. The bar pattern speed (Ω) and the damping factor (λ) were allowed to vary. The resultant model velocity fields for a range of values of Ω and λ, when using the best-fit Bertola model of the stellar velocity field, are shown in Fig. 11. The most notable effect of varying Ω is the change in the radii of the resonances. Gas orbits change abruptly when crossing these resonances; the effect of increasing damping (increasing λ) is to smooth out these large swings in the orbits. Most of the panels in Fig. 11 show the characteristic “butterfly” pattern expected from bar perturbations. However, for this butterfly pattern to fall within the central ∼2″ as observed, i.e., to explain the innermost high-velocity features, one requires extremely high (∼300 km s⁻¹ kpc⁻¹) bar pattern speeds. Alternatively, the intrinsic rotation curve requires to rise slower or flatten at lower velocities. We must note that the uncertainty in the distance to NGC 1566 (see Sect. 1) plays a significant role in the bar pattern speeds used here. If a distance of 20 Mpc is used for NGC 1566 then the bar pattern speeds we list here would halve, so that less extreme bar pattern speeds could replicate the observed resonances. In any case, even if the resonance radii are matched, the pattern of the model velocities are significantly different from the observed CO velocity field (and the larger scale Hα velocity field from Pence et al. (1990): specifically at higher pattern speeds the strongest perturbations inside the inner resonance are in PA ∼ 100°, offset from the PA of our posited outflow, and beyond the inner resonance the kinematic axis of is highly curved, starting at PA ∼ 0° and then curving to the observed PA of the galaxy.

Fig. 11. Bar-perturbed velocity fields obtained with our epicyclic perturbation models applied to our best fit Bertola stellar rotation curve (see text), on varying the bar pattern speed (left to right; Ω in units of km s−1 kpc−1) and the dimensionless damping parameter (top to bottom; λ). The FOV of each image is 12″ × 12″ and major tick marks are shown every 1″. All velocities follow the same color bar shown on the top.

For illustration, we compare the predictions of the perturbation model which uses the Bertola best fit model to the stellar velocity field as the intrinsic rotation curve, and parameters Ω = 120 km s⁻¹ kpc⁻¹ and λ = 0.2 (the model shown in the third row, third column of Fig. 11) with our observed pv diagrams in Fig. 12. While this model does not well fit the observed velocity field, it uses a pattern speed argued for in C14 (based on corotation placed at the bar end) and a damping parameter within the range of values typically invoked for other well studied galaxies (between 0 and 0.5, Wada 1994; Fathi et al. 2005), and is thus a good reference point.

Fig. 12. Position-velocity diagrams, as in Fig. 9, but here the overplotted lines show the velocity of the corresponding pixel in the CO moment 1 map (dashed black lines; i.e., the flux-weighted average velocity at a given pixel), the best-fit velocity field from Diskfit (solid orange line), and for comparison the λ = 0.2 and Ω = 120 km s−1 kpc model from our linear perturbation analysis (bottom right panel of Fig. 10; solid green line). To better trace the highest velocity components in the CO moment 1 (velocity) map, we used a cutoff of 8 mJy beam−1 per channel (∼4σ) to create the moment 1 map which is overplotted here and used as the input map to Diskfit. The PA of the “slit” along which the pv diagram was extracted is marked above each panel. Bottom right panel: the moment 0 map of CO J:2-1 together with the positions of the slits used to create the pv-diagrams in the upper row panels.

As a further illustration, the bar-perturbed velocity fields for Ω = 120 km s⁻¹ kpc, λ = 0.2, and for both options of the intrinsic rotation curve (gas- and stellar-rotation curve models) are shown in bottom panels of Fig. 10. The model which uses an intrinsic rotation curve similar to that of the gas (bottom left) exhibits a resonance at ∼4″, similar to that obtained by Diskfit (left top), but presents less pronounced nuclear distortions as compared to the Diskfit model. Using the slower rising Bertola (stellar velocity) model as the intrinsic rotation model (bottom right) changes the position of the resonance to (∼$2\stackrel{″}{.}6$), but also gives lower velocity distortions along the minor axis, or rather the higher velocities seen in the observed velocity field (∼1″) are further out in the model (∼3″): therefore, to spatially matching these velocity distortions requires higher bar pattern speeds or a slower rise in the intrinsic axisymmetric rotation curve. Apart from the mismatch in resonance radii, these two panels also clearly illustrate the mismatch between the observed and modeled velocity fields noted above; specifically, the misalignment of kinematic axes inside the resonance (related to the posited outflow), and the large curvature in the kinematic axis beyond the inner resonance.

Overall, we are unable to convincingly fit the observed CO kinematics with perturbations produced by the large scale bar. Our linearized epicyclical bar perturbation models, which use realistic values for the intrinsic rotation curves and the bar pattern speed (with the caveat of the uncertainty in the distance to NGC 1566), are able to reproduce the amplitudes of the inner perturbations. However, the resonances are produced further out than observed, and the velocity changes are not as sharp as observed. Higher bar pattern speeds, perturbations by an inner bar in a different PA, or different intrinsic rotation curves, would be required. Diskfit reproduces reasonably many of the observed features in the velocity map, and at first glance provides a reasonable explanation for the perturbations observed, even if the amplitude of these perturbations is not as high as observed. However, we are wary of the results of Diskfit for two main reasons. First, Diskfit does not provide feedback on the underlying physical parameters of the resultant model, and thus, e.g., we are unable to evaluate whether the bar pattern speed used is physical and second, we have applied Diskfit to about a dozen galaxies for which we have disturbed optical emission line kinematics over the inner 5″ of the galaxy and almost always found relatively good fits (Schnorr-Muller, priv. commun.), even though our detailed multi-component analysis either found the perturbations to be due to bars (Schnorr-Müller et al. 2017a) or outflows and/or streaming inflows (most other cases, eg., Schnorr-Müller et al. 2017b). In fact, Spekkens & Sellwood (2007) obtained a good fit to the velocity field in NGC 2976, but to conclude that the perturbations were due to the bar, they confirm their existence at the PA predicted by the model, according previous photometry. This consistently good performance of Diskfit makes it more difficult to believe that the fits are truly consistent and physically motivated rather than empirical best fits to distorted velocity fields. We emphasize that we are not stating that bar-related perturbations do not exist in the velocity field, rather we argue that bar-related perturbations are not the unique and dominant driver of the observed nuclear perturbations in the CO velocity maps, and it is most likely that the nuclear perturbations are produced by an AGN-driven outflow.

3.7. Modeling observed velocities: CO J:2-1 streaming?

The presence of putative outflows and/or bar related perturbations (previous sections) makes it difficult to search for signatures of streaming motions in the inner few arcsec in velocity maps or even residual velocity maps (e.g., Fig. 7). That is, since the velocity maps show the intensity weighted average velocity of the spectrum corresponding to each spatial pixel, they are likely dominated by the “outflow” signature in the inner ∼3″. Further, when the intensity of the gas in rotation dominates that of the gas in inflow, velocity residual maps will not show strong indications of the inflow. It is thus important to examine the velocity profile of each pixel or aperture in order to separate outflows, rotation, and streaming inflows. Ideally, one requires a complete velocity field to analyze the azimuthal average of radial gas velocities at each radius. While this is often possible in the case of ionized gas, molecular gas is often detected only over a limited range of azimuths at each radius. In the case of NGC 1566, our CO velocity maps are “complete” out to a radius of ∼3″, beyond which the velocity filling factor is ∼5–40%. Under the assumption that these detected CO regions dominate the CO flux at their respective radii, the detected clumps or arms can still be used to constrain the presence of streaming flows. Given the above, we model streaming inflows with a very simple toy model in which the inflow is assumed to have a constant radial inflow velocity (which we fix to 50 km s⁻¹ after initial inspection of the results). This velocity is then projected and added to the projected radial velocity expected from our rotational model. We first examine the spectral profiles in apertures along the inner spiral arms (left panel of Fig. 13). The CO J:2-1 spectra extracted from these apertures are shown in the right panel of Fig. 13, together with the average radial velocities expected from our models of rotation, outflows, and streaming (and combinations thereof).

Fig. 13. CO spectral profiles in apertures along the inner spiral arms. Left panel: location of each circular aperture ($0\stackrel{″}{.}2$ in radius) overlaid on the CO flux map. For reference, the center of the galaxy, major axis, and the near and far sides of the disk, are indicated. Apertures are numbered 0–8 with aperture 0 being the farthest aperture on the E arm, aperture 4 the nuclear one, and aperture 8 the farthest aperture on the W arm. Right panel: extracted CO J:2-1 spectra; thick lines are used for the spectra corresponding to apertures from the W arm (Apertures 5–8), and different colors are used for each spectrum, following the color bar on top of the panel. Symbols with the corresponding color (plotted at the y value of the peak flux density of the spectrum) denote the radial velocities expected in that aperture for our rotation model (hourglasses), our outflows + rotation model (asterisks), and our rotation plus radial streaming inflow model (open circles). In some apertures, adding outflows and/or streaming inflows does not change the predicted radial velocity; this is a result of projection effects and/or the fact that our outflow model has zero velocity beyond ∼2″ from the nucleus.

Analyzing the CO J:2-1 spectra in the right panel of Fig. 13 we see the following: (a) The nuclear spectrum (light green) is clearly double peaked: the peaks at V_rad ≈ 65 km s⁻¹ are attributed to putative nuclear outflows (previous sub-sections), and the highest velocities seen are V_rad ≈ 100 km s⁻¹. It should be noted that there is a plateau of CO J:2-1 emission at lower velocities, potentially from gas rotating in the disk (recall that our rotation models predict velocities of V_rad ≈ 0 − 50 km s⁻¹within this aperture). (b) The off-nuclear apertures show profiles with varied shapes and widths, and the off-nuclear apertures which intersect the galaxy major axis clearly show multiple velocity components. (c) For the W spiral arm (spectra plotted in thick lines) there is a large mismatch between the spectral profiles and the expectations from our rotation model only (hourglass symbols in Fig. 13). In order of increasing distance from the nucleus along the W spiral arm, gas in the first aperture rotates faster than predicted; the gas profile in the second aperture has a strong red shoulder at a velocity consistent with rotation, while the profile peak is offset ∼30 km s⁻¹ to the blue; gas in the third aperture rotates slower than predicted; and gas in the farthest aperture is centered at zero velocity since the aperture lies on the minor axis of the galaxy. (d) For the E spiral arm (spectra plotted in thin lines) the profiles are more centered on the predictions of our rotation only model. In order of increasing distance from the nucleus along the E spiral arm, gas in the first aperture (which includes the edge of the strong CO knot ∼1″ from the nucleus to NE) lies close to the prediction of the rotation model but a clear blue shoulder is seen; gas in the second aperture shows a clear double-peaked profile with the expected rotation velocity lying in the middle of the two peaks; gas in the third aperture also shows a broad blue shoulder. (e) Including our decelerating outflow model (i.e., using the predictions of rotation plus outflows; asterisks in the figure) the model predictions change significantly only for the two off-nuclear apertures closest to the nucleus. Here the aperture to the S (solid yellow spectrum in the figure) fits the prediction satisfactorily, i.e., as if almost all gas is in outflows, but to the N (cyan spectrum), while the aperture profile shows an extra blue wing in the correct velocity direction for outflows the magnitude of the offset does not fit well with our model, and the gas seems to be dominantly in rotation rather than outflows. (f) Using a model which sums our streaming toy model to pure rotation (open circles in the figure), we see that the apertures in the W arm are inconsistent with the predictions of streaming inflows: the mismatch between the spectrum peak and the prediction increases when changing from rotation only to rotation plus radial inflow. In the E arm, however, the profiles are in general as consistent with the streaming inflow+rotation model as with the rotation only model, i.e., while the profile peaks are consistent with rotation, the prominent shoulders on these profiles are roughly consistent with streaming inflows. This is also clearly seen in the skewness map (Fig. 6 where almost the full inner spiral arm structure shows a blue “skewness”, independent of being on the near- or far-side of the galaxy disk.

We also show the spectra in apertures along the PA of the large-scale bar in Fig. 14; here we use larger apertures ($0\stackrel{″}{.}5$ radius) to obtain a higher signal-to-noise. Here the peak of the spectral profiles are consistent with the predictions of rotation+outflow and a shoulder is seen roughly at the predicted velocity of rotation plus radial streaming inflows.

Fig. 14. Same as Fig. 13, but for apertures oriented along the bar. Apertures are numbered 0 (southernmost) to 4 (northernmost), and apertures here are each $0\stackrel{″}{.}5$ in radius.

To test for streaming inflows along spiral patterns further from the nucleus, we also examined the spectra from apertures along the spiral patterns SE to S of the nucleus which connect to the inner spiral arm discussed above, and spectra in apertures tracing spiral structure to the N and NW of the nucleus, which connect to the W inner spiral. We do not show these spectra as the essential results can be seen in the skewness map of the CO J:2-1 line (Fig. 6, left panel): here blue (red) colors represent pixels where the CO spectral profile is skewed in the sense of having excess emission towards the blue (red) of the weighted mean velocity at that pixel. To avoid contamination by noise the skewness was calculated using the spectral profile down to 10% of the peak flux. If most of the gas follows regular rotation and a smaller fraction of gas participates in a radial streaming inflow then we expect blue (red) skewness on the far (near) side of the galaxy. The inner 2″ shows profiles skewed consistently to the blue. show that the gas is mainly consistent with rotation. Two spiral arm sections at ∼10″ from the nucleus (roughly at the putative location of the Inner Lindblad Resonance (ILR) of the large-scale bar (Comerón et al. 2010, C14) to the SSE and to the N show a skewness consistent with inflows. All spectra in these arms show non-symmetric profiles which are likely from multiple velocity components. In fact velocity differences between the peak of the profile and the shoulders to the red match well the predictions of the offset (from rotation) velocity expected from streaming inflows.

3.8. Molecular mass in the inner kiloparsec

Molecular gas mass is typically estimated from the CO luminosity (Solomon & Vanden Bout 2005) using $M_{\mathrm{mol}}{M}_{\odot }={\alpha }_{\mathrm{CO}}\times L_{\mathrm{CO}}\prime$ where:

$$\begin{array}{c}\hfill L_{\mathrm{CO}}^{\prime }=3.25\times {10}^7\times S_{\mathrm{line}}\mathrm{\Delta }\nu \frac{D_L^2}{{(1+z)}^3{\nu }_{\mathrm{obs}}^2}·\end{array}$$(2)

Here, $L{\prime }_{\mathrm{CO}}$ has units of K km s⁻¹ pc², S_lineΔν is the integrated flux density of the CO J:1-0 line in Jy km s⁻¹, D_L is the luminosity distance in Mpc, z is the redshift, and ν_obs is the observed frequency in GHz. There remains significant debate on the value of α_CO, and we use the value α_CO= 4.3 M_⊙ (K km s⁻¹ pc²)⁻¹ as suggested by Bolatto et al. (2013) for the Galaxy and other nearby spiral galaxies which are not extreme starbursts. It should be noted that $L_{\mathrm{CO}}^{\prime }$ is directly proportional to the surface brightness in K units, and therefore the $L_{\mathrm{CO}}^{\prime }$ ratio of two CO J transitions gives the ratio of their surface brightness temperatures. Furthermore, $L_{\mathrm{CO}}^{\prime }$ is constant for all J levels if the molecular gas emission comes from thermalized optically-thick regions, i.e., the brightness temperature and line luminosity are independent of J and rest frequency for a given molecule (Solomon & Vanden Bout 2005).

Since we observed the CO J:2-1 line, we require to convert L′_{COJ : 2 − 1} to L′_{COJ : 1 − 0}, a conversion which depends on the physical conditions of the gas. Bajaja et al. (1995) have observed the CO J:1-0 and CO J:2-1 lines in NGC 1566 at low resolution using the SEST telescope, and C14 have presented ALMA CO J:3-2 observations of NGC 1566. The nuclear CO J:2-1 / CO J:1-0 intensity ratio found by Bajaja et al. (1995) is 1 in temperature units: as expected from thermalized optically-thick gas. C14 compared the Bajaja et al. (1995) CO J:2-1 integrated flux densities with their ALMA-derived CO J:3-2 integrated flux densities and inferred that the latter was missing some flux; they thus also assumed that the gas is thermalized and optically-thick in their calculation of molecular gas masses.

To further constrain the observed flux ratios, we downloaded the CO J:3-2 data of C14 from the ALMA Science Archive⁹ and created a Moment 0 map. For thermalized optically-thick gas one expects that the integrated flux density (in units of Jy km s⁻¹) of each CO J transition varies as ν². We find that in the inner ∼3″ radius (the inner disk), the integrated CO J:2-1 flux density (251 ± 25 Jy km s⁻¹) is 1/2.3 times that of CO J:3-2. Over a larger field, 12″ × 12″ in size, the total flux density of the CO(J:2-1) line (406 ± 41 Jy km s⁻¹) is half that of CO J:3-2. We thus, as in C14, assume thermalized optically-thick gas, i.e., L′_{COJ : 2 − 1}/L′_{COJ : 1 − 0} = 1. Under this assumption, using a value of α_CO listed above, and D_L = 10 Mpc, the molecular gas mass in the inner spiral arms is (6.6 ± 0.7)×10⁷ M_⊙, and the total molecular gas mass in the central 12″ × 12″ region is (1.1 ± 0.1)×10⁸ M_⊙.

Estimating the molecular gas mass in the nuclear outflows of NGC 1566 is more difficult. Recall that outflowing molecular gas is clearly detected out to ∼2″ (Sect. 3.3), and perhaps out to 5″ (Sect. 3.7). To estimate the mass and momentum of the outflowing gas in the unresolved nucleus we extract the nuclear CO J:2-1 emission spectrum in a circular aperture of $0\stackrel{″}{.}2$ in radius (similar to our synthesized beam area; Fig. 15). This spectrum appears to be made of three distinct components: a strong blue Gaussian representing emission from the blue outflow ((2.0 ± 0.2)×10⁵ M_⊙), a weaker red Gaussian representing emission from the red outflow ((1.0 ± 0.1)×10⁵ M_⊙) (Sect. 3.3) and an intermediate velocity region (green line in the figure; radial velocities between ±40 km s⁻¹; (0.8 ± 0.1) × 10⁵ M_⊙). This intermediate velocity region could originate in (a) gas in a spherical outflow in the nucleus since variation in the projection angles to the line of sight will give emission at all velocities; (b) unresolved emission from gas in solid body rotation-only; (c) emission from dispersion-dominated gas in the nucleus. With the masses calculated above, and using the (deprojected) mean flux-weighted velocity of each of the above components under the assumption that the outflows are in the plane of the disk we can estimate the outflow momentum in the unresolved nucleus, obtaining (−19.3 ± 2.8) × 10⁶ M_⊙ km s⁻¹and (12.1 ± 1.7) × 10⁶ M_⊙ km s⁻¹for the blue and red outflows in the inner (unresolved; ≲0. pc) nucleus, respectively.

Fig. 15. CO J:2-1 spectrum (black) in a circular aperture centered on the nucleus with radius $0\stackrel{″}{.}2$. The spectrum was divided by eye into three parts: the pure outflow components (blue and red) and an intermediate velocity region (green; see text) between the blue and red outflow peaks.

Since the outflows extend beyond the inner (unresolved) nucleus we can extend this outflow mass analysis to larger scales. From the bottom right panel of Fig. 9, we can distinguish that the putative outflows along the minor axis show two stages: (1) an initial stage between radii of 0 and ∼1″ which show velocities close to the peak velocity, and (2) a final stage with monotonically decreasing velocities (down to ∼0) between 1″ and $1\stackrel{″}{.}5$ from the nucleus. For these two stages (but excluding the innermost $0\stackrel{″}{.}2$ considerated as the nuclear outflows above) we obtain total outflow mass values of (9.3 ± 0.9) × 10⁶ M_⊙ and (11.6 ± 1.2) × 10⁶ M_⊙, respectively.

In the absence of clear signatures of inflowing gas we do not attempt to estimate an exhaustive streaming inflow rate. Instead we show two examples of spectra in apertures for which we found the strongest signatures of non-circular motions which could be explained by our toy streaming model (see Sect. 3.7) and use these to roughly estimate the mass involved in the inflow. In the E inner spiral we use the spectra of apertures 1 and 2 of Fig. 13. These spectra are shown in Fig. 16, with the red overlay denoting the velocities over which the emission is from gas potentially participating in the streaming inflows. For each of these apertures we find masses of ∼(1.1 ± 0.4)×10⁵ M_⊙ potentially participating in a streaming inflow. We performed the same exercise with four apertures in the outer extension of the E spiral arm (not shown) and find a median mass of ∼(4 ± 0.4)×10⁴ M_⊙ potentially participating in a streaming inflow.

Fig. 16. CO J:2-1 spectra extracted from apertures 1 and 2 of Fig. 13 (the E inner spiral arm) are shown in black. The red overlay marks the velocities over which emission comes from gas potentially participating in streaming inflows.

3.9. Mass of the SMBH in NGC 1566

To estimate the black hole mass in NGC 1566, we use the empirical correlation between SMBH mass and stellar velocity dispersion σ_⋆ (Ferrarese & Merritt 2000; Gebhardt et al. 2000; Tremaine et al. 2002; Gültekin et al. 2009; Kormendy & Ho 2013). It should be noted that this M − σ_⋆ relationship can also be reliably applied to nearby AGNs (Nelson et al. 2004; Woo et al. 2010; Graham et al. 2011). While there are several versions of the M − σ_⋆ relation we use that of Gültekin et al. (2009) for their entire sample for several reasons; they demonstrate that previous studies are biased by considering culled samples according whether the SMBH sphere of influence is resolved (see Sect. 4 of their paper). Therefore we get:

$$\begin{array}{c}\hfill \left(\frac{M_{\mathrm{BH}}}{M_{\odot }}\right)={10}^{(8.12\pm 0.08)}·{\left(\frac{\sigma }{200}\right)}^{(4.24\pm 0.41)}\end{array}$$(3)

where σ, is the central velocity dispersion in km s⁻¹. To obtain the latter, we fit the integrated spectrum from our GMOS/IFU datacube with the program pPXF described in Cappellari & Emsellem (2004). The spectrum from 5700 Å to 6300 Å was fit in pPXF using SSPs templates from Bruzual & Charlot (2003), obtaining a central stellar velocity dispersion of σ_⋆ = 116 ± 9 km s⁻¹. This value is consistent with that obtained using σ_⋆ from Bottema (1992) for the innermost region in NGC 1566, and higher than the bulge stellar velocity dispersion used in other previous studies (van der Kruit & Freeman 1984; Nelson & Whittle 1995; Woo & Urry 2002).

Our measured value of σ implies an estimated super massive black hole mass of M_BH = 1.3 ± 0.6 × 10⁷ M_⊙.

3.10. Bolometric luminosity and accretion, inflow, and outflow rates

Given the SMBH mass estimated above, for which the Eddington Luminosity is L_Edd = (4.0 ± 1.8) × 10¹¹ L_⊙, we would like to ascertain the bolometric luminosity of the AGN L_Bol as well the Eddington ratio (${\text{l}}_{\text{Edd}}=\frac{L_{\text{Bol}}}{L_{\text{Edd}}}$).

Several differing values for L_Bol, for NGC 1566 have been obtained in previous studies. Here we use two different approaches to estimate L_Bol: from the nuclear [O III] luminosity, and from the nuclear hard X-ray luminosity. The nuclear 2–10 keV luminosity as measured by XMM and scaled to our adopted distance (10 Mpc) is L_X = (7 ± 3)×10³³ W (Levenson et al. 2009). Using the hard X-ray to bolometric luminosity conversion of Ulvestad & Ho (2001) (L_Bol = 6.7 × L_{X(2 − 10keV)}) we obtain L_Bol = (4.69 ± 2)×10³⁴ W. The nuclear [O III] flux (Moustakas et al. 2010) at our adopted distance gives L_{[O III]} = (8.2 ± 1.8) × 10³¹ W. Using the scaling of Heckman et al. (2004) modified as recommended in Dumas et al. (2007) (L_Bol = 90 × L_{[O III]}) we obtain L_Bol = (7.42 ± 1.62)×10³³ W. These two methods thus give relatively consistent values for L_Bol and we thus adopt the mean value of L_Bol = (2.7 ± 1.3)×10³⁴ W which implies l_Edd ≈2.0 × 10⁻⁴, i.e., a relatively low-efficiency regime for the SMBH, considering the fact that in the most active galaxies, gas is accreted onto the SMBH in a efficient regime with ratios between 0.01-1 (Khorunzhev et al. 2012). We can now estimate the mass accretion rates as follows:

$$\begin{array}{c}\hfill \dot{m}=\frac{L_{\mathrm{Bol}}}{c^2\eta }\end{array}$$(4)

where η is the accretion efficiency which in nearby galaxies with geometrically thin, optically thick accretion disks, is typically taken to be 0.1 (Soltan 1982; Fabian & Iwasawa 1999; Yu & Tremaine 2002; Davis & Laor 2011). We thus get an accretion rate of ṁ = (4.8 ± 2.3) × 10⁻⁵ M_⊙ yr⁻¹.

We note that very discrepant values of L_Bol were obtained by Woo & Urry (2002) who integrated the flux in the spectral energy distribution (SED) of NGC 1566 using data from NED, after averaging multiple datapoints in the same band and correcting for dust. The value of L_Bol they obtained, scaled to our adopted distance, is L_Bol = 7.2 × 10³⁶ W. This implies l_Edd = 0.05 and an accretion rate of ṁ = 1.3 × 10⁻² M_⊙ yr⁻¹, 3 orders of magnitude higher than the values obtained above. Two potential reasons for this large discrepancy are (a) flux variability in the hard X-ray, see, e.g., Landi et al. (2005); (b) the use of large (galaxy-wide), rather than nuclear, apertures for the data in NED. In Sect. 3.8 we roughly estimated the mass of the molecular gas potentially participating in streaming inflows using six apertures along the SSE to E spiral arm. This can be used to estimate a rough mass inflow rate. The radius of each aperture ($0\stackrel{″}{.}2$) corresponds to a linear diameter of 0. pc. A streaming inflow velocity of ∼50 km s⁻¹ would imply an aperture crossing time of ∼3.8 × 10⁵ yr, and thus a mass inflow rate of 0.1 M_⊙ yr⁻¹ along this spiral arm.

A similar procedure can be used to estimate the outflow rates. We use the molecular masses deduced for the blue and red components of the outflows in the nuclear ($0\stackrel{″}{.}2$ radius) aperture (see Sect. 3.8 and Fig. 15). In Sect. 3.3 we have argued that the velocity of the outflows in the nucleus is ∼180 km s⁻¹ in the plane of the disk, which gives a crossing time of around 5.42 × 10⁴ yr for the nuclear $0\stackrel{″}{.}7$ aperture used. The nuclear molecular mass outflow rates are thus 3.7 M_⊙ yr⁻¹ for the blueshifted component and 1.9 M_⊙ yr⁻¹ for the redshifted component. We note that these outflows do not appear to escape from the nucleus, but instead “pile up” in the inner gas ring (see the CO residual map and pv diagrams). We can also estimate the kinetic power associated with cold molecular outflow following, e.g., Harrison et al. (2014) and Lena (2015):

$$\begin{array}{c}\hfill {\dot{E}}_{\mathrm{out}}=\frac{{\dot{M}}_{\mathrm{out}}}{2}(v_{\mathrm{outflow}}^2+3{\sigma }^2).\end{array}$$(5)

Assuming a nuclear velocity dispersion of 60 km s⁻¹(see Fig. 6) and, as mentioned above, v_outflow = 180 km s⁻¹, we obtain a total outflow kinetic power (over both blue and red outflow components) of ${\dot{E}}_{\mathrm{out}}=7.62\times {10}^{33}\mathrm{W}$. The ratio between the outflow kinetic power and the AGN bolometric luminosity is thus $\frac{{\dot{E}}_{\mathrm{out}}}{L_{\mathrm{Bol}}}\approx 0.28$. This is significantly higher than the values obtained in previous studies of nearby active galaxies which found $\frac{{\dot{E}}_{\mathrm{out}}}{L_{\mathrm{Bol}}}$ ranging between 0.1 and 10% (e.g., Storchi-Bergmann et al. 2010; Müller-Sánchez et al. 2011, 2016; Harrison et al. 2014; Lena 2015) or even lower (e.g., Barbosa et al. 2009, less than 0.01%). However all of these studies have measured outflow masses using ionized gas, which is expected to be a minor fraction of the total gas content and also, using measurements on larger spatial scales than those considered here ($0\stackrel{″}{.}4$ or ∼24 pc).

4. Summary and conclusions

We have analyzed the kinematics in the inner kiloparsec of the nearby active galaxy NGC 1566, using ALMA observations of CO J:2-1 along with GMOS/IFU data of ionized gas emission lines and stellar absorption lines. Our results allow us to conclude the following:

NGC 1566 presents a cold molecular dense disk in the inner 144 pc, along with a clear two-arm spiral structure in the inner 96 pc. These structures are also seen in our ionized gas ([N II]) images, and in previous studies using CO J:3-2 and optical/IR emission lines. The inner spiral arms and dense disk appear to have some continuity with more extended CO J:2-1 spiral arms which extend out of the inner disk to larger scales, and which coincide with dust lanes seen in HST images.

Ionized gas and stars are detected over the full IFU FOV at high signal-to-noise ratios. The [N II] emission is peaked at the nucleus and is also strong at the known star-forming region to the SW, about $1\stackrel{″}{.}5$ from the nucleus.

The superior spectral resolution (∼2.6 km s⁻¹) and image fidelity in this new CO J:2-1 datacube provides improved constraints on the nuclear kinematics. Furthermore, our use of pv diagrams (rather than only intensity weighted velocity maps) allows the full exploitation of our velocity resolution. The molecular gas kinematics of the inner disk is dominated by rotation with peak velocities of ∼±140 km s⁻¹. The residual velocity field shows clear signs of non-rotational motion especially in the innermost ∼2″ region. We argue that the strongest deviations are the result of nuclear outflows, though the strong two-arm inner spiral structure and large-scale bar also play a role. The CO rotational curve over the inner ∼3″ is asymmetric, which could be the consequence of a warped disk. PV diagrams at PAs between −15° and 15° show some discrepancies from our rotation + outflow model at radii ∼$1\stackrel{″}{.}5$ from the nucleus. These discrepancies could be explained by perturbations due to the barred potential and/or streaming velocities.

We argue for the presence of nuclear (inner ∼2″) molecular gas outflows in the plane of the galaxy disk, without ruling out the presence of bar- or spiral-related perturbations (see Sect. 3.6 for a detailed analysis of the potential bar influence). The arguments for a nuclear outflow include the following: (a) several previous authors have claimed kinematic and morphological evidence of outflows in the NLR, which likely intersect the galaxy disk given the observed geometries; (b) the nucleus shows a double-peaked profile with FWZI 200 km s⁻¹, higher than that seen in the lower fidelity maps of C14, and which imply inclination corrected outflow velocities of up to 180 km s⁻¹ in the plane of the disk. A large angle between the outflow axis and the line of sight is unlikely since this would imply extremely high true outflow velocities. Attributing these velocity features to other perturbations in the plane of the disk and along the minor axis requires true radial velocities around ∼80–130 km s⁻¹ in the nuclear region where rotation velocities are expected to be ≤40 km s⁻¹. In Sect. 3.6 we analyzed the large-scale bar perturbation and its implications, showing that it does not reliably produce the morphology and the large perturbations seen in the observed velocity field; (c) the pv diagram along the minor axis (bottom right panel of Fig. 9) connects high-velocity components to the zero-velocity components seen at r≈$1\stackrel{″}{.}8$ on both sides of the nucleus. The regions with the greatest velocity deceleration are correlated with the regions brightest in CO and thus richest in molecular gas, indicating that the outflows decelerate due to mass loading, which is clearly seen in the inner arcsec to the N and in the r ∼ 1 − 2″ range to the S; (d) the pv diagrams (Fig. 9) show velocity deviations which are consistent with radial outflows not just in the minor axis, but also in the plane of the disk over all PAs and over apertures at distances of several synthesized beams from the nucleus (see Sect. 3.5); and (e) the consistence with the evidence of a nuclear spherical (or bipolar) outflow in our ionized gas kinematics. We also argue that the molecular outflow is primarily detected within the galaxy disk. The supporting arguments for this include the low velocity dispersion of the molecular gas, the absence of evidence for radio jet-related outflows, the posited deceleration of the outflow, which would be consistent with the high density of gas in the galaxy disk. Other potential scenarios as warped disk or non-coplanar disk cannot be constrained by us due to the limited resolution, the sparse velocity field and the lack of a reliable circular rotation model.

While the stellar and ionized gas kinematics predominantly traces rotation in the galaxy disk, we find evidence of two components in the ionized gas: a narrow (σ ∼ 60 km s⁻¹) component detected over almost the full field of view and which traces rotation in the disk, and a broad component detected in the inner ∼3″ radius. We postulate that the broad component of ionized gas is part of a nuclear spherical outflow. The broad component of the [N II] emission line does not participate in the rotation of the galaxy disk. Instead it is preferentially blueshifted, as expected since nuclear dust obscures the receding side of the spherical outflow. Its velocity dispersion, velocity, and the correlation with strong and weaker dust features in the nucleus all strengthen the arguments for the spherical ionized gas outflow. Even so, we do not rule out the possibility of a bipolar outflow which currently is impeded by a poor resolution in the GMOS data. We are unable to test for deceleration in this outflow.

We have constrained and analyzed the circular vs. perturbed kinematics of the CO gas in the disk using Kinemetry, Diskfit, and linear perturbation theory. Large radial velocity perturbations are clearly required. The Diskfit model better fits the observed data as compared to the pure rotation model, but still does not attain the highest velocities seen in the inner arcsecs and fails to predict the high peak velocities seen along the minor axis in pv diagram (bottom middle panel in Fig. 12). The relatively good fits obtained by Diskfit do not necessarily imply that the perturbations are bar related. Our concerns here include the lack of definition of the physical parameters used (e.g., pattern speed), and our findings that Diskfit typically provides good fits to perturbed kinematics of other similar datasets, even in the absence of a bar. For our linear epicyclic perturbation modeling (in the presence of a bar or m = 2 mode), which provides a more physically-based model compared to Diskfit, we used two input rotation models, one based on the best-fit gas rotation model and one on the best-fit stellar rotation model. In both cases we varied the bar pattern speed and damping factor. While the perturbed velocity fields show the characteristic butterfly pattern expected from bar perturbations, very high pattern speeds (≳300 km s⁻¹ kpc⁻¹) are required to cause a resonance close to r ≈ 1″ and thus explain the observed high-velocity features along the minor axis. We note however that if the true distance to NGC 1566 is greater, the problem of high pattern speeds is mitigated. Even if the resonance radii are matched, the velocity structures seen in the models are not aligned with or as sharp as those observed.

We are unable to definitely prove the existence of streaming inflows based on the kinematics of the molecular or ionized gas. This is in part due to the line profiles being complex with evidence of multiple velocity components, and also due to an asymmetry in the velocity profiles of the near and far sides of the disk, which could be due to disk warping or the non-axisymmetric potential. We present and analyze specific apertures along a spiral arm where the spectral profiles are similar to those expected from our toy streaming inflow model and from these we estimate a potential streaming inflow rate of 0.1 M_⊙ yr⁻¹.

We estimate the molecular mass of the unresolved nuclear outflow (the innermost $0\stackrel{″}{.}2$ aperture) as (3.8 ± 0.4)×10⁵ M_⊙ and its momentum as (−19.3 ± 2.8) × 10⁶ M_⊙ km s⁻¹and (12.1 ± 1.7) × 10⁶ M_⊙ km s⁻¹for the blue and red outflows, respectively. Summing all gas believed to be participating in the nuclear outflow (out to ∼$1\stackrel{″}{.}5$ or ∼72 pc from the nucleus) we find a mass of 2.1 × 10⁷ M_⊙. Given the nuclear velocities, the implied outflow rates from the nuclear $0\stackrel{″}{.}2$ region are 3.7 M_⊙ yr⁻¹ for the blueshifted component and 1.9 M_⊙ yr⁻¹ for the redshifted component. We emphasize that these outflows appear to decelerate within the inner 100 pc and thus the gas is not lost to the galaxy nucleus.

We have used the results of three methods to estimate L_Bol: from the nuclear [O III] luminosity, from the nuclear hard X-ray luminosity, and from a SED fit. The first two methods give consistent results: a mean value of L_Bol = (2.7 ± 1.3) × 10³⁴ W, which implies l_Edd of ∼2.2 × 10⁻⁴, indicating a relatively low-efficiency regime for the SMBH, and an accretion rate of ṁ = (4.8 ± 2.3) × 10⁻⁵ M_⊙ yr⁻¹, significantly lower than the posited nuclear molecular outflow rate.

A direct comparison between the molecular outflow kinetic power and the AGN bolometric luminosity gives a ratio of $\frac{{\dot{E}}_{\mathrm{out}}}{L_{\mathrm{Bol}}}\approx 0.28$, a value significantly higher than the typical values (between 0.01 and 10%) found from previous studies of ionized gas in nearby AGNs. This supports the idea that the ionized gas is a minor fraction of the total gas content.

Acknowledgments

We thank the anonymous referee for the very constructive comments that improved this work. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2012.1.00474.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. Based on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the National Research Council (Canada), CONICYT (Chile), Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina), and Ministério da Ciência, Tecnologia e Inovação (Brazil). This research has made use of the services of the ESO Science Archive Facility. N.N. and R.S. acknowledge support from CONICYT (PIA ACT172033, FONDECYT 1171506, and BASAL AFB-170002). R.A.R. acknowledges support from FAPERGS (project No. 2366-2551/14-0) and CNPq (project No. 470090/2013-8 and 302683/2013-5). We acknowledge the use of the HyperLeda database (http://leda.univ-lyon1.fr). This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. IRAF is distributed by the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. G.O. acknowledges the support provided by CONICYT (Chile) through FONDECYT postdoctoral research grant No. 3170942. R.S. acknowledges support from CONICYT Beca/Nacional-Doctorado 21120516.

References

All Figures

Fig. 1. Wide field images of NGC 1566. Left panel: (ESO-DSS image from the UK Schmidt Telescope) shows the full galaxy, while the right panel (ID:13364, PI:Calzetti. HST image taken with the F555W filter) shows the inner morphology, and highlights the inner spirals arms and the ∼1′ bar in PA ∼ 0. In both panels the fields of view of the datasets used in this work are shown in yellow (GMOS/IFU), green (inner ALMA FOV: 12″ × 12″), and red (full ALMA FOV: 27″ × 27″) squares.

In the text

Fig. 2. Moment maps of the CO J:2-1 emission in NGC 1566. Left to right panels: 0th (integrated flux), 1st (velocity) and 2nd (velocity dispersion) moments. Top row: moment maps derived from the highest spatial resolution maps (over a 12″ × 12″ FOV) to best emphasize the nuclear features. The moment 0 (left) panel has a r.m.s. of 0.17 Jy beam−1 km s−1, and pixels with flux density lower than 0.51 Jy beam−1 km s−1 (3σ) in the moment 0 image were “masked” in all panels of the row by setting them to a value which results in a white color in the panel. Bottom row panels: moments derived from lower spatial resolution (but higher signal-to-noise) maps and show a larger 27″ × 27″ FOV to emphasize the larger scale spiral arms. The moment 0 (left panel) has a r.m.s. of 0.16 Jy beam−1 km s−1, and pixels with flux density lower than 0.48 Jy beam−1 km s−1 (3σ) in the moment 0 images were “masked” in all panels.

In the text

Fig. 3. Top row: moment maps of the [N II] 6583 Å emission line in NGC 1566. From left to right, the panels show the 0th (integrated flux), 1st (velocity) and 2nd (velocity dispersion) moments. The nuclear position (marked with a cross) was derived from the peak flux of the continuum in Gemini-GMOS datacube, and the solid lines indicate the major axis (thick line) and the minor axis (thin line). In the moment 0 map (left panel), the brown contours show the integrated flux (moment 0) of the CO J:2-1 line. Bottom row: from left to right the structure map, and maps of the stellar velocity and stellar velocity dispersion. All maps have N to the top and E to the left (see compass). The structure map was created from a HST F606W filter on which unsharp masking was used to emphasize dust features. It is shown at the same size and orientation as the other panels.

In the text

Fig. 5. Amplitudes of the Fourier components obtained from the Kinemetry analysis of the CO velocity map, as a function of distance from the nucleus. The solid red line represents the coefficient of the cosθ term (pure circular rotation in a disk) and the other solid lines show the additional components (associated to perturbations) in the Fourier decomposition, following the colors specified in the inset. Only odd Fourier components were used. The vertical lines delineate radii at which we have abrupt changes in the fraction of pixels at a given radius which have values in the velocity map. At radii smaller than the left vertical line, this fraction is 1. Between the two vertical lines the fraction drops linearly from 1 to 0.5, and beyond 4″ the fraction is relatively steady at ∼0.3. The horizontal dashed line delineates zero velocity.

In the text

Fig. 4. Comparison of the (projected) rotation models and observed velocities extracted along the major axis, in the central kpc of NGC 1566. The solid body model is shown with a dotted line, our ModC2014 model (see text) is shown with a dashed line, the circular (“gascirc”) model obtained by the Kinemetry fit to the CO J:2-1 velocity field is shown with the solid line, and the Bertola model fit to the stellar velocity field is shown with a dashed-dot line. CO J:2-1 and stellar velocities extracted along the major axis are shown with crosses and stars, respectively. Blue is used for the NE (approaching) side of the galaxy, and red for the SW (receding) side of the galaxy. The zero velocity on the y-axis corresponds to 1485 km s−1; at the nucleus, the velocity of the CO J:2-1 line is offset from this by −12.8 km s−1; a consequence of the asymmetric double-peaked profile of the nucleus, most likely caused by the effects of nuclear outflows and/or bar related perturbations.

In the text

Fig. 6. Top panel: maps of the CO J:2-1 velocity dispersion (left panel) and velocity residuals after subtraction of the ModC2014 rotation model (right panel) for the inner FOV. Contours in both panels show the moment 0 (integrated flux) map of the CO J:2-1 emission, ranging from 0.4 mJy beam−1 km s−1 to 5.4 mJy beam−1 km s−1. Bottom left panel: moment 3 (skewness) map of the CO J:2-1 emission for the larger FOV, shown in color following the color bar: blue colors represent spectra with excess emission towards the blue side of the weighted mean velocity. The CO Moment 0 (total flux) map is overlaid with black contours. Bottom right panel: as in the top right panel but for the larger FOV. Pixels with velocity less than −90 km s−1 in the right panels are shown in white. Contours in bottom panels show the moment 0 (integrated flux) map of the CO J:2-1 emission, ranging from 0.4 mJy beam−1 km s−1 to 8.4 mJy beam−1 km s−1.

In the text

Fig. 7. Nuclear velocity field of NGC 1566 in the CO J:2-1 (top row) and [N II] 6583 Å (bottom row) emission lines, in contours and color following the color bar above each panel. The left column shows the expected radial velocity field from the ModC2014 model (see text), the middle column shows the observed radial velocity, and the right column shows the residual velocities (observed − model). A cyan dashed squares were drawing inside bottom panels for showing the FOV of top panels. In these figures we use a systemic velocity of 1485 km s−1and the nuclear position (marked with a cross) is determined by the 230 GHz continuum emission peak. The major (PA = 45°) and minor axes are shown in solid lines. One arcsecond corresponds to 0. pc.

In the text

	Fig. 8. Velocity and dispersion maps of the [N II] emission line obtained from the two-Gaussian fit version of profit.pro. Top panels: from left to right, the velocity and velocity dispersion maps of the narrow and broad components, respectively. Bottom panels: velocity dispersion maps of the narrow (left panel) and broad (right panel) components.
In the text

Fig. 9. Position-velocity diagrams of the CO J:2-1 emission in NGC 1566 along several PAs are shown in color following the color bar above each panel. The PA of the “slit” over which the pv diagram was extracted is indicated above each panel, as are the PAs corresponding to the major and minor axis of the galaxy, and of the large-scale bar. Negative offsets on the x-axis correspond to the PA listed above the panel, i.e., positive offsets are along the 180° plus the listed PA. The black cross indicates the position of the 230 GHz continuum peak (presumed to be the galaxy center) and the systemic velocity of the CO J:2-1 line (1485 km s−1). To guide the eye horizontal and vertical lines delineate ±2 arcsec from the nucleus and ±50 km s−1 from the systemic velocity. The dashed white and purple lines are the solid body rotation model and the ModC2014 rotation model, respectively. The solid black line shows the expectation of adding our outflow model to the ModC2014 model (see Sect. 3.3). The pv diagrams, were created from a “hybrid” datacube: the inner ∼12″ × 12″ square region centered on the nucleus is taken from a high-resolution (Briggs weighting, Robust = −2) map with spatial resolution $0\stackrel{″}{.}52$ × $0\stackrel{″}{.}35$ and an r.m.s. noise of 0.1 mJy beam−1 per channel (up to 2 mJy beam−1 per channel in channels with strong signal), while the rest of the cube is from a higher signal-to-noise (Briggs weighting, Robust = 2) map with spatial resolution $0\stackrel{″}{.}6$ × $0\stackrel{″}{.}5$ and an r.m.s. noise of 1.2 mJy beam−1 per channel (up to 4 mJy beam−1 per channel in channels with strong signal).

In the text

Fig. 10. Top panel: Diskfit model, with perturbations from a bar at PA 0°, fit to the CO velocity field (left panel) and the resulting residual CO velocity field (observed − model) (right panel). Bottom panel: illustrative bar-perturbed velocity fields resulting from our linearized epicyclic perturbation models (see text) and using λ = 0.2 and Ω = 120km s−1 kpc. The left (right) panel shows the results when setting the intrinsic rotation curve as that from C14 (our Bertola model fit to the stellar velocity field).

In the text

Fig. 11. Bar-perturbed velocity fields obtained with our epicyclic perturbation models applied to our best fit Bertola stellar rotation curve (see text), on varying the bar pattern speed (left to right; Ω in units of km s−1 kpc−1) and the dimensionless damping parameter (top to bottom; λ). The FOV of each image is 12″ × 12″ and major tick marks are shown every 1″. All velocities follow the same color bar shown on the top.

In the text

Fig. 12. Position-velocity diagrams, as in Fig. 9, but here the overplotted lines show the velocity of the corresponding pixel in the CO moment 1 map (dashed black lines; i.e., the flux-weighted average velocity at a given pixel), the best-fit velocity field from Diskfit (solid orange line), and for comparison the λ = 0.2 and Ω = 120 km s−1 kpc model from our linear perturbation analysis (bottom right panel of Fig. 10; solid green line). To better trace the highest velocity components in the CO moment 1 (velocity) map, we used a cutoff of 8 mJy beam−1 per channel (∼4σ) to create the moment 1 map which is overplotted here and used as the input map to Diskfit. The PA of the “slit” along which the pv diagram was extracted is marked above each panel. Bottom right panel: the moment 0 map of CO J:2-1 together with the positions of the slits used to create the pv-diagrams in the upper row panels.

In the text

Fig. 13. CO spectral profiles in apertures along the inner spiral arms. Left panel: location of each circular aperture ($0\stackrel{″}{.}2$ in radius) overlaid on the CO flux map. For reference, the center of the galaxy, major axis, and the near and far sides of the disk, are indicated. Apertures are numbered 0–8 with aperture 0 being the farthest aperture on the E arm, aperture 4 the nuclear one, and aperture 8 the farthest aperture on the W arm. Right panel: extracted CO J:2-1 spectra; thick lines are used for the spectra corresponding to apertures from the W arm (Apertures 5–8), and different colors are used for each spectrum, following the color bar on top of the panel. Symbols with the corresponding color (plotted at the y value of the peak flux density of the spectrum) denote the radial velocities expected in that aperture for our rotation model (hourglasses), our outflows + rotation model (asterisks), and our rotation plus radial streaming inflow model (open circles). In some apertures, adding outflows and/or streaming inflows does not change the predicted radial velocity; this is a result of projection effects and/or the fact that our outflow model has zero velocity beyond ∼2″ from the nucleus.

In the text

	Fig. 14. Same as Fig. 13, but for apertures oriented along the bar. Apertures are numbered 0 (southernmost) to 4 (northernmost), and apertures here are each $0\stackrel{″}{.}5$ in radius.
In the text

	Fig. 15. CO J:2-1 spectrum (black) in a circular aperture centered on the nucleus with radius $0\stackrel{″}{.}2$. The spectrum was divided by eye into three parts: the pure outflow components (blue and red) and an intermediate velocity region (green; see text) between the blue and red outflow peaks.
In the text

	Fig. 16. CO J:2-1 spectra extracted from apertures 1 and 2 of Fig. 13 (the E inner spiral arm) are shown in black. The red overlay marks the velocities over which emission comes from gas potentially participating in streaming inflows.
In the text