© ESO 2019

1. Introduction

The detection of the hidden broad-line region of NGC 1068 (Antonucci & Miller 1985), the archetypical type 2 active galactic nucleus (AGN), supported the idea of the so-called unified model. A dusty molecular torus obscures the view of the central engine in type 2 AGN along directions near the equatorial plane, while near polar directions, it allows the detection of the broad-line region, thus classifying the AGN as a type 1. This torus was proposed to be both geometrically and optically thick (e.g., Antonucci 1993). To explain the observed properties of NGC 1068 and other active galaxies, the theoretical models initially favored a relatively compact torus (Pier & Krolik 1992; Nenkova et al. 2002, 2008a,b). Subsequently, it was suggested that the torus is part of a larger and more diffuse structure that extends for up to 100 pc (Pier & Krolik 1993; Granato & Danese 1994; Maiolino & Rieke 1995; Efstathiou & Rowan-Robinson 1995; Young et al. 1996). However, only the inner warm parts of the torus would be responsible for the near- and mid-infrared emission (Schartmann et al. 2008). Additionally, most local Seyfert galaxies with sufficiently high signal-to-noise ratio mid-infrared interferometric observations show emission along the polar direction on parsec scales (Tristram et al. 2009; Hönig et al. 2013; López-Gonzaga et al. 2016). Theoretical arguments indicate that this phenomenon could be associated with an outflowing torus, as already captured by dusty disk-wind models (Elitzur & Shlosman 2006; Hönig & Kishimoto 2017) and radiation-driven outflow models (Wada 2012; Wada et al. 2016). Other alternatives to the optically and geometrically thick torus, including nuclear warped disks (see Jud et al. 2017, for Circinus), are discussed in detail by Lawrence & Elvis (2010). See also Ramos Almeida & Ricci (2017) for a recent review.

The excellent high angular resolution achieved with the Atacama Large Millimeter Array (ALMA) is now revealing that relatively large (diameters of up to 20−50 pc) molecular tori or disks appear to be common in the few local Seyferts and low-luminosity AGN that have been observed so far (García-Burillo et al. 2014, 2016; Gallimore et al. 2016; Imanishi et al. 2018; Alonso-Herrero et al. 2018; Izumi et al. 2017, 2018; Combes et al. 2019). Moreover, these nuclear tori or disks appear to be decoupled from their host galaxies both morphologically and kinematically. This is in good agreement with previous results that found random orientations of the radio jets or ionization cones of local AGN with respect to their host galaxies (see, e.g., Nagar & Wilson 1999; Kinney et al. 2000; Fischer et al. 2013). Studying the mechanisms responsible for this decoupling as well as the physical scales where they operate is thus crucial for proceeding in our understanding of the general question of the fueling and nuclear obscuration of active galaxies.

This is the second paper of a series (first paper, Alonso-Herrero et al. 2018) based on several ALMA programs to observe the cold molecular gas in the nuclear and circumnuclear regions of an ultra-hard X-ray selected sample of nearby Seyfert galaxies. Our parent sample is drawn from the X-ray Swift/BAT 14 − 195 keV all-sky 70-month catalog (Baumgartner et al. 2013). Our main goal is to understand the connections between the cold and hot molecular gas, the AGN torus, and nuclear and circumnuclear star formation activity in local active galaxies. In this paper we present a detailed study of the cold molecular gas emission of the nearby Seyfert 1.5 galaxy NGC 3227. We use a distance of 14.5 ± 0.6 Mpc as derived from AGN lags (Yoshii et al. 2014) and thus 1″ ∼ 73 pc. The intrinsic 14 − 195 keV X-ray luminosity is approximately 2 × 10⁴² erg s⁻¹ (Ricci et al. 2017) for our assumed distance. This galaxy is interacting with NGC 3226 and presents evidence for large-scale streaming motions detected in neutral hydrogen (Mundell et al. 1995a). It is also located in a group of 13−14 galaxies (Davies et al. 2014). Finally, NGC 3227 is classified as SAB(s)a pec in the Third Reference Catalogue of Bright Galaxies (de Vaucouleurs et al. 1991). However, Mulchaey et al. (1997) modeled near-infrared observations of this galaxy, which revealed a large -scale bar.

Schinnerer et al. (2000) used the Plateau de Bure interferometer (PdBI) to obtain CO(1–0) and CO(2–1) observations of NGC 3227 with angular resolutions in the ∼0.5 − 1.5″ range. They detected apparent counterrotation in the nuclear region of this galaxy and interpreted this kinematic signature as evidence of the warping of a thin inner molecular disk starting at an approximately outer radius of 1″. Davies et al. (2006), on the other hand, based on the high-velocity dispersion of the hot molecular gas traced by the H₂ 2.12 μm 1–0 S(1) emission line, argued that the nuclear kinematics could be explained with a nuclear thick disk with turbulent motions. In this paper we present ALMA Cycle 4 observations of the molecular gas CO(2–1) and CO(3–2) transitions and associated continuum of the nuclear region of NGC 3227, obtained with angular resolutions ranging from 0.21″ to 0.085″. The paper is organized as follows. Section 2 presents the ALMA observations and ancillary data used in this paper. In Sects. 3 and 4 we discuss the properties of the (sub-) millimeter (submm) continuum and molecular gas CO(2–1) and CO(3–2). We model the CO(2–1) molecular gas kinematics in Sect. 5. In Sect. 6 we derive the cold molecular gas mass and column density in the nuclear region of NGC 3227. Finally, in Sect. 7 we present our conclusions.

2. Observations and data reduction

2.1. ALMA observations in bands 6 and 7

We obtained band 6 ALMA observations of NGC 3227 between March and September 2017 using the 12 m array in two configurations (compact with baselines between 15 and 460 m and extended with baselines between 20 and 3700 m) and the Atacama Compact Array (ACA) with baselines between 9 and 50 m. These observations were part of the project 2016.1.00254.S (PI: A. Alonso-Herrero). We defined two spectral windows of 1.875 GHz bandwidth (with 3.9 MHz ∼ 5 km s⁻¹ channels), at the observed frequency of the CO(2–1) transition (229.8 GHz) and the other at an observed frequency of ∼231 GHz (1.3 mm) to measure the submm continuum. We also used band 7 ALMA data of NGC 3227 observed between August and September 2017 with the 12 m array in a single configuration with baselines between 21 and 3700 m observed for the 2016.1.01236.S (PI: M. Vestergaard) project. A spectral window was centered on the CO(3–2) transition (344.4 GHz) and another three spectral windows at 342.6, 354.6, and 356.5 GHz to measure the continuum (approximately 850 μm).

We calibrated and imaged the data using the ALMA reduction software CASA (v.5.1, McMullin et al. 2007). For the 12 m array data (Band 6 and 7), either J0854+2006 or J1058+0133 were used for both the bandpass and amplitude calibrations, and J1025+1253 was used for the phase calibration. For the band 6 ACA data, J1058+0133 was used for the bandpass calibration, Ganymede for the flux calibration, and J1041+0610 for the phase calibration. The Ganymede flux was estimated using the Butler-JPL-Horizons 2012 model. The consistency of the flux calibration was verified by comparing the observed amplitudes as a function of the uv distance for all the configurations. Because the amplitudes of the overlapping baselines are similar, we can conclude that the flux calibration is consistent. The J0854 and J1058 fluxes are 2.5 and 3.1 Jy at 229 GHz. Finally, we note that we used the same phase calibrator for bands 6 and 7 except for the ACA band 6 data (J1041+0610). However, the ACA data are not expected to significantly affect the emission that is detected on the scales discussed in this work, and therefore all the data are referenced to the same phase calibrator.

We subtracted the corresponding continua from the CO(2–1) and CO(3–2) spectral windows directly in the visibility data by fitting the continuum emission with a constant in the line-free channels. Then, we combined and cleaned the data using the CASA CLEAN task. The output frequency reference frame was set to the kinematic local standard of rest (LSRK). For the band 6 continuum images, we used a continuum spectral window that we cleaned with the CASA CLEAN task using Briggs weighting (Briggs 1995), with robustness parameters of b = 0.5 and b = −0.5. These result in a lower angular resolution and higher sensitivity, and in a higher angular resolution and slightly decreased sensitivity, respectively. We summarize the beam properties and sensitivities of the continuum maps of bands 6 and 7 in Table 1. For the band 7 data, we combined the line-free channels of the four spectral windows also using the CASA CLEAN task and a natural weight (b = 2). For the CO(2–1) and CO(3–2) transitions we produced cleaned data cubes with weights b = 0.5 and natural, respectively, and the resulting beam sizes and sensitivities are listed in Table 2. These two tables show that the angular resolutions of our ALMA observations lie between 0.085 and 0.21″. For the assumed distance to NGC 3227, they probe physical resolutions in the range of 7−15 pc.

From the data cubes of bands 6 and 7, we produced maps of the CO(2–1) and CO(3–2) integrated intensity, mean-velocity field, and velocity dispersion. For this we used the GILDAS¹ MOMENT task and set the threshold limit to 3σ and 5σ. In Fig. 1 (left panel) we show the integrated CO(2–1) molecular line map created from the band 6 data cube with b = 0.5 and the 3σ detection threshold. We show an approximate field of view (FoV) of 10″ × 10″, which corresponds to a circumnuclear region of 730 pc × 730 pc. In this and the following figures, the coordinates are in the International Coordinate Reference System (ICRS) frame.

Fig. 1. Left panel: ALMA CO(2–1) integrated molecular line map (zeroth-order moment map) created using a 3σ threshold. The image displays the central 10″ × 10″ (730 pc × 730 pc) and was created from the b = 0.5 weighted data cube. We show the image and contours with a linear scale. The CO(2–1) contours are 0.0925, 0.515, 0.9375, 1.36, 1.7825, 2.205, 2.6275, and 3.05 Jy km s−1 beam−1. The star indicates the AGN location, based on the peak of the unresolved 1.3 mm continuum emission (see text). Right panel: HST/ACS pseudo-Hα+[N II] line map (highlighted in light colors) plus extinction map (highlighted in dark colors) constructed as an F658N/F814W map (in arbitrary units) and showing the same FoV as the ALMA CO(2–1) image. The broadband F814W image is saturated at the nuclear position.

2.2. Ancillary archival HST images

We downloaded from the Hubble Legacy Archive (HLA) images taken with the Advanced Camera for Surveys (ACS) on board the Hubble Space Telescope (HST) using the narrowband filter F658N, which contains the Hα+[N II] emission lines, and the broad-band filter F814W (proposal ID: 9293). The HLA images are fully reduced and drizzled to a pixel size of 0.05″. We adjusted the HST/ACS astrometry for NGC 3227 by fixing the coordinates of the bright nuclear point source to that of the peak of the unresolved ALMA band 6 continuum at 1.3 mm (see Sect. 3.1). To trace the emission from the NGC 3227 ionization cone and circumnuclear H II regions (see Gonzalez Delgado & Perez 1997) at high angular resolution, we constructed a pseudo-Hα+[N II] line map by dividing the F658N image by that of the nearby broadband continuum traced by the F814W image (Fig. 1, right panel). This map also highlights the regions with high extinction. They are seen as dust lanes toward the west, southwest, and southeast of the nucleus (marked in the figure), which are due mostly to material on the near side of the galaxy (see Martini et al. 2003; Davies et al. 2014).

In Fig. 1 we show the F658N/F814W map for the same FoV as the ALMA CO(2–1) integrated line emission map. As inferred from optical imaging, line ratios, and emission line kinematics of this galaxy (Mundell et al. 1995b; Gonzalez Delgado & Perez 1997), the line emission to the north, northwest, and northeast is mostly photoionized by the AGN, and it likely mostly traces one side of the ionization cone of this galaxy (Fischer et al. 2013), whereas the emission at ∼4″ to the southwest of the AGN is mostly due to star formation activity in the disk of the galaxy. Further indication for AGN photoionization to the north and northeast is the detection of outflowing ionized gas at PA ∼ 10°, which reaches velocities of up to 900 km s⁻¹ (see Barbosa et al. 2009). Moreover, modeling of the ionized gas kinematics finds a position angle (measured east from the north) for the ionization cone of PA_cone = 30° (Fischer et al. 2013). This value is also in good agreement with the value derived from optical polarization, assuming that the polarization is due to polar scattering (Smith et al. 2004).

3. Nuclear continuum emission

3.1. Morphology

In Fig. 2 we show the images of the 1.3 mm and 850 μm continua of NGC 3227 for the nuclear 1″ × 1″ region (73 pc × 73 pc). The continuum at both wavelengths appears clearly extended. The 1.3 mm continuum presents a complex morphology with an unresolved source and probably several extended components (see below). The 1.3 mm continuum map derived from the b = 0.5 weighted data cube shows an extended component along the northeast-southwest direction, that is, at a PA ≃ 30 − 40°. This direction coincides with the axis of the projected Hα+[N II] emission (see the right panel of Fig. 1), which mostly traces emission on one side of the cone in the ionization cone (see, e.g., Mundell et al. 1995b; Gonzalez Delgado & Perez 1997; Barbosa et al. 2009; Fischer et al. 2013). The 1.3 mm continuum map derived with the b = −0.5 weight, which has a higher angular resolution (see Table 1 and middle panel of Fig. 2), more clearly shows that in addition to the emission in the direction of the ionization cone, another component lies nearly perpendicular to it (see below).

Fig. 2. ALMA continuum images in band 6 at ∼1.3 mm from the b = 0.5 weighted data cube (left) and b = −0.5 weighted data cube (middle) and in band 7 at ∼850 μm from the natural weight data cube (right). The FoV is 1″ × 1″, which corresponds to approximately the central 73 pc × 73 pc region. In all images the first contour corresponds to the 2.5σ level and the color bars indicate the flux in units of 10−4 Jy beam−1. The hatched ellipses show the corresponding beam size and orientation in each image (see Table 1), and the star indicates the approximate location of the peak of the unresolved 1.3 mm continuum emission (see text).

For simplicity, we modeled the 1.3 mm continuum in the uv plane using only two components, namely a point source and a Gaussian source, to encompass all the extended emission. We left the positions of the two sources free and only fixed the width of the extended component to FWHM = 0.35″ due to the limited signal-to-noise ratio of this component. For the unresolved component, we obtained the following coordinates (ICRS frame): RA(ICRS) = 10^h23^m30.574^s and Dec(ICRS) = +19°51′54.278″. These coordinates are within the errors coincident with source C detected at 1.7 and 5 GHz by Bontempi et al. (2012) using radio observations obtained with the high angular resolution European VLBI Network (EVN). Because of the high brightness temperature and spectral index of source C, these authors identified it as the position of the AGN of NGC 3227. The 1.3 mm extended component does not peak at exactly the same position, but it appears to be offset by approximately 0.07″ to the southeast of the AGN position. This offset appears to be real because it is a factor of a few above the positional uncertainties of the phase calibrator used during the data reduction. The derived fluxes at 1.3 mm are 0.51 ± 0.15 mJy and 1.92 ± 0.21 mJy for the unresolved and resolved components, respectively.

The 850 μm continuum map (Fig. 2, right panel) shows that most of the nuclear emission at this wavelength has a disk-like morphology extending for about 0.6″ (∼45 pc) and oriented at PA ∼ −20° to −30°, which is perpendicular to the axis of the ionization cone. There is also fainter emission in the direction of the ionization cone (see also the top panel of Fig. 3). The 850 μm continuum morphology shows a remarkable resemblance to the SINFONI Brγ narrow-component map presented by Davies et al. (2006; their Fig. 3, with a similar FoV as our continuum maps). These authors interpreted this Brγ emission as mostly originating in star formation activity, and it is probably associated with the nuclear stellar disk. There is also evidence for nuclear ongoing or recent star formation activity from the detection of polycyclic aromatic hydrocarbon emission on scales of ∼1″ − 0.5″ (Imanishi 2002; Rodríguez-Ardila & Viegas 2003; Davies et al. 2007; Esquej et al. 2014; Alonso-Herrero et al. 2016). Moreover, Davies et al. (2006) also suggested that the MERLIN 6 cm radio continuum emission (Mundell et al. 1995b), which has an extended component in approximately the same direction as the nuclear disk, might be due to the superposition of many supernova remnants rather than due to a radio jet (see also Chapman et al. 2000).

Fig. 3. Top panel, color and white contours: map of the 850 μm continuum emission (in arbitrary units) degraded to the beam size and orientation of the 1.3 mm continuum map with b = 0.5 weight (see Table 1). This was necessary to construct the spectral index map. The map is shown for the nuclear 1″ × 1″ region. We mark the approximate orientations of the ionization cone axis and the nuclear star forming disk. Bottom panel: map of the spectral index α as derived from the ALMA continuum observations in bands 6 and 7 using only pixels detected at > 2.5σ. The FoV is the same as in the top panel. The map has the same angular resolution as the band 6 map with the b = 0.5 weight (see Table 1). The color bar indicates the values of α (see text). The contours are the band 6 continuum emission at 1.3 mm as in the left panel of Fig. 2, and the star marks the AGN position.

3.2. Spectral index

To investigate the nuclear (sub-) millimeter continuum emission of NGC 3227 further, we constructed a map of the 1.3 mm to 850 μm spectral index, defined as f_ν ∝ ν^α, with the continuum images of ALMA bands 6 and 7. We used GILDAS to convolve the 850 μm continuum map to the beam size and orientation of the 1.3 mm continuum map for the b = 0.5 weight (0.196″ × 0.151″ at PA_beam = 48°). We show the result in the bottom panel of Fig. 3. For the spectral index map we only used pixels where both continua were detected at a 2.5σ level or higher. With the adopted definition of the spectral index, positive values in the range α = +1 to α = +2 would indicate a thermal dust component (from the dusty molecular torus or dust heated by star formation or both). Pure synchrotron radiation (from a jet or supernova remnants) has α = −0.8 (for non-beamed emission), and thermal free-free emission (associated for instance with H II regions) takes values of α ≃ 0 for optically thin emission and α ≃ +2 for optically thick emission. More complicated scenarios, which have been discussed for NGC 1068, include synchrotron emission with free-free absorption and electron-scattered synchrotron emission (Hönig et al. 2008; Krips et al. 2011; Pasetto et al. 2019). These processes result in positive spectral indices. Finally, optically thick synchrotron emission with α ∼ 0 can also be relevant to low-luminosity AGN (Koljonen et al. 2015).

In Fig. 3 (bottom panel) we show the spectral index map of the nuclear region of NGC 3227. The spectral index map displays a range of values. These variations are real because they take place on scales larger than the map beam. This suggests that the (sub-) millimeter emission in this region is a combination of different mechanisms. To the south and southwest of the AGN position out to projected distances of 0.2″, the spectral index appears to be mostly negative. This might indicate that synchrotron radiation dominates in this region. This direction coincides with that of the south part of the cone, which is obscured in the optical due to extinction from the host galaxy², as inferred from the modeling of the optical ionized gas kinematics (Fischer et al. 2013).

To the north and northeast of the AGN, the observed spectral index could reflect a combination of thermal processes and synchrotron emission. One possible interpretation for the extended 1.3 mm emission at PA ∼ 30° to the northeast would be a combination of polar dust emission and emission produced in a radio jet. In the mid-infrared, this galaxy shows extended 8 − 12 μm emission towards the north, northeast, and south of the AGN (scales of approximately 2″, see Alonso-Herrero et al. 2016; Asmus et al. 2016; García-Bernete et al. 2016), which might be tracing dust emission in the ionization cone. When we interpret the nonthermal part of the extended 1.3 mm emission as a radio jet, the misalignment with the ionization cone emission and radio would be small. It is worth mentioning that Mundell et al. (1995b) noted that they also detected 18 cm radio emission perpendicular to PA_radio = −10°, which could be associated with the ALMA extended emission at 1.3 mm. This extended radio emission at a position angle similar to that of the ionization cone is also visible at 6 cm (see Fig. 9 of Davies et al. 2006).

Along the nuclear disk at PA ∼ −20°. and more importantly, toward the northwest of the AGN, the spectral index takes positive values. This could be associated with the ongoing or recent star formation activity in the nuclear disk (Davies et al. 2006, 2007) via thermal processes such as dust emission or free-free emission due to H II regions or both.

3.3. Nuclear radio to submm SED

In addition to the ALMA spectral index map, we also compiled radio to (sub-) millimeter small-aperture measurements to construct the nuclear spectral energy distribution (SED) shown in Fig. 4. For the radio part we used the small-aperture (∼0.1 − 0.2″) MERLIN 18 cm and 6 cm data points of the southern source of Mundell et al. (1995b), which we associate with the AGN position. From the ALMA 1.3 mm and 850 μm continuum maps, we measured the flux densities through a square aperture 0.2″ (∼15 pc) on a side. We also measured the ALMA continuum flux densities for a larger square aperture 0.5″ on a side to approximately encompass the emission from the nuclear disk detected at 850 μm (see Fig. 2). Table 3 lists the ALMA continuum measurements. We compare them with the Very Large Array (VLA) 5 GHz total flux density of Nagar et al. (2005) and the 3 mm flux density measured within approximately 1″ using the PdBI by Sani et al. (2012). To identify the AGN, we show the EVN radio flux density for source C of Bontempi et al. (2012), which has a size of a few milliarcseconds. We also show the fit to the nuclear infrared emission of NGC 3227 (typical angular resolutions 0.2 − 0.3″) using the Nenkova et al. (2008a); Nenkova et al. (2008b) clumpy torus models derived by García-Bernete et al. (2019).

Fig. 4. Nuclear continuum SED from centimeter wavelengths to the infrared of NGC 3227. The small-aperture (∼0.1 − 0.2″) MERLIN data points are for the southern source of Mundell et al. (1995b). The large-aperture (> 0.2″ to ≤1″) VLA and PdBI data points are from Nagar et al. (2005) and Sani et al. (2012), respectively. All the ALMA continuum data points are from this work. The very small aperture (a few milli-arcseconds) points are the EVN radio fluxes from Bontempi et al. (2012). The fit (black dashed line, with the ±1σ uncertainty plotted with the shaded gray region) is made to the near- and mid-infrared emission using the Nenkova et al. (2008a,b) clumpy torus models (see García-Bernete et al. 2019, for more details). We also plot the observed mid-infrared data point for reference although the fitted data also included near-infrared photometry and mid-infrared spectroscopy. The magenta dashed line is the fit to the MERLIN radio points extrapolated to the ALMA wavelengths.

The derived spectral index fitted to the MERLIN nuclear values, α_MERLIN = −1.1, is similar (see also Mundell et al. 1995b) to that derived between 6 and 20 cm by Edelson (1987) for NGC 3227 and other Seyfert galaxies. The simplest interpretation is that it mostly reflects synchrotron emission at these frequencies. When we extrapolate this fit to the frequencies of ALMA bands 6 and 7 (dashed magenta line in Fig. 4), the contribution from optically thin synchrotron emission to the ALMA bands is small. However, the core emission might be self-absorbed. It is also possible that the low-frequency spectral index of −1.1 may not hold all the way to the millimeter wavelengths (Behar et al. 2018). Similarly, the predicted emission from the fitted Nenkova et al. (2008a); Nenkova et al. (2008b) clumpy torus model at 1.3 mm and 850 μm would only account for a small fraction of the observed emission at these frequencies, even for the 0.2″, aperture. This suggests that even at 850 μm within the inner ∼14 pc of NGC 3227, there is a significant contribution from star formation that may produce thermal and nonthermal emission in this region. Davies et al. (2006) estimated that approximately 75% of the narrow component Brγ emission might originate in star formation in this nuclear region (see also Davies et al. 2007). Alternatively, there could be a more extended cold disk component that is still heated by the AGN and is is not captured by the clumpy torus models.

4. CO molecular gas emission

4.1. Circumnuclear region

Figure 1 (left panel) shows the integrated CO(2–1) intensity map (zeroth-order moment map) covering a FoV of 10″ × 10″ with a beam size of 0.214″ × 0.161″ at PA_beam = 42.0°. As previously shown by Schinnerer et al. (2000) using PdBI observations at lower angular resolution, the circumnuclear cold molecular gas shows an asymmetric ring-like morphology with a projected diameter of approximately 4 − 5″ (∼350 pc). This circumnuclear ring is also coincident with a stellar ring with low-velocity dispersion, which Barbosa et al. (2006) interpreted as a ring with recent star formation activity. The circumnuclear ring is connected with the host galaxy disk through a large-scale molecular bar that extends for 10−15″, to the northwest and southeast of the galaxy nucleus, with the former part being brighter and in the direction of the companion galaxy (Meixner et al. 1990; Schinnerer et al. 2000).

Our ALMA CO(2–1) observations further resolve the circumnuclear ring into several clumps, none of which coincides with the AGN location derived from the peak of the unresolved 1.3 mm continuum (see Sect. 3.1 and below). The comparison of the CO(2–1) map with the HST F658N/F814W map in the central 10″ × 10″ almost suggests an anti-correspondence between the ionized gas (bright regions in this map) and the molecular gas, except for the bright H II region located approximately 4″ southwest of the AGN. This may be explained by the assumption that most of the cold molecular gas is in the disk of the galaxy (see Sect. 5), while the ionized emission in the circumnuclear region to the northeast traces outflowing material in the ionization cone (Mundell et al. 1995b; Barbosa et al. 2009; Fischer et al. 2013). The regions of high extinction to the west and southwest on the near side of the galaxy, seen as dark lanes in the F658N/F814W map and in the color maps presented by Chapman et al. (2000), are all well traced by the CO(2–1) emission. This coincidence between molecular gas and dust lanes is expected and has been observed in other Seyfert galaxies (García-Burillo et al. 2005; Combes et al. 2013, 2014; Alonso-Herrero et al. 2018; Izumi et al. 2018).

4.2. Nuclear region

Figure 5 zooms in on the nuclear 1″ × 1″ region of NGC 3227. The left panels display in color the CO(2–1) and CO(3–2) integrated line emission. These maps show a complex morphology with several emission peaks. The higher angular resolution CO(3–2) map reveals that none of these peaks is associated with the AGN position. The peaks of the CO(2–1) emission appear to be distributed along PA_gas = −40°, extending for about 1.5″ (see also Fig. 1). The morphology of the extended CO(2–1) emission is similar to that seen in the 850 μm map when degraded to the same resolution (top panel of Fig. 3). The CO(3–2) morphology in the nuclear region is similar to that of the hot molecular gas traced by the VLT/SINFONI adaptive optics near-infrared H₂ 1–0 S(1) line at 2.122 μm (Davies et al. 2006). We reproduce this map in Fig. 6 (left). It has a FWHM angular resolution of 0.085″, which is similar to that of the ALMA CO(3–2) map. This indicates that near the AGN position, these CO transitions trace warmer gas than in other environments (e.g., star-forming regions), as was also found in NGC 1068 by Viti et al. (2014).

Fig. 5. Zoom-in on the CO(2–1) and CO(3–2) emission in the nuclear 1″ × 1″ region in the top and bottom panels, respectively. These moment maps were created with the GILDAS MOMENT task using a 5σ detection threshold and show the integrated emission (zeroth-order moment map) in the left panels and the mean velocity field (first-moment map) in the right panels. The hatched ellipses show the corresponding beam size and orientation of each image (see Table 2). The contours in all panels are the ALMA continuum emission at 1.3 mm as in the left panel of Fig. 2, and the star marks the AGN position.

The CO(3–2) intensity map resolves the inner 0.5″ (∼37 pc) emission into two main regions, with the northwest clump being brighter of the two. In projection, the AGN is located approximately in between these two CO(3–2) nuclear peaks. We note that even though the AGN does not coincide with any of the CO(3–2) clumps, there is molecular gas emission at the AGN location, as we show in Sect. 6. This nuclear molecular gas emission is also elongated, similarly to that of the nuclear disk 850 μm and 1.3 mm continuum emission. Furthermore, the superposition of the 1.3 mm emission contours on the CO(3–2) maps (Fig. 5) reveals that the faint 1.3 mm continuum emission also traces the innermost molecular gas clumps. There is also CO(3–2) emission in the direction of the polar dust and radio jet to the northeast of the AGN. However, to the southwest of the AGN, there is a region in which some of the gas might have been evacuated. This region shows evidence of shocks (Schönell et al. 2019). Based on the high-velocity dispersion of the hot molecular gas in this region, Davies et al. (2006) interpreted the hot molecular gas emission as due to a highly turbulent rotating disk (see also Hicks et al. 2009).

5. CO(2–1) and CO(3–2) kinematics

In this section we study the molecular gas kinematics of the nuclear and circumnuclear regions of NGC 3227 (central 10″ × 10″) using our ALMA CO(2–1) and CO(3–2) observations. Figure 7 (left panel) shows that in the central ∼4″ × 4″ (a few hundred parsecs) the molecular gas follows circular motions. The CO(2–1) velocity field in this figure is noticeably similar to that of the 1–0 S(1) H₂ line (Davies et al. 2014). However, at radial distances greater than 2 − 3″, the CO(2–1) and H₂ velocity fields also show clear deviations from circular motions (see Sect. 5.1). The stellar kinematics, on the other hand, shows relatively well-ordered motions on these scales (Barbosa et al. 2006; Davies et al. 2014), although with high-velocity dispersion (Davies et al. 2006) in the nuclear region. Within the inner 1″, the CO(2–1) and CO(3–2) transitions show very complex kinematics (see the mean velocity fields in the right panels of Fig. 5). The motions do not appear to follow the large-scale velocity field, except at the edges of the FoV. To the southeast (approximately 0.2″) of the AGN position and along the direction of the nuclear disk lies a region with apparent counterrotation. In the direction of the ionization cone axis there are also clear deviations from circular motions. The hot molecular gas kinematics (see the VLT/SINFONI map in the right panel of Fig. 6) traced with the 1–0 S(1) H₂ line (Davies et al. 2006) again shows a remarkable resemblance with the CO(2–1) and CO(3–2) kinematics. Davies and collaborators also noted the high-velocity dispersions and interpreted the observations as due to a nuclear thick disk.

Fig. 6. VLT/SINFONI H2 1–0 S(1) intensity (left) and velocity (right) maps from Davies et al. (2006) showing a FoV similar to that of Fig. 5. The maps were derived by fitting a Gaussian to the emission line. The black dot marks the approximate location of the peak of the near-infrared continuum. The velocities are relative to this peak. To facilitate the comparison with the CO(2–1) and CO(3-2) emission, we show the same velocity range as in the right panels of Fig. 5.

To understand the kinematics and the strong deviations from circular motions, we modeled the CO(2–1) data cube with ^3DBAROLO (Di Teodoro & Fraternali 2015). This code was developed to fit 3D tilted-ring models to emission-line data cubes using simple rotating disks. We first started by fitting a simple ^3DBAROLO model of a rotating disk allowing the following parameters to vary: the systemic velocity (v_sys), the disk inclination (i_disk) and PA³ of the kinematic major axis (PA_MAJ), the circular velocity, and the velocity dispersion (σ_gas). For the kinematic center we used the position of the AGN as that of the unresolved component from the two-component fit to the 1.3 mm continuum emission (see Sect. 3.1). We first fit the central few arcseconds where clear circular motions are observed and derived an average value of v_sys = 1140 km s⁻¹. With the systematic velocity fixed, we fitted the central 10″ × 10″ and derived average values of the major axis PA and inclination of the host galaxy disk of PA_MAJ = 152° and i_disk = 52°, respectively. These are in good agreement with those inferred from the stellar kinematics over similar physical scales (Barbosa et al. 2006; Davies et al. 2006, 2014; Hicks et al. 2009). The next step was to fix the values of the galaxy disk PA and inclination and the systemic velocity and fit only the velocity dispersion and circular velocity. The resulting ^3DBAROLO model together with the CO(2–1) mean velocity field (first-order moment map) are presented in Fig. 7. In the following subsections we use this simple rotating disk model as our base model to first search for departures from circular motions and then to add nuclear components to try to reproduce the observed kinematics.

Fig. 7. 3DBAROLO maps of the observed CO(2–1) mean velocity field (first-order moment, left panel), the fitted model of a rotating disk with fixed PA, inclination (see Table 4 and text), systemic velocity, and kinematic center (middle panel), and the residual mean-velocity field (right panel). The FoV is 10″ × 10″, as in Fig. 1. The color bars indicate line-of-sight velocities in units of km s−1. The contours for the velocities in the left and middle panels are in steps of 50 km s−1. Right panel: the contours are the CO(2–1) intensity (zeroth-order moment map) on a linear scale.

5.1. Streaming motions

We constructed a CO(2–1) residual mean velocity field map (right panel of Fig. 7) by subtracting the ^3DBAROLO rotating disk model from the observed mean velocity field. On circumnuclear scales, the model reproduces the mean velocity field, except along the minor kinematic axis. There are strong velocity deviations from circular motions at approximately projected r >  2 − 3″ along this axis, as well as in the nuclear region (see below). The velocity residuals, which are of about 100 km s⁻¹ or more, are positive to the northwest and west (projected radial distances greater than ∼200 pc) and negative to the southeast (projected radial distances greater than ∼300 pc). Davies et al. (2014) observed similar velocity residuals on scales of a few arcseconds from the AGN when they modeled the kinematics of the near-infrared H₂ line at 2.12 μm with a rotating disk.

We can also compare the observed CO(2–1) position-velocity (p − v) diagrams along the major and minor kinematic axes of the host galaxy disk (shown in Fig. 8 in gray scale and black contours) with the fitted ^3DBAROLO model of a rotating disk (red contours, top panel). Along the kinematic major axis, the ^3DBAROLO model with the rotating disk fits the rotation curve of the host galaxy disk for projected radial distances ∼1 − 5″ fairly well. However, it fails to reproduce it at radial distances from the nucleus of less than approximately 1″, where we find emission from forbidden velocities for a simple rotating disk. Along the kinematic minor axis of the host galaxy disk, the ^3DBAROLO rotating disk model again only reproduces the observations outside the nuclear regions. However, at radial distances greater than 1″, the non-circular motions to the southwest and northeast parts of the galaxy are clearly visible, where we identified the strong mean-velocity residuals in Fig. 7. Along this position angle the inflow appears to be stronger to the southwest of the AGN.

Fig. 8. Grey-scale and black contours: observed CO(2–1) p − v diagrams along the kinematic major axis (left panels) and minor axis (right panels) of the host galaxy disk. We only show fluxes above 5σ. In all the panels the red contours represent the 3DBAROLO models as follows (see also Table 4): the first row is a rotating disk, the second row is a rotating disk with a slight nuclear warp, the third row is a rotating disk + a warp model consistent with the orientation of the ionization cone, the fourth row is a rotating disk with a nuclear radial component, and the fifth row is a rotating disk with a nuclear warp + radial component. The model contour levels are identical to the contour level of the observations. The vertical dotted black line marks the AGN position, and the horizontal dotted black line shows the derived systemic velocity. Top panels: we also highlight nuclear regions of forbidden velocities for a purely rotating disk. Along the minor axis the blue boxes indicate the nuclear outflow and the magenta ones regions with streaming motions due to the large-scale stellar bar.

Fig. 8. continued.

Figure 9 shows the same observed p − v diagrams along the kinematic major and minor axes of the host galaxy disk for the CO(3–2) transition at higher angular resolution (see Table 2). Qualitatively, we observe the same evidence for non-circular motions along the kinematic minor axis of the host galaxy, especially to the southwest of the AGN and as close as 0.5−0.7″ (projected distances 35−50 pc) from the AGN location (Fig. 9, right panel). The clear kinematic decoupling of the nuclear emission in the p − v diagram along the major kinematic axis is also obvious at projected radial distances smaller than ∼0.7″.

Fig. 9. Observed CO(3–2) p − v diagrams along the kinematic major axis (left) and minor axis (right) of the host galaxy disk. We only show fluxes above 3σ. The black contours are fluxes in a linear scale. The black dotted vertical line marks the AGN position, whereas the horizontal red lines mark the approximate velocity virial range of the simple rotating disk model (see top panel of Fig. 8).

As discussed in the Introduction, NGC 3227 presents a large -scale stellar bar (Mulchaey et al. 1997) with a measured position angle of PA_bar  =  150° and a length greater than 80″. Based on the bar properties, Davies et al. (2014) used a hydrodynamical simulation developed by Maciejewski (2004) to illustrate the effects of a strong stellar bar in the gas kinematics of NGC 3227. They showed that the simulated line-of-sight velocity field reproduced both the circular motions and the streaming motions at the outer part of the ring well, as observed in their H₂ map and also our CO(2–1) map. The streaming motions are expected in and near the dust lanes that enter the circumnuclear ring, which is exactly what we observe in our maps (see Figs. 1 and 7).

Here we evaluate if the perturbation from a large-scale bar alone is responsible for the CO(2–1) non-circular motions observed on different physical scales in NGC 3227. To do so, we used the same Maciejewski (2004) bar simulation (model termed 0S20r) as was used by Davies et al. (2014) and adapted it to the disk geometry fit with ^3DBAROLO. Then we smoothed it to the angular resolution of our ALMA CO(2–1) data (assuming a circular Gaussian with a FWHM of 0.2″) and scaled it to the observed velocities and derived v_sys. To compare with the observed CO(2–1) p − v diagrams, we extracted plots of the line-of-sight velocity against the distance from the center along the major and minor axes from the simulated velocity field of the gas. The comparison of Figs. 10 and 8 (top panels) shows that the bar simulation model is able to reproduce the line-of-sight velocity along the major axis as well as the streaming motions in regions along the minor axis at ∼1.5 − 3.5″, to the northeast and southwest of the AGN (marked with the magenta rectangles in Fig. 8). At these radial distances we are inside the large-scale bar corotation, therefore the streaming motions are associated with material inflowing near the leading edges of the bar. The bar model cannot reproduce the amplitude of the forbidden velocities along the major axis (marked with the blue rectangles in the left top panel of Fig. 8), however.

Fig. 10. Gas line-of-sight velocity versus distance from the center along the galaxy disk kinematic major (left panel) and minor (right panel) axes. The gas line-of-sight velocity field corresponds to the 0S20r simulation of a strong bar (Maciejewski 2004), and adapted to our ALMA resolution as well as to the NGC 3227 disk geometry, and observed CO(2–1) velocities. In each panel the cut passing through the center is plotted with the dot symbols. We also extracted two additional cuts on both sides of the center at ±0.1″ which are represented with the shaded area. These plots can be directly compared with the observed CO(2–1) p − v diagrams in the top panels of Fig. 8.

Along the minor axis, the bar model cannot reproduce the strong non-circular motions in the nuclear region (r <  1″). This finding is typical for the gas flow in barred potentials, which becomes close to circular inwards from the nuclear ring (e.g., Piner et al. 1995). In addition to the inflow along the bar, the bar triggers nuclear spiral shocks in model 0S20r, which can be revealed in the line-of-sight velocity residuals (e.g., Storchi-Bergmann et al. 2007; Davies et al. 2009), and which are manifested by the wiggles at radial distances r <  1″, in the nuclear region along the model minor axis p − v diagram (right panel of Fig. 10). However, the amplitude of these residuals is not sufficient to explain the observed velocities along the minor axis of NGC 3227.

5.2. Nuclear warp?

Schinnerer et al. (2000) observed similar non-circular motions in the nuclear region of NGC 3227 using lower resolution CO(2–1) observations. They first attempted to model them with a nuclear bar but were unable to reproduce their observed p − v diagrams in the inner ∼1″. They explored a second possibility using a circumnuclear disk that was warped at approximately radial distances of 0.7−1″, from the nucleus and was orthogonal to the host galaxy plane in the innermost region. They found that this model reproduced the CO(2–1) observations better and that the geometry would agree with the geometry required from the kinematics of the ionized gas of NGC 3227. Fischer et al. (2013) derived an angle between the axis of the ionization cone and the normal to the host galaxy disk of β = 76° (Fischer et al. 2013). This would therefore also imply a drastic gradient of the axis of the rotating molecular gas as it approaches the nuclear collimating structure.

In this subsection we briefly explore with ^3DBAROLO and our high angular resolution ALMA data the possibility of a nuclear warp. To do so, we took the initial model fit with ^3DBAROLO and refit the nuclear region (radial distances smaller than 1″) to allow for a slight geometric warp, in the sense of a small tilt between the inclination of the host galaxy and the nuclear region. This is because we do not know the scales on which this change in tilt occurs, that is, whether on the scales probed by our ALMA CO(2–1) data or much farther in. As before, we fixed the kinematic center and v_sys. For the nuclear warp we started with initial values of i_nuc = +30° and PA_{MAJ − nuc} = 320°. The change of almost 180° in the PA with respect to that of the host galaxy disk (i.e., the twist) is necessary to fit the counterrotation seen in the nuclear region and to conform with the orientation of the nuclear CO(2–1) clumps (see Fig. 5). Using these initial values, we allowed ^3DBAROLO to vary the nuclear and disk PA and inclination by ±10°. The second row of Fig. 8 shows that in the p − v diagram along the kinematic major axis, a nuclear warp starts populating some of the forbidden velocities.

We ran another ^3DBAROLO model allowing for a more dramatic warp to make it consistent with the inclination of the collimating structure (i_nuc = −30°) for the ionization cone. In other words, we assumed that the required warp takes place on the physical scales probed by our ALMA CO(2–1) data set. In this case we did not allow the nuclear PA and inclination to vary but kept them fixed (see Table 4 for a summary). The model p − v diagrams are rather similar (third row of Fig. 8) to those of the slight warp. This suggests that with this data set we cannot determine the physical scales where this warp would be taking place. We note that the model with the slight warp appears to produce a marginally better fit to the nuclear region along the major axis. This is because the fitted velocity dispersion in the nuclear region for the slight warp extends to values σ_gas ≃ 50 − 100 km s⁻¹ (but see Sect. 5.4), while the drastic warp model only fits values up to σ_gas ≃ 75 km s⁻¹. The effects of the increased velocity dispersion in the models with the warp are also visible in the p − v diagrams along the minor axis of the model when compared to the simple rotating disk. Even so, the nuclear warps cannot reproduce the strong non-circular motions to the northeast and southwest along the minor axis (which are also visible in the hot H₂ gas velocity map in Fig. 6) in the nuclear region.

5.3. Nuclear non-circular motions

The observed CO(3–2) p − v diagram along the minor axis of the nuclear disk (PA_{nuc − MIN} = 50°) in Fig. 11 clearly shows non-circular motions. The deviations from circular motions reach line-of-sight velocities of (v − v_sys)∼ + 200 km s⁻¹ (i.e., redshifted) to the northeast of the AGN out to a projected distance r = 0.5″, and (v − v_sys)∼ − 150 km s⁻¹ (i.e., blueshifted) to the southwest of the AGN at a similar distance. Therefore, these radial motions are observed to extend for approximately 1″ (∼70 pc) along the minor axis of the nuclear disk. These excess velocities with respect to a purely rotating nuclear disk are in agreement with the nuclear velocity residuals observed in hot molecular gas using the near-infrared H₂ line by Davies et al. (2014).

Fig. 11. CO(3–2) p − v diagrams along the major and minor kinematic axes of the nuclear disk. We only show fluxes above 3σ. The black contours are fluxes on a linear scale. The arrows show the approximate projected size of the 850 μm nuclear continuum emission.

To illustrate the effects of strong non-circular motions in the nuclear regions as observed along the kinematic minor axis, we ran another ^3DBAROLO model again using the initial rotating disk model and including a radial velocity component in the nuclear region. Based on the deviations discussed above, we fixed the radial velocity component in the nuclear region (r <  1″) to V_RAD = 180 km s⁻¹ (see Table 4). We allowed the rotational velocity and the velocity dispersion to vary and left the other parameters fixed. The fitted model is shown in the fourth row of Fig. 8. This model reproduces the p − v diagram along the minor axis and in the nuclear region fairly well, as is also more clearly visible in the right panel of Fig. 11. We note that this is not the best fit to the data because inflowing motions near the nuclear region are also observed to the southwest and northeast of the AGN (see above). However, with ^3DBAROLO we cannot simultaneously include radial velocity components with different signs at a particular radial distance. Thus we included the dominant radial component in the nuclear region, which produces a better fit to the nuclear part of the p − v diagram while keeping the fitted velocity dispersion relatively low (approximately σ_gas = 20 − 30 km s⁻¹).

As a final ^3DBAROLO model, we combined a nuclear warp with a radial component with the parameters listed in Table 4. This last model (Fig. 8, 5th row) visually produces a better fit to p − v diagrams along both the major and minor axes than just having the nuclear warp or the radial component alone. Additionally, it does not need high nuclear gas velocity dispersions at the transition region of the warp.

5.4. Nuclear molecular outflow

In the previous subsection we showed that non-circular motions along the minor axis of the galaxy are certainly detected within the central 1″, and more importantly, along the minor axis. When we assume that these motions are in the plane of the putative nuclear disk (at i_nuc = −30°), then these motions would be interpreted as inflow. However, the inferred radial velocities (up to 150 − 200 km s⁻¹/sin i_nuc) would be unrealistically high for a nuclear inflow. For instance, the hydrodynamical simulations we ran to fit the observations of the low-luminosity AGN NGC 1097 required radial flows of about 20 − 50 km s⁻¹ (Davies et al. 2009). Even in the presence of a nuclear bar, the inflow velocities are expected to be a fraction of the circular velocity (Haan et al. 2009), so that it is unlikely that these radial velocities would be associated with nuclear inflows.

Alternatively, we can assume that the radial motions occur in the plane of the host galaxy disk because they appear farther away from the nuclear disk, as seen from the projected vertical size of the 850 μm continuum map (see Fig. 11). In this case the non-circular motions would imply a nuclear molecular gas outflow that might be due to molecular gas in the disk of the galaxy being swept by the AGN wind. This has been observed in other nearby AGN, for instance, NGC 1068 (García-Burillo et al. 2014) and IC 5063 (Morganti et al. 2015). The observations would imply an outflow on physical scales of 70 pc with a maximum velocity of 150 − 200 km s⁻¹/sin i_disk ∼ 190 − 250 km s⁻¹. This velocity is consistent with the molecular outflows observed in a few Seyfert galaxies with similar AGN bolometric luminosities (see Fiore et al. 2017, for a compilation), although these are generally detected much further away from the AGN (several hundred parsecs). On nuclear scales, they are similar to those identified by Davies et al. (2014) with the near-infrared H₂ line in NGC 3227 and other Seyferts.

The CO(3–2) total intensity map of the nuclear region (Fig. 5) shows that the emission is not symmetric around the AGN location. Especially along the minor axis, the kinematics reveals a cavity like structure (Fig. 11) that is also visible in the CO intensity maps and even better from the hot molecular gas map in Fig. 6. It thus appears that molecular gas and possibly dust are more efficiently hollowed out, although not completely emptied, by the AGN wind and/or jet in the region immediately southwest of the AGN. This has also been observed in NGC 1068 (see Müller Sánchez et al. 2009; García-Burillo et al. 2014, and in prep.). This would explain why the spectral index between bands 6 and 7 (see Fig. 3) in this region appears to be mostly synchrotron-like emission, whereas the region northeast of the AGN shows evidence of both thermal and nonthermal emission. We also note that the nuclear outflow coincides with a region with strong evidence for shocks or outflow based on near-infrared [Fe II] observations (Schönell et al. 2019).

6. Nuclear molecular gas mass and column density

We measured the CO(2–1) and CO(3–2) line intensities over the nuclear region with a square aperture of 0.2″ × 0.2″ centered at the AGN position (see Table 3). We used the CO(2–1) measurement to estimate the molecular gas mass in this region. We took a CO(1–0)/CO(2–1) brightness temperature ratio of one, assuming a thermally excited and optically thick gas (Braine & Combes 1992), and used the relation of Sakamoto et al. (1999) with a Galactic CO-to-H₂ conversion factor X  =  2 × 10²⁰ cm⁻² (K km s⁻¹)⁻¹. We obtained a molecular gas mass in the inner 15 pc of approximately 5 × 10⁵ M_⊙.

Sani et al. (2012) detected HCN emission using observations with an approximate 1.2″ × 0.7″ resolution. The HCN transition peaks at the AGN position of NGC 3227, thus indicating that there is dense molecular gas. In this case, we would need to assume a brightness temperature ratio for the CO(1–0) and CO(2–1) transitions that is more appropriate for these gas conditions in the nuclear region. Viti et al. (2014) measured a (2–1) to (1–0) brightness temperature ratio of approximately 2.5 at the AGN location in NGC 1068. They also demonstrated that the observed suite of molecular gas transitions at the AGN location of NGC 1068 can be modeled with dense ( >  10⁵ cm⁻³) and hot (T >  150 K) gas. If the gas conditions in the nuclear region of NGC 3227 are similar to those in NGC 1068, then the derived mass of molecular gas (as well as the column density, see below) in the nuclear region should be taken as an upper limit. This is just due to the higher CO(2–1) to CO(1–0) ratio assumed for the calculation.

The nuclear molecular gas mass of NGC 3227 is similar to that of the low-luminosity AGN NGC 1097 (Izumi et al. 2017) but lower than the typical values measured by Combes et al. (2019) for low-luminosity AGN and other Seyferts (Alonso-Herrero et al. 2018; Izumi et al. 2018). However, it agrees remarkably well with the torus gas mass derived from SED fits of near- and mid-infrared emission (see García-Bernete et al. 2019). When we use a square aperture 0.5″, which would encompass the nuclear disk identified from the 850 μm continuum (Sects. 3.1 and 3.2), we derive a mass in molecular gas of 3 × 10⁶ M_⊙, which agrees well with the typical values for other nearby AGN. When we use the CO(3–2) measurement for this 0.5″, aperture and the expressions used by Combes et al. (2019), we derive 7 × 10⁶ M_⊙, which includes helium. For reference, the molecular gas mass estimates in the central 0.8″ are in the (2 − 20)×10⁷ M_⊙ range (Davies et al. 2006; Hicks et al. 2009), using the dynamical mass based on the near-infrared H₂ line and assuming a gas mass fraction of 10%.

From the CO(2–1) measurement in the nuclear region, we derive an average H₂ column density toward the AGN location (inner ∼15 pc) of approximately N(H₂) = 2 − 3 × 10²³ cm⁻² or a factor of two lower if there is a significant amount of dense molecular gas in the nuclear region. This value would appear to be high considering that NGC 3227 is optically classified as a Seyfert 1.5. However, X-ray, UV, and optical studies have shown much evidence for different absorbers toward the AGN line of sight, including a dusty warm absorber, lukewarm absorbers, and cold absorbers (see, e.g., Komossa & Fink 1997; Crenshaw et al. 2001, and references therein). Moreover, this complex X-ray absorption is also variable (see, e.g., Turner et al. 2018). Recently, Beuchert et al. (2015) derived column densities for the ionized absorbers of N_H ∼ 5 − 16 × 10²² cm⁻². They also discussed that the absorbers could be part of an overall clumpy medium located in the outermost dust-free broad-line region or the inner part of the dusty torus.

7. Discussion and conclusions

We have presented observations of ALMA bands 6 and 7 of the molecular gas using the CO(2–1) and CO(3–2) transitions and their associated continua of the nuclear and circumnuclear regions of the nearby Seyfert 1.5 galaxy NGC 3227. The ALMA observations have angular resolutions between 0.085 and 0.21″, which for this galaxy correspond to physical resolutions in the range 7−15 pc. With these observations we have studied the molecular gas morphology and kinematics in detail, as well as the (sub-) millimeter continuum emission. Figure 12 summarizes our results.

Fig. 12. Summary of results for NGC 3227. The large RGB image to the left was produced using the ALMA CO(2–1) intensity in blue, the HST/ACS F814W image in green, and the HST/ACS F658N image in red. The FoV is approximately 10″ × 10″ as in Fig. 1. The white straight line marks the PA of the kinematic major axis of the large-scale disk, while the dashed magenta lines delineate the projected size of the ionization cone as seen in the right panel of Fig. 1. The square approximately marks the central region with a size of 1″ × 1″ that is displayed in Figs. 2 and 5. The small figure to the right is a cartoon (not to scale) representing the nuclear disk mostly traced by the 850 μm continuum and the 1.3 mm extended emission perpendicular to the nuclear disk. In both panels the arrows and their colors represent the non-circular motions observed with the CO(2–1) and CO(3–2) transitions.

From a morphological point of view, the ALMA continuum emission at 1.3 mm and 850 μm in the innermost ∼70 pc shows an unresolved component that is likely to be associated with the AGN itself, as well as extended components. One extended component is along the projected direction of the ionization cone axis (PA_cone ∼ 30 − 40°) and the other along the position of the nuclear star-forming disk (PA ∼ −30°) that is identified with near-infrared observations (Davies et al. 2006). We show them in the cartoon in the right panel of Fig. 12. In the direction of the ionization cone axis, the spectral index between the two ALMA bands shows mostly positive values toward the northeast of the AGN and negative values to the southwest. For the latter we can conclude that the ALMA continuum emission is dominated by synchrotron emission that is likely associated with a putative radio jet. Along the nuclear disk extension, we derive spectral indices from slightly negative to positive values, which are consistent with a combination of thermal emission from cold dust and free emission from H II regions and nonthermal emission due to the AGN.

Our analysis of the nuclear radio-to-infrared SED of NGC 3227 shows that the ALMA continuum at 1.3 mm and 850 μm could be produced by a combination of different physical processes, as have been found for other nearby Seyferts (Krips et al. 2011; García-Burillo et al. 2014; Pasetto et al. 2019). In particular, the extrapolation to the far-infrared of the Nenkova et al. (2008a); Nenkova et al. (2008b) clumpy torus model fit to the infrared emission would only account for a small fraction of the observed 850 μm flux measured over 0.2″, (15 pc). This could be due to a significant contribution from star formation to the submm emission in the nuclear region.

The exquisite ALMA angular resolution reveals the distribution of the molecular gas in the nuclear and circumnuclear regions, with a ring-like morphology (diameter of ∼350 pc) that was previously identified by Schinnerer et al. (2000). The overall (FoV ∼ 10″ × 10″) kinematics of the CO(2–1) line shows circular motions consistent with those derived from the stellar kinematics (Davies et al. 2006; Barbosa et al. 2009). We also detect strong streaming motions (inflow) near the outer edges of the ring to the southeast and northwest of the AGN (at projected radial distances ∼3 − 4″) that can be explained with the presence of a large-scale bar (Mulchaey et al. 1997; Davies et al. 2014). These streaming motions appear to continue inward to regions as close as 0.5 − 0.7″ from the AGN.

In the innermost regions, the CO(2–1) and CO(3–2) maps show a number of clumps with a morphology similar to that of the hot molecular gas traced by the near-infrared H₂ line (Davies et al. 2006). None of the CO or H₂ peaks coincides with the AGN location (see Figs. 5 and 6), which is assumed to correspond with the position of the unresolved 1.3 mm emission. However, we also showed that the extended 1.3 mm and 850 μm emission in the nuclear region traces the distribution of the molecular gas well. The nuclear region (inner ∼1″  =  73 pc) has extremely complicated kinematics with strong deviations from pure rotating motions along the kinematic major and minor axes. We explored the possibility of a nuclear warp that would be required to explain the angle between the axis of the cone and the normal to the galaxy plane of β  =  76° (Fischer et al. 2013). Although we cannot conclusively demonstrate the presence of a nuclear warp, a significant change in tilt between the host galaxy and the region collimating the ionization cone needs to occur in the inner 1″.

Both the CO(2–1) and the CO(3–2) transitions show strong non-circular motions along the minor kinematic axis in the nuclear region (Figs. 8 and 11). They reach maximum velocities of (v − v_sys)∼190 − 250 km s⁻¹ (inclination-corrected values) that are blueshifted to the southwest of the AGN and redshifted to the northeast (see also the cartoon in the right panel of Fig. 12). These motions were well reproduced with ^3DBAROLO by adding a nuclear radial velocity component to the host galaxy rotating disk. Assuming that they are in the plane of the host galaxy disk, we can interpret them as due to a nuclear molecular outflow over scales of approximately 70 pc. The molecular gas outflow velocities are lower than those derived from the ionized gas using the Hα and [S III]λ9069 emission lines, which in this region reach maximum values of ±500 km s⁻¹ and ∼ − 900 km s⁻¹, respectively (Fischer et al. 2013; Barbosa et al. 2009). This is understood if the molecular outflow in NGC 3227 is due to gas in the host galaxy entrained by the faster AGN wind.

Morphologically and kinematically, the AGN wind appears to have excavated the molecular gas and possibly the dust in regions close to the AGN. This created a cavity-like structure, which is more apparent to the southwest of the AGN (see Fig. 5 and especially Fig. 6) and is similar to the cavity observed in NGC 1068 (see Müller Sánchez et al. 2009; García-Burillo et al. 2014, and in prep.). We can estimate the nuclear outflow rate on both sides of the AGN along the minor axis taking expression given in Eq. (4) of García-Burillo et al. (2014): dM/dt = 3 × v_out × M_mol/(R_out × tan(i_disk)) for a multi-conical outflow that is uniformly filled by the outflowing clouds. We measured the CO(2–1) fluxes in square apertures of 0.2″, on a side just outside of the central AGN position and estimate molecular gas masses (following the same assumptions as in Sect. 6) of 7 × 10⁵ M_⊙ and 8 × 10⁴ M_⊙ to the northeast and southwest of the AGN, respectively. For the outflow velocity we took a typical value of v_out = 100 km s⁻¹ and a size for the outflow of R_out = 35 pc. This results in nuclear molecular outflow rates of dM/dt = 5 M_⊙ yr⁻¹ and dM/dt = 0.6 M_⊙ yr⁻¹ to the northeast and southwest of the AGN. These values are about an order of magnitude lower than the molecular outflow rate at the circumnuclear disk of NGC 1068 (García-Burillo et al. 2014), but are comparable to the rate measured for NGC 1433 (Combes et al. 2013).

At the AGN position and with our angular resolution, this region is not entirely devoid of molecular gas. We estimated a mass in molecular gas of 5 × 10⁵ M_⊙ (for a square aperture 0.2″ ∼ 15 pc on a side) and an equivalent average column density of N(H₂) = 2 − 3 × 10²³ cm⁻² at the AGN position. However, if the molecular gas at the AGN location is dense and relatively hot, as in NGC 1068, the derived molecular gas mass and column densities would be approximately a factor of two lower. The lower values are due to the higher CO(2–1) to CO(1–0) ratio that is observed in the nuclear region of NGC 1068 compared to the typical ratios observed in the disks of galaxies. Nevertheless, although NGC 3227 is classified as a Sy1.5, there is evidence of a gas-rich nuclear or circumnuclear environment. This agrees with the weak relationship between the circumnuclear (1 − 3 kpc) and nuclear (tens of parsecs) molecular gas mass and the X-ray column density in nearby AGN (see Rosario et al. 2018; Combes et al. 2019, respectively) and implies that the molecular gas has a low filling factor. As a quick estimate, the gas mass and size of the region give an average filling factor of molecular gas, assuming a density of ∼10⁶ cm⁻³, of 6 × 10⁻³, or a ∼80% chance that a given line of sight will not intersect a molecular clump. In this sense, it is not surprising that the AGN is partially unobscured.

In summary, the nuclear CO(2–1), CO(3–2) and dust morphology of NGC 3227 do not resemble a classical compact torus. Rather, the dust and molecular gas emission is in a nuclear disk that is extended over several tens of parsecs. The orientation of the nuclear disk is consistent with being part of the collimating structure of the ionization cone. We find evidence that the dust might be heated not only by the AGN, but also by the nuclear ongoing or recent star formation activity (Davies et al. 2006; Esquej et al. 2014). In the CO(2–1) and CO(3–2) transitions we do not observe a bright and isolated nuclear structure or torus, but several clumps that are well connected with the circumnuclear ring in the host galaxy disk. In this situation, it is difficult to say where the torus starts, although such molecular gas morphologies are also common in other AGN of similar or lower luminosity (Izumi et al. 2017, 2018; Alonso-Herrero et al. 2018; Combes et al. 2019).

Acknowledgments

We are grateful to the anonymous referee for constructive comments that helped us improve the manuscript. We thank Luis Colina, Alvaro Labiano, Leonard Burtscher and Daniel Asmus for useful discussions. A.A.-H. and I.G.-B. acknowledge support from the Spanish Ministry of Science, Innovation and Universities through grant AYA2015-64346-C2-1-P, which was party funded by the FEDER program. A.A.-H. also acknowledges CSIC grant PIE201650E36. A.A.-H., S.G.-B. and A.U. acknowledge support from the Spanish Ministry of Science, Innovation and Universities through grant PGC2018-094671-B-I00, which was party funded by the FEDER program. S.G.-B. and C.R.A. acknowledge the Spanish Plan Nacional de Astronomía y Astrofísica under grant AYA2016-76682-C3-2-P. M.P.-S. and D.R. acknowledge support from STFC through grant ST/N000919/1. T.D.-S. acknowledges support from ALMA-CONICYT project 31130005 and FONDECYT regular project 1151239. S.F.H. is supported by European Research Council Starting Grant ERC-StG-677117 DUST-IN-THE-WIND. C.R.A. also acknowledges the Ramón y Cajal Program of the Spanish Ministry of Science and Technology through project RYC-2014-15779. S.R. and M.V. gratefully acknowledge support from the Independent Research Fund Denmark via grant numbers DFF 4002-00275 and 8021-00130. L.K.H. acknowledges funding from the INAF PRIN-SKA 2017 program 1.05.01.88.04. C.R. acknowledges support from the CONICYT+PAI Convocatoria Nacional subvencion a instalacion en la academia convocatoria año 2017 PAI77170080. D.R. acknowledges the support of the Science and Technology Facilities Council (STFC) through grant ST/P000541/1. AU acknowledges support from the Spanish MINECO grants ESP2015-68964 and AYA2016-79006. This paper makes use of the following ALMA data: ADS/JAO.ALMA.2016.1.00254.S and JAO.ALMA.2016.1.01236.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST 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 made with the NASA/ESA Hubble Space Telescope, obtained from the data archive at the Space Telescope Science Institute. STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555.

References

All Tables

All Figures

Fig. 1. Left panel: ALMA CO(2–1) integrated molecular line map (zeroth-order moment map) created using a 3σ threshold. The image displays the central 10″ × 10″ (730 pc × 730 pc) and was created from the b = 0.5 weighted data cube. We show the image and contours with a linear scale. The CO(2–1) contours are 0.0925, 0.515, 0.9375, 1.36, 1.7825, 2.205, 2.6275, and 3.05 Jy km s−1 beam−1. The star indicates the AGN location, based on the peak of the unresolved 1.3 mm continuum emission (see text). Right panel: HST/ACS pseudo-Hα+[N II] line map (highlighted in light colors) plus extinction map (highlighted in dark colors) constructed as an F658N/F814W map (in arbitrary units) and showing the same FoV as the ALMA CO(2–1) image. The broadband F814W image is saturated at the nuclear position.

In the text

Fig. 2. ALMA continuum images in band 6 at ∼1.3 mm from the b = 0.5 weighted data cube (left) and b = −0.5 weighted data cube (middle) and in band 7 at ∼850 μm from the natural weight data cube (right). The FoV is 1″ × 1″, which corresponds to approximately the central 73 pc × 73 pc region. In all images the first contour corresponds to the 2.5σ level and the color bars indicate the flux in units of 10−4 Jy beam−1. The hatched ellipses show the corresponding beam size and orientation in each image (see Table 1), and the star indicates the approximate location of the peak of the unresolved 1.3 mm continuum emission (see text).

In the text

Fig. 3. Top panel, color and white contours: map of the 850 μm continuum emission (in arbitrary units) degraded to the beam size and orientation of the 1.3 mm continuum map with b = 0.5 weight (see Table 1). This was necessary to construct the spectral index map. The map is shown for the nuclear 1″ × 1″ region. We mark the approximate orientations of the ionization cone axis and the nuclear star forming disk. Bottom panel: map of the spectral index α as derived from the ALMA continuum observations in bands 6 and 7 using only pixels detected at > 2.5σ. The FoV is the same as in the top panel. The map has the same angular resolution as the band 6 map with the b = 0.5 weight (see Table 1). The color bar indicates the values of α (see text). The contours are the band 6 continuum emission at 1.3 mm as in the left panel of Fig. 2, and the star marks the AGN position.

In the text

Fig. 4. Nuclear continuum SED from centimeter wavelengths to the infrared of NGC 3227. The small-aperture (∼0.1 − 0.2″) MERLIN data points are for the southern source of Mundell et al. (1995b). The large-aperture (> 0.2″ to ≤1″) VLA and PdBI data points are from Nagar et al. (2005) and Sani et al. (2012), respectively. All the ALMA continuum data points are from this work. The very small aperture (a few milli-arcseconds) points are the EVN radio fluxes from Bontempi et al. (2012). The fit (black dashed line, with the ±1σ uncertainty plotted with the shaded gray region) is made to the near- and mid-infrared emission using the Nenkova et al. (2008a,b) clumpy torus models (see García-Bernete et al. 2019, for more details). We also plot the observed mid-infrared data point for reference although the fitted data also included near-infrared photometry and mid-infrared spectroscopy. The magenta dashed line is the fit to the MERLIN radio points extrapolated to the ALMA wavelengths.

In the text

Fig. 5. Zoom-in on the CO(2–1) and CO(3–2) emission in the nuclear 1″ × 1″ region in the top and bottom panels, respectively. These moment maps were created with the GILDAS MOMENT task using a 5σ detection threshold and show the integrated emission (zeroth-order moment map) in the left panels and the mean velocity field (first-moment map) in the right panels. The hatched ellipses show the corresponding beam size and orientation of each image (see Table 2). The contours in all panels are the ALMA continuum emission at 1.3 mm as in the left panel of Fig. 2, and the star marks the AGN position.

In the text

Fig. 6. VLT/SINFONI H2 1–0 S(1) intensity (left) and velocity (right) maps from Davies et al. (2006) showing a FoV similar to that of Fig. 5. The maps were derived by fitting a Gaussian to the emission line. The black dot marks the approximate location of the peak of the near-infrared continuum. The velocities are relative to this peak. To facilitate the comparison with the CO(2–1) and CO(3-2) emission, we show the same velocity range as in the right panels of Fig. 5.

In the text

Fig. 7. 3DBAROLO maps of the observed CO(2–1) mean velocity field (first-order moment, left panel), the fitted model of a rotating disk with fixed PA, inclination (see Table 4 and text), systemic velocity, and kinematic center (middle panel), and the residual mean-velocity field (right panel). The FoV is 10″ × 10″, as in Fig. 1. The color bars indicate line-of-sight velocities in units of km s−1. The contours for the velocities in the left and middle panels are in steps of 50 km s−1. Right panel: the contours are the CO(2–1) intensity (zeroth-order moment map) on a linear scale.

In the text

Fig. 8. Grey-scale and black contours: observed CO(2–1) p − v diagrams along the kinematic major axis (left panels) and minor axis (right panels) of the host galaxy disk. We only show fluxes above 5σ. In all the panels the red contours represent the 3DBAROLO models as follows (see also Table 4): the first row is a rotating disk, the second row is a rotating disk with a slight nuclear warp, the third row is a rotating disk + a warp model consistent with the orientation of the ionization cone, the fourth row is a rotating disk with a nuclear radial component, and the fifth row is a rotating disk with a nuclear warp + radial component. The model contour levels are identical to the contour level of the observations. The vertical dotted black line marks the AGN position, and the horizontal dotted black line shows the derived systemic velocity. Top panels: we also highlight nuclear regions of forbidden velocities for a purely rotating disk. Along the minor axis the blue boxes indicate the nuclear outflow and the magenta ones regions with streaming motions due to the large-scale stellar bar.

In the text

	Fig. 8. continued.
In the text

	Fig. 9. Observed CO(3–2) p − v diagrams along the kinematic major axis (left) and minor axis (right) of the host galaxy disk. We only show fluxes above 3σ. The black contours are fluxes in a linear scale. The black dotted vertical line marks the AGN position, whereas the horizontal red lines mark the approximate velocity virial range of the simple rotating disk model (see top panel of Fig. 8).
In the text

Fig. 10. Gas line-of-sight velocity versus distance from the center along the galaxy disk kinematic major (left panel) and minor (right panel) axes. The gas line-of-sight velocity field corresponds to the 0S20r simulation of a strong bar (Maciejewski 2004), and adapted to our ALMA resolution as well as to the NGC 3227 disk geometry, and observed CO(2–1) velocities. In each panel the cut passing through the center is plotted with the dot symbols. We also extracted two additional cuts on both sides of the center at ±0.1″ which are represented with the shaded area. These plots can be directly compared with the observed CO(2–1) p − v diagrams in the top panels of Fig. 8.

In the text

	Fig. 11. CO(3–2) p − v diagrams along the major and minor kinematic axes of the nuclear disk. We only show fluxes above 3σ. The black contours are fluxes on a linear scale. The arrows show the approximate projected size of the 850 μm nuclear continuum emission.
In the text

Fig. 12. Summary of results for NGC 3227. The large RGB image to the left was produced using the ALMA CO(2–1) intensity in blue, the HST/ACS F814W image in green, and the HST/ACS F658N image in red. The FoV is approximately 10″ × 10″ as in Fig. 1. The white straight line marks the PA of the kinematic major axis of the large-scale disk, while the dashed magenta lines delineate the projected size of the ionization cone as seen in the right panel of Fig. 1. The square approximately marks the central region with a size of 1″ × 1″ that is displayed in Figs. 2 and 5. The small figure to the right is a cartoon (not to scale) representing the nuclear disk mostly traced by the 850 μm continuum and the 1.3 mm extended emission perpendicular to the nuclear disk. In both panels the arrows and their colors represent the non-circular motions observed with the CO(2–1) and CO(3–2) transitions.

In the text