Free Access
Issue
A&A
Volume 529, May 2011
Article Number A45
Number of page(s) 17
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/201014488
Published online 29 March 2011

© ESO, 2011

1. Introduction

Gas accretion onto supermassive black holes (SMBHs) is believed to be the cause of nuclear activity in galaxies. SMBHs are accepted as a common component in most galaxies with a significant massive bulge (e.g. Ferrarese & Ford 2005, and references therein). However, only about half (43%) of the local galaxies host active galactic nuclei (AGN) in the Seyfert, LINER or transition-object categories (Ho et al. 1997). This discrepancy between the presence of SMBHs and nuclear activity must therefore be sought in the possibility and efficiency of gas transportation to the central regions. Gas transport toward the centers of galaxies can only happen when the gas is able to lose its angular momentum. Two categories of dynamical mechanisms can cause inflow. Gravitational mechanisms such as galaxy-galaxy interactions on large scales or non-axisymmetries (i.e. spiral density waves, large-scale stellar bars or nuclear ovals) within the galaxy potential exert gravitational torques. Viscosity torques and shocks caused by turbulence in the interstellar medium (ISM) are a hydrodynamical mechanism for losing angular momentum. When present, gravitational torques are more efficient (e.g. Combes 2004).

Observations of the inner regions of galaxies are necessary to understand how the interplay of mechanisms results in gas transportation all the way down to the nucleus. The dynamic timescales in the central region of galaxies are short. Therefore, the gas distribution in these regions needs to be mapped with high angular resolution. The NUclei of GAlaxies (NUGA) project (García-Burillo et al. 2003a) has been obtaining high-resolution (0.5″ − 1″) detailed mapping of the molecular gas kinematics in 12 nearby (D = 4 − 40 Mpc) low luminosity active galactic nuclei (LLAGN) with the IRAM PdBI (Plateau de Bure Interferometer) and 30 m telescope. This sample spans the whole sequence of nuclear activity types. In the central kiloparsec most of the gas is in the molecular phase, making CO lines optimal tracers of the gas dynamics. The spatial resolution ( < 100 pc) of this survey allows one to observe the gaseous distribution over an impressive spatial range. This has already led to the identification of a wide range of morphologies in the nuclear regions of these galaxies, including lopsided disks (NGC 4826: García-Burillo et al. 2003b; NGC 3718: Krips et al. 2005; NGC 5953: Casasola et al. 2010), bars and spirals (NGC 4569: Boone et al. 2007; NGC 2782: Hunt et al. 2008; NGC 6574: Lindt-Krieg et al. 2008; NGC 4579 García-Burillo et al. 2009) and rings (NGC 7217: Combes et al. 2004; NGC 3147: Casasola et al. 2008; NGC 1961: Combes et al. 2009).

Large-scale stellar bars are believed to be efficient in driving gas towards the inner Lindblad resonance (ILR) (i.e. Regan & Teuben 2004; Sakamoto et al. 1999; Sheth et al. 2005). There they induce the formation of spiral or ring structures (Prendergast 1983; Athanassoula 1992a,b; Englmaier & Shlosman 2000; Maciejewski 2004; Combes 2004). Martini & Pogge (1999) and Regan & Mulchaey (1999) have proposed that spiral shocks induced by the large-scale stellar bar can generate further gas inflow across the ILR. Numerical simulations have shown that spiral structure can extend across the ILR into the nuclear region, if the sound speed is high enough (Englmaier & Shlosman 2000) or if the velocity dispersion in the ISM is large enough (Maciejewski 2004).

Alternatively, other mechanisms on smaller scales could take over gas transportation towards the nucleus, i.e. the bars-within-bars scenario (Shlosman et al. 1989). Viscosity torques can only become the dominant mechanism for inflow in the innermost regions ( < 200 pc: Combes 2004).

In this paper we study the nearby spiral galaxy NGC 6951. Its distance of 24.1 Mpc allows for high spatial resolution (1″ = 117 pc). There are several arguments that gas inflow is currently occurring in this galaxy. A large-scale stellar bar has been detected in the near-infrared (Buta et al. 2003; Block et al. 2004) with a bar radius of 26″ ( ≈ 3.0 kpc) and a position angle of 84° (Mulchaey et al. 1997). This galaxy also shows a pronounced starburst (SB) ring at 5″ ( = 580 pc) radius in Hα (Pérez et al. 2000; Rozas et al. 2002) and radio (Saikia et al. 2002) emission. Haan et al. (2009) have shown that the inner region of NGC 6951 is HI depleted, implying the ISM must be dominated by its molecular phase. But 12CO(1−0) and 12CO(2−1) emission associated with the SB ring has been found (Kohno et al. 1999; García-Burillo et al. 2005), as well as HCN(1−0) emission (Krips et al. 2007). NGC 6951 has also been classified as “a high excitation LINER and a possibly high nitrogen abundant Seyfert 2” galaxy (Pérez et al. 2000). In the central 1″ of NGC 6951 Krips et al. (2007) and García-Burillo et al. (2005) have detected HCN(1−0) and 12CO(2−1) emission, respectively, indicating that further inflow of gas beyond the SB ring must be or has been occuring recently.

We present 12CO(1−0) and 12CO(2−1) observations made with the IRAM PdBI and 30 m telescope. Previous papers by García-Burillo et al. (2005) and Haan et al. (2009) presented the PdBI-only data. Here the 12CO(1−0) and 12CO(2−1) data cubes have been combined with the 30 m observations. The addition of the 30 m observations provides the full CO emission present, sampled on all scales in the inner 3 kpc of NGC 6951.

This paper has three aims. The first is to investigate how the addition of the 30 m data changes the maps and the gravitational torque results derived by García-Burillo et al. (2005) and Haan et al. (2009) and to perform a more detailed investigation of the kinematics of the molecular CO gas than has been done before. The second is to quantify the influence of the large-scale stellar bar within the inner 3 kpc using a parametric kinematic model. The third is to study the impact of the nuclear stellar oval on the gas flow inside the circumnuclear gas ring.

In Sect. 2 we present the observations from the IRAM PdBI and 30 m telescope and their reduction. Section 3 contains a discussion of the changes due to the addition of the 30 m data and a presentation of the spatial and kinematic properties of the CO emitting gas as seen in the PdBI+30 m data cubes. In Sect. 4 we detail the large-scale stellar bar model, with the results that follow from the model being presented in Sect. 5. Finally, in Sect. 6 we compare our result with previous gravitational torque studies of NGC 6951 and discuss observational evidence for inflow to the nucleus. We conclude in Sect. 7.

2. CO observations

2.1. IRAM PdBI observations

The IRAM PdBI observations in ABCD configuration were carried out between June 2001 and March 2003 using the 6-antenna array in dual-frequency mode. Only the D-configuration observations were executed with 5 antennas. The correlator was centered at 114.726 GHz and 229.448 GHz at 3 mm and 1 mm, respectively, corresponding to a heliocentric velocity of 1425 km s-1. The bandwidth covered changed between the CD and AB configurations, although all four configurations covered at least the central  ± 200 km s-1. The flux calibration used CRL618, 0932+392, 3C 273 and/or 3C 345. Bandpass correction was derived from observations of a strong quasar at the beginning of the track. 1928+738 and 2010+723 served as the phase calibrators, allowing for correction of atmospheric effects. Phase corrections derived for the 3 mm receiver were applied to the 1 mm band resulting in a better phase correction at 1 mm. The data were reduced using standard routines in the GILDAS software package1.

For both data sets CLEANed data cubes were produced with natural and robust2 weighting. Cleaning was done down to the 2σ noise level within a fixed polygonal area that was defined based on the zeroth moment map for all channels with line emission. The rms noise for these data cubes is listed in Table 1.

2.2. IRAM 30 m observation

30 m observations of the central 132″ by 66″ were obtained on December 24 and 25, 1997. The 3 mm and 1 mm receivers were tuned to 114.730 GHz and 229.460 GHz. A bandwidth of about 1200 km s-1 [600 km s-1] was covered by 512 channels with a width of 2.6 km s-1 [1.6 km s-1] at 3 mm [1 mm]. The spacing between individual grid points was 11″, i.e. half the size of the 3 mm beam. The integration time per scan was usually 4 min, and both polarizations were simultaneously observed for each frequency. Typical system temperatures during the observations were 350 K and 630 K for the 3 mm and 1 mm receivers. The data reduction was done using the GILDAS/CLASS software package. Each scan was inspected and those few with extremely high system temperatures or other instrumental effects were discarded. The baseline was corrected in the individual spectra by fitting a first order polynomial through channels outside the expected line emission. After this correction all spectra for an individual position were averaged together using a noise weight.

2.3. Short spacing correction

The 30 m observations were used to compute the short spacing correction (SSC) and recover the large-scale low-level flux. The 30 m observations were reprojected to the field center and frequency of the PdBI observations. The bandwidth coverage of the 30 m observations was resampled to match the velocity axis of the interferometric observations. A combined data cube was produced using the task “UV-short” in GILDAS. The single dish weight scaling factor was 6.55 × 10-3 for 12CO(1−0) and 7.74 × 10-4 for 12CO(2−1).

Two sets of final CLEANed data cubes were produced using natural and robust weighting. The resolution of the natural [robust] weighted 12CO(1−0) and 12CO(2−1) data is given in Table 2. The data cubes have 512 × 512 pixels, with a pixel scale of 0.25″/pixel [0.10″/pixel] and velocity bins of 10 [5] km s-1 for the 12CO(1−0) [12CO(2−1)] data. The rms noise per channel in the SSC 12CO(1−0) observations is 2.5 mJy/beam and 2.2 mJy/beam with natural and robust weighting, respectively. For the SSC 12CO(2−1) observations these values are both 7.8 mJy/beam. CLEANing was done down to the 2σ noise level with the assistance of individual polygons defined for each channel with line emission present.

The noise values of the SSC data cubes are, with the exception of the 12CO(2−1) natural weighted data cube, below the noise levels of the PdBI only data. The beam sizes increase slightly (~8%) due to the added short spacings.

Table 1

Overview of PdBI data.

Table 2

Overview of PdBI+30 m data.

2.4. HST/NICMOS observations

We retrieved from the HST archive the NICMOS F110W and F160W images of NGC 6951. The reduction was carried out using the “best” calibration files, and the van der Marel algorithm (e.g. Böker et al. 1999) to reduce the “pedestal” effect. Sky values were assumed to be zero, which is generally a good assumption for these kinds of NICMOS images (Hunt & Malkan 2004). The images were calibrated, converted to magnitude scale, and subtracted to obtain a  [F110W]  −  [F160W]  ≈ J − H color map. This color image and the starburst ring it reveals will be discussed in Sect. 3.1.

3. Properties of the CO-emitting gas

3.1. Morphology and H2 masses

thumbnail Fig. 1

Top: integrated 12CO(1−0) emission in natural a) and robust b) weighting for the SSC data. The zeroth moment maps have been primary beam corrected. The 12CO(1−0) emission has been integrated from –220 km s-1 to 250 km s-1. Contours run from 5σ in steps of 10σ. The rms value is: a) 1σ = 0.025 Jy beam-1 km s-1 and b) 1σ = 0.022 Jy beam-1 km s-1. Bottom: integrated 12CO(2−1) emission in natural c) and robust d) weighting for the SSC data. The 12CO(2−1) emission has been integrated from –200 km s-1 to 200 km s-1. Contours run from c) 5σ in steps of 20σ, d) 3σ in steps of 10σ. The rms value is: 1σ = 0.039 Jy beam-1 km s-1 for both maps. The red cross indicates the position of the dynamical center in all panels (Table 3). The beam sizes are shown in the lower left corners and correspond to the values listed in Table 2.

Open with DEXTER

The intensity maps of the CO(1−0) and CO(2−1) observations presented in Fig. 1 have been constructed from the CLEANed data cubes using the software GIPSY3 (van der Hulst et al. 1992). These zeroth moment maps are computed as the pixel-wise sum of emission above a fixed threshold. Here we chose the 3σ noise level. Spurious signals are filtered out by imposing the constraint that the emission be above this threshold in at least two consecutive velocity channels. All signals that satisfy these requirements are added to the intensity map. To the CO(1−0) emission zeroth moment maps we have applied a primary beam correction.

The intensity distribution of the observed CO(1−0) and CO(2−1) emitting gas has all the components we would expect from a gas distribution driven by a large-scale stellar bar. The morphology of the CO(1−0) emission map (Fig. 1, top, left) from larger radii inward is comprised of the following components. At the edge of our field we see a spatially unresolved resonance ring at a radius of  ~ 30″ ( = 3.5 kpc). This is most likely at the location of the ultra harmonic 4:1 resonance inside corotation (e.g. Sakamoto et al. 2000). This ring is connected to straight gas lanes, with an approximately horizontal (east-west) orientation. Their position corresponds well to the orientation of the large-scale stellar bar, which has a position angle of 84° (Mulchaey et al. 1997). The gas lanes also coincide with the dust lanes at the leading edges of the stellar bar (Pérez et al. 2000, Fig. 4a of their paper). In the robust weighted map of the CO(1−0) emission (Fig. 1, top, right), which is more sensitive to emission on smaller scales, the straight gas lanes are no longer detectable. This indicates that the straight gas lanes consist of more diffuse, low intensity gas.

Going further inward, we find the straight gas lanes connect to a “twin peaks” morphology and a spiral pattern. This was also observed by Kohno et al. (1999). It is likely that the “twin peaks” are due to the crowding of gas streamlines in a barred potential (Kenney et al. 1992). The orbit crowding leads to a buildup of gas at those locations where orbits change from x1 to x2, i.e. with orientation along the bar major axis changing to orientation along the bar minor axis for elliptical orbits. The peaks inside the spiral arms have a distance of  ~ 6″ ( = 700 pc) from the nucleus. Slightly within that radius (~5″, 580 pc) a circumnuclear SB ring has been detected in Hα (Pérez et al. 2000; Rozas et al. 2002) and radio emission (Saikia et al. 2002).

The natural and robust weighted intensity CO(2−1) maps (Fig. 1, bottom) show a distribution very similar to the one seen in the CO(1−0) maps. Some of the differences, however, stem from the smaller field-of-view (FOV) of the CO(2−1) data. The outer ring at the edge of the CO(1−0) FOV is not visible, nor are the straight gas lanes. Due to the higher resolution, the two prominent peaks seen in the CO(1−0) maps now break up into multiple maxima. In the natural weighted map the northern peak is now joined by a second one slightly (φ ~ 30°) offset to the west. The southern spiral structure shows an elongated ridge with a far less distinct center. In the natural weighted CO(2−1) map the peaks and ridge are still unresolved structures.

thumbnail Fig. 2

12CO(2−1) robust weighted intensity map (contours from 3σ in steps of 10σ steps, 1σ = 0.039 Jy beam-1 km s-1); overlayed on a J − H color map (greyscale) from HST. White corresponds to higher extinction.

Open with DEXTER

thumbnail Fig. 3

Position-velocity diagram along a PA of 113° through the dynamical center. Contours are at 2σ, in 1σ steps, with the –2σ and –1σ contours given in grey. A gas bridge is visible between the northern part of the ring (lower left here) and the nucleus (denoted with a cross).

Open with DEXTER

thumbnail Fig. 4

Emission line ratio CO(2−1)/CO(1−0) in temperature units. This ratio was computed using the natural weighted CO(2−1) and the robust weighted CO(1−0) data cubes, smoothed to the resolution of the former. The high ratio south-east of the nucleus inside the ring is not significant.

Open with DEXTER

We start to resolve two tightly wound gas spiral arms, as well as peaks of dense gas within the ring only in our highest resolution robust weighted CO(2−1) maps. In the higher density tracer HCN(1−0) a similar effect is seen (Krips et al. 2007). The peaks found in that high resolution image, although barely resolved, sit almost 90° down the spiral arms. This seems to indicate some density gradient in the gas ring, starting at the northern and southern peaks and going downstream. Figure 2 shows the CO(2−1) integrated intensity map as contours superimposed on the NICMOS J − H color map (Sect. 2.4). The figure shows that the CO spiral arms coincides nicely with the SB ring seen in the red J − H color. The extinction, AV, can be derived from this color image assuming a J − H color of  ~ 0.7 for a normal (unreddened) composite stellar population, and (AJ − AH)/AV = 0.092 (as in Cardelli et al. 1989, with RV = 3.1). The AV’s so derived for the SB ring are quite red, ranging from  ≲ 1 to more than 6 in the dustiest, most compact regions.

Table 3

Properties of NGC 6951.

Finally, at the dynamical center (Table 3) we find barely resolved molecular gas emission (M ≈ 106 M) at the 3σ level, as already reported by García-Burillo et al. (2005). Central molecular gas emission has also been detected in HCN at a much more prominent level (Krips et al. 2007). The gas bridge between the ring and the nucleus, identified by García-Burillo et al. (2005) as lying at PA 134°, is also observed in our data cubes. Here we find the bridge to be present between PAs of 78° and 158°, but the most intense emission is along a PA 113°, as shown in the position-velocity diagram in Fig. 3.

The morphology discussed here was for the most part also observed in the previously published PdBI-only maps. The exception is the spatially unresolved resonance ring at the edge of our field, which was not seen at all. The straight gas lanes are more prominent in the SSC maps presented here. The consequences for the gravitational torques will be discussed in Sect. 6.1.

Table 4

CO line fluxes.

Table 5

MH2 masses.

In Tables 4 and 5 we list the integrated CO line fluxes and the derived H2 masses for the different components described above. The regions over which the flux has been measured are indicated in Fig. 1b. These values have been derived with a CO-to-H2 conversion factor XCO of 2.2 × 1020 cm-2 [K km s-1]-1 (Solomon & Barrett 1991) and have not been corrected for helium abundance. As the CO line ratio (Sect. 3.2) is fairly constant in the center, the assumption of a single conversion factor seems valid for our purpose. There are claims that the CO-to-H2 conversion factor is a factor 3 − 4 lower in gas-rich centres of spiral galaxies (e.g. Weiß et al. 2001). A lower conversion factor would lower our mass estimates correspondingly. However, for consistency with other papers in this series, we use the Galactic value. The total integrated 12CO(1−0) flux within a 40″ field of view is 314  ±  55 Jy km s-1. The 12CO(1−0) flux we find, is similar to the measurements of 334  ±  12 Jy km s-1 within 65″ by Kohno et al. (1999, obtained with the Nobeyama Millimeter Array and the 45 m telescope) and 350  ±  41 Jy km s-1 within 45″ by Young et al. (1995, with the FCRAO single dish telescope). The corresponding H2 mass is 1.6 × 109   M. For the 12CO(2−1) emission we measure an integrated flux of 452  ±  65 Jy km s-1 within a field of view of 14″, which corresponds to an MH2 of 4.6 × 108   M if we assume ICO(2−1)/ICO(1−0) = 0.8 (see Sect. 3.2). As can be seen from Table 5 the total measured mass in CO is completely contained in the components mentioned. The areas within which flux is measured are the same for natural and robust weighted maps. In both CO(1−0) and CO(2−1) lines we see that the natural weighted maps are more sensitive to low-level diffuse components of the CO gas. The areas over which the flux of the components has been measured, are marked in Fig. 1.

3.2. Line ratio

The CO(2−1)/CO(1−0) line ratio has been used to convert the observed CO(2−1) flux into equivalent CO masses in Sect. 3.1. We have derived the line ratio (Fig. 4) in the following manner. Both maps, CO(1−0) and CO(2−1), have been short-spacing corrected and sample down to the same minimum uv radius. The robust weighted CO(1−0) data cube was smoothed to the resolution of the CO(2−1) natural weighted cube. Then the zeroth moment map of this smoothed CO(1−0) cube was constructed as before. The CO(2−1) zeroth moment map was regridded to the pixel scale of the CO(1−0) map (from 0.10″ to 0.25″). The fluxes of both maps were converted into temperature and the ratio taken.

The ratio is almost constant along the spiral arms, with the map displaying an average ratio of 0.8. The high ratios found at the edges of the ring are insignificant due to low S/N in the CO maps.

3.3. Kinematics

Evidence of the dominant influence of the large-scale stellar bar is visible in the velocity field. CO(1−0) line emission has been detected in the velocity range of  −220 km s-1 to 250 km s-1 relative to the heliocentric velocity of 1425 km s-1 in the natural weighted data cube. For the CO(2−1) emission, the velocity range is slightly smaller; from –200 km s-1 to 195 km s-1. That is in part due to the higher rms noise in this data cube, which affects detection of the signal in the channel maps at the highest relative velocity offsets. All channel maps with significant emission are shown in Figs. 5 and 6. As the velocity increases, the line emission shifts from the north-west to the south-east. This is consistent with the major kinematic axis of this galaxy, measured by Haan et al. (2009) as 138° using HI data (Table 3). The channels close to systemic velocity (−50 to 100 km s-1) show two maxima, from both sides of the CO ring, as well as extended arms from the bar-driven straight gas lanes.

thumbnail Fig. 5

Channel maps of the naturally weighted 12CO(1−0) data cube. The size of each channel map is 44″ by 44″. The contours are at  − 3σ,  − 2σ, 3σ, 5σ, 10σ, 15σ and 25σ, with 1σ = 2.5 mJy/beam. The velocity relative to systemic velocity of the galaxy (vsys = 1425 km s-1) is indicated in the upper left corner. The dynamical center is indicated by a cross in each channel map. The synthesized beam is given in the lower left corner of each channel and the dirty beam is shown in the lower right panel.

Open with DEXTER

The CO(1−0) velocity map is shown in the top left panel of Fig. 7. The iso-velocity contours in the center are almost perpendicular to the major kinematic axis. As the radius increases the iso-velocity contours bend, forming the “S”-shape distinctive of velocity fields dominated by large-scale bars.

The dispersion (see the velocity dispersion map in the right panel of Fig. 7) reaches values of up to 70 km s-1. These values are high for a gas disk and have led us to further investigate the kinematics of the CO. The data cubes show that the observed CO line emission arises from two distinct components in velocity space at several positions, most prominently where the dispersion maps show high values. This complex velocity structure means that we are not seeing a truly high local velocity dispersion in a single component, but rather the projection of two velocity components within the same beam. The assumption of a single component when we determine the velocity dispersion is clearly wrong for some positions. We see these double peaks in both CO(1−0) cubes and the natural weighted CO(2−1) cube.

In order to quantify the double emission peaks in detail, we fitted double Gaussians along the spectral axis at each spatial pixel of the CO(2−1) natural weighted data cube (Fig. 8). The choice of this cube was made based on the higher spectral resolution of the (2−1) cube with respect to the CO(1−0) cubes. We used the function “XGAUFIT” from GIPSY for the fitting. At the spatial pixels where we have a double peak, we separate the two components based on the large-scale bar model we derive in Sect. 4.2. In Fig. 7 (middle, left) we show the central value of the fitted Gaussian for pixels requiring only a single Gaussian fit, and the central value of the Gaussian fit closest to the bar model for the pixels with a double Gaussian fit (component “V1”). This results in a good representation of the velocity field of the disk of the galaxy. In Fig. 7 (bottom, left) we present the central values of the second Gaussian (component “V2”). For the most part, the double peaks are present in the northern spiral arm/ring region, with the exception of the complex in the south-south-west close to the nucleus. The velocity difference, Δv, between the two components is between 40 km s-1 and 120 km s-1.

thumbnail Fig. 6

Channel maps of the naturally weighted 12CO(2−1) data cube. The size of each channel map is 25″ by 25″. The contours are at  − 3σ,  − 2σ, 3σ, 5σ, 10σ, 15σ and 25σ, with 1σ = 7.8 mJy/beam. All other notation is as in Fig. 5.

Open with DEXTER

Most of the two components might be interpreted in terms of x1 and x2 orbits. x1 orbits are parallel to the large-scale bar (building up the straight gas lanes), and x2 orbits are perpendicular to the bar (supporting the ring). At the resolution of the data cube (beam size 1.72′′ by 1.56′′), the two orbit families might blend together. We find double emission peaks in the CO(1−0) cubes and the natural weighted CO(2−1) cube, where we are unable to resolve the end of the one spiral arm from the other spiral arm. Double velocity components would be the natural consequence of the emission of CO gas on the two orbits families being blended together spatially in our data cubes. In comparison, in the robust weighted CO(2−1) cube we do not find evidence of double peaks in velocity. We find little difference in the width (velocity dispersion) of the Gaussians (Fig. 7middle/bottom, right), except in the northern region. There we see a narrower component, connected to “V1”, with a dispersion of about 30 km s-1, and a wider component, connected to “V2”, with a value of around 40 − 50 km s-1. This region is where we have the northern peak of the “twin peaks” morphology, which persists even in the high resolution CO(2−1) intensity map (in the south we see a ridge). The high dispersion here might be explained by the gas being shocked, which leads to more real, local turbulence/dispersion.

The resolution argument does not explain the double emission peaks found inside the spiral arms, south-south-west of the nucleus. In this region, based on the CO(2−1) map, we expect very little emission, let alone emission arising from two distinct (Δv ~ 80 km s-1) components. The velocity of the first component agrees well with what is expected from the large-scale bar velocity field. The second component’s central velocity (~70−90 km s-1) seems to connect it with the southeastern part of the ring. There is little difference in dispersion between the two components. The strength of the second Gaussian is nearly equal to the first component. This double emission region is not spatially coincident with the gas bridge mentioned earlier. In Sect. 6.2 we will discuss the significance of both the bridge and this second kinematic component in the nuclear region for nuclear fueling.

thumbnail Fig. 7

Top left: first moment map of the CO(1−0) emission. Contours: –200 to 175 km s-1 with 50 km s-1 steps, 0 km s-1 thick gray contour. Kinematic major axis (white dashed line at PA = 138°), orientation of the large-scale bar (dashed-dotted line, PA = 84°). Top right: second moment map of the 12CO(1−0) emission. Contours: 0 to 70 km s-1 in 10 km s-1 steps. The dynamical center is marked by a black cross. Middle/Bottom panels: decomposed velocity field of the CO(2−1) natural weighted line emission. Black spots correspond to blanks due to bad fits at some spatial pixels. Middle left: the central value of the Gaussian fitted for each spatial pixel. If two Gaussians where fitted, the value of the Gaussian closest to the large-scale bar model velocity is given. The white contours represent the CO(2−1) intensity map (5σ in 40σ steps). Bottom left: for the spatial pixels where a double Gaussian was fitted, we plot the central value of that second Gaussian here, i.e. the value further away from the bar model velocity. The size and PA of the nuclear stellar oval and PA of the gasbridge, that will be discussed in Sect. 6.2, are highlighted with a white ellipse and dashed lines. Middle/Bottom right: velocity dispersion corresponding to respective velocity components. The dispersion is the Gaussian fitted σ at each position. Contours from 0 to 70 km s-1. The region with significant differences in velocity dispersion (discussed in Sect. 3.3) is highlighted with a white ellipse.

Open with DEXTER

4. Kinematic modeling

4.1. Motivation

NGC 6951 has a prominent large-scale stellar bar with a radius of 3.0 kpc at a position angle of 84° (Mulchaey et al. 1997). We suggested in Sect. 3 that the gas morphology and kinematics in NGC 6951 can be well explained in terms of this barred potential. The spiral pattern and the high intensity twin peak morphology are both consistent with gas streamlines under the influence of a large-scale bar.

To further investigate this claim, we model our CO observations with a gravitational potential where the only non-axisymmetric component is a large-scale bar. This cannot be done from observations directly since this would assume exact prior knowledge of the non-axisymmetric components of the gravitational potential. At the same time this modeling will also allow us to derive gravitational torques in this region for an independent comparison with the gravitational torque results of García-Burillo et al. (2005) and Haan et al. (2009) (Sect. 6.1).

thumbnail Fig. 8

Six example spectra from the robustly weighted CO(2−1) data cube where we fitted a double Gaussian to the velocity axis. The location of each spectrum is indicated with red lines in the intensity map (middle; same area as Fig. 1 bottom, left panel).

Open with DEXTER

Table 6

Bar potential model parameters.

4.2. Description of the model

We modeled the inner 5 kpc (diameter) of the observed CO(1−0) gas distribution and its underlying gravitational potential using the DALIA modeling software (Boone et al. 2006) that was already used to model the NUGA target NGC 4569 by Boone et al. (2007).

The tool constructs a kinematic model of gas particles in the presence of a barred potential. The potential is built up from two components, an axisymmetric component representing the disk, and a weak m = 2 perturbation: the “bar”. The axisymmetric component has a logarithmic shape and is defined by a characteristic length (rp) and velocity (vp). (1)The bar has the same radial profile as the axisymmetric component, but is tapered by a sine in azimuth and can be oriented in any preferred direction relative to the galaxy’s kinematic major axis (Table 6: perturbation azimuth). It is further defined by the relative amplitude of the perturbation with respect to the axisymmetric component (i.e. the bar strength, ϵ), and its pattern speed. (2)The total kinematic model potential can therefore be described as follows: (3)The model is populated with a distribution of gas particles, described in terms of the radially varying parameters column density, velocity dispersion and scale height. Two dissipation terms, one acting radially and one acting azimuthally, are set to reproduce the dissipative behavior of the gas particles. Finally the model can be inclined and positioned corresponding to the observations, and a rotation sense is set.

Table 7

Bar potential model radial density parameters.

thumbnail Fig. 9

Top, left: deprojected orbit pattern in the inner 3.0 kpc of the bar model. Top, right: resonance diagram. The rotation frequency Ω (red), (blue), (cyan) and the bar pattern speed (yellow dashed). The iILR at 160 pc, oILR at 1150 pc and CR at 4000 pc (outside the displayed range) are indicated with black dots. Overlayed is the rotation curve, which approaches 200 km s-1. The vertical dashed line indicates the outer edge of the models radial mass distribution. Bottom, left: integrated emission of the bar model with contours of CO(1−0) intensity map overlaid. The model emission has been computed at the CO(1−0) resolution. Contours run from 5σ in steps of 20σ (1σ = 0.025 Jy beam-1 km s-1). The red cross indicates the position of the center. Bottom, right: velocity field of the bar potential model with contours of the CO(1−0) velocity field. Contours run from –200 to 200 km s-1 in 25 km s-1 steps.

Open with DEXTER

This model has several limitations as discussed by Boone et al. (2007). First, it is singular at corotation. Buta et al. (2003) and Block et al. (2004) place NGC 6951’s bar radius at 3.0 kpc and indeed, as can be seen from Fig. 9 (top, right), our model does not extend to corotation. Second, closed orbits are computed, which are not self-consistent with the inclusion of two dissipative terms. Therefore, the locations of resonances obtained from this model should be taken as a first-order estimate only.

All parameters in the model can in principle be freely chosen. Our interest here is in obtaining a good set of parameters that reproduces well the observed CO gas distribution. For this reason, we adopted first guesses for the model parameters based on known values of the system and subsequently slightly varied them so that the model better fits our observations. The radial gas mass distribution (Table 7) was deduced from the CO(1−0) emission maps (Fig. 1, top) and can be scaled with the total observed H2 mass. For the velocity dispersion we take a single value, 30 km s-1, which is reasonable for these radii and the assumption of a pressure supported gas disk (Jogee et al. 2005). The rotation sense is set to clockwise, as the dust lanes along the bar are assumed to be leading. The dissipation rates were initially left at default values and only changed after the other parameters had been fixed. Table 6 lists the values of the potential parameters with which the model best reproduces our observations. The best fit was estimated by comparing the modelled channel maps (Fig. 10) and intensity map (Fig. 9bottom, left) against their observed counterparts by eye. In Appendix A we show the agreement of our best fit model with the CO observations at several different positions in the data cubes.

4.3. Best fit model

Figure 9 shows the orbits (top, left) and rotation curve (top, right) for our best model. The orbit pattern shows the change from x1 to x2 orbits, the main periodic orbits in a bar potential, starting at r = 1.5 kpc. The existence of two ILRs, denoted iILR (inner) and oILR (outer) in our model, can be seen from the double intersection of the curve with the pattern speed. Corotation is beyond 2.5 kpc, which means the model does not diverge in the range we are studying here.

The presence of two ILRs in NGC 6951 has also been found by both Rozas et al. (2002) and Pérez et al. (2000), who put them at 180 pc and 900 pc, and 180 pc and 1100 pc, respectively. Our CO observations cannot be fit by the model without 2 ILRs. Their locations are determined by the free parameters of the model; the pattern speed and the shape of the potential (defined by its characteristic length and velocity). From the resonance diagram in Fig. 9 we find our ILRs to be at 160 pc and 1150 pc. These radii are close to the estimates reported by the earlier papers. The orientation of the bar perturbation (130°) in the model is also close to the value for the bar reported by Mulchaey et al. (1997). This degree of consistency shows that the model is a good representation of the observations, and thus we conclude that the CO gas kinematics are indeed dominated by the large-scale stellar bar.

The final model moment maps and channel maps are presented in Figs. 9 (bottom) and 10. The model has been computed at the resolution of the CO(1−0) data. A comparison of Fig. 10 with the observations (Fig. 5) shows good agreement. The model reproduces the straight gas lanes as well as the two peaks of emission in the ring in the individual channels from  − 30 km s-1 to 40 km s-1. The channels with emission in the model are only slightly fewer in number than in the observed channel maps. In the integrated emission map of Fig. 9 (bottom, left), the two spiral arms and the twin peaks are also well reproduced in this model. It is very encouraging that this bar potential model with its simple kinematics and intrinsic limitations fits the observations so well. Therefore, we are confident when we use the model velocity field to separate the two velocity components we detected in the data cubes.

In the channels,  − 30 km s-1 to 40 km s-1, the orientation of the two emission peaks is slightly more north-south in the model channel maps than in the observed one. From Fig. 7 we already saw that when we decompose the velocity field with two Gaussian fits, the velocity contours of the “bar”-component also become more north-south oriented.

In Sect. 3.1 we reported the location of higher density gas traced by CO(2−1) and HCN further downstream in the spiral arms than what is seen in the CO(1−0) line, which is mainly tracing lower-density gas. High emission peaks in both the CO(2−1) and HCN gas distribution occur away from the “twin peaks” that dominate the CO(1−0) distribution. We attemped to reproduce this pattern by also computing a best model at the resolution of the CO(2−1) data. The model gas distribution still agrees with the observations in its main characteristics; spiral arms and the twin peaks. However, the details, such as the southern intensity ridge and the second peak in the northern spiral, are no longer well reproduced. We point out that this model therefore reproduces a morphology closer to the CO(1−0) emission than either of the higher gas density maps mentioned above. A plausible explanation could be that the CO(2−1) and HCN distributions reflect the evolution of the molecular gas inside the circumnuclear ring itself, after the large-scale bar has brought the molecular gas there. As the molecular gas streams within the ring, cloud collapse can occur and/or continue, leading to higher molecular gas densities away from the contact points between ring and gas lanes.

thumbnail Fig. 10

Channel maps of our best model fit at the resolution of the CO(1−0) data. The FOV of each channel map is 44″ by 44″. The velocity indicated in the upper left corner is relative to the central channel (0 km s-1).

Open with DEXTER

5. Model derived results

The model described in Sect. 4.2 can also be used as the basis for a gravitational torque computation. The model’s large-scale bar potential is known at each position, since the model is purely analytical. This makes it easy to derive the gravitational torque on the gas particles at each point, following the equation from Boone et al. (2007): (4)here ϵ is the bar strength, vp the characteristic velocity and rp the characteristic length as set in the model. φ is the azimuthal angle with respect to the bar position angle.

thumbnail Fig. 11

Left: net gravity torque derived from the best fit model. Right: gas accumulation rate as a function of radius. The orbital characteristic radius corresponds to the radius of circular orbits that a given gas particle would maintain if unperturbed by the model’s large-scale bar.

Open with DEXTER

Since the orbits for the gas particles are known, we define the net torque at a given radius as the time-average integrated torque over the orbit with this characteristic radius. We recall that the characteristic radius of an orbit is the radius this orbit would have without a bar (the orbit would then be circular). As mentioned in Sect. 4.2, this approach is not completely self-consistent because when there is angular momentum change the orbits cannot be perfectly closed. Therefore, the gravitational torques computed here should be considered as a first-order approximation to the case in which the orbits are nearly closed. This is a valid approximation if the angular momentum loss is small. Using the simple argument that is of order , where v is the characteristic velocity of the potential v = 200 km s-1, we can see that the relative amplitude of the angular momentum loss per orbit is of order  ~0.0125−0.05. Therefore, it takes at least 20 orbital times for the gas to fully lose its angular momentum.

Figure 11 (left) shows the net gravity torque for each orbit in the model. On the horizontal axis we plot the characteristic radius of each orbit. The net gravity torque is negative over the entire radial range studied here. From 1.9 kpc down to 1.0 kpc the torque stays somewhat constant at values of  − 2000(km s-1)2. There is a slight upturn at larger radii. This is the region of the straight gas lanes. The large value indicates that the bar is easily able to transport gas inward to the gas spiral arms. The torques then increase at an almost constant rate down to a radius of 400 pc. We can understand this from Fig. 9 (top, left): the orbit position angle with respect to the bar stays almost constant (at its maximum value) and the orbit elongation decreases with decreasing radius. A second, smaller, minimum is present at r = 150 pc, at the iILR. The net gravity torque becomes zero within the 50 pc radius. A close inspection at the model orbits in the inner 300 pc shows that the orbit orientation changes again from x2 to x1.

The mass accumulation rate that we can calculate from the gravitational torques also depends, in contrast to the gravitational torque per se, no longer solely on the analytical potential, but also on the radial gas mass distribution in the model. Figure 11 (right) shows the mass accumulation rate for this model. The variation in the mass accumulation rate indicates that gas must accumulate on some orbits and be depleted from others. That the gravitational torques drive the gas inward is also evident from this curve. At radii larger than 800 pc gas is being depleted and the accumulation rate is negative. These are the radii where the gas is located in the straight gas lanes. From 400 pc to 800 pc radius, the accumulation rate becomes positive. These are the radii at which we find the gas spiral arms and the gas ring, which are separated into two peaks in the figure. In the radius range 400 − 1000 pc we already saw from the torque plot that the torques decrease with decreasing radius. As the radius decreases the torques are less strong, so less gas is being moved inwards at each radial bin and the gas can accumulate. Inside 400 pc we see a sharp positive peak, at 150 pc. This is the radius at which we found a local minimum in the gravitional torques. The feature is confined to only a very limited radial range, but indicates that some accumulation may still happening inside the ring.

Mass accumulation within 1.5 kpc in radius, i.e. integrated over this radial range, is occuring at a rate of  ~ +2.0 M yr-1 and mostly in the two spiral arms. From our observations we know that the total amount of molecular gas within this radius is about 1.6 × 109   M. Therefore, at our derived accumulation rate, all matter in the spiral arms must have been brought there within the last 0.8 Gyr. However this is a lower limit to the age of the ring, since it assumes that no gas has been turned into stars and that there is always enough gas available at larger radii to be driven inward.

Errors in the gravitational torque computation can come from six model parameters; the characteristic length and velocity of the potential, the bar pattern speed, the bar strength and the dissipation factors. The first three influence the locations of the resonances as well as the values of the net gravity torques. A 10% difference in either of these parameters can cause up to a  ~ 20% difference in the locations of the gravitational torque minima. The change in value is 5% for the characteristic length and 20% for the characteristic velocity and pattern speed. The net gravity torque scales linearly with bar strength. The dissipation factors only influence the strength of the gravitational torques. A 10% difference in these factors leads to a decrease of  ≈ 25% in the net gravity torques. However, such changes to the model parameters lead to significantly worse models. Therefore, the uncertainties in the results we derive here for the net gravity torque and the mass accumulation rate are within the relative errors given here.

6. Discussion

6.1. Gravitational torques

We can now compare different gravitational torque measurements we obtain using either the “torque map” method introduced by García-Burillo et al. (2005, hereafter: GB05; see also Haan et al. 2009, hereafter: H09), and the analytical computation possible for the large-scale bar model following Boone et al. (2007).

First, we investigate how the inclusion of the 30 m data changes the results obtained with the “torque map” method. Using the code PyPot of H09 we have derived the gravitational torques for our PdBI+30 m data cube. PyPot uses NIR high-resolution images to evaluate the stellar potential. The gravitational torque curves for PdBI-only (dashed line) and PdBI+30 m (solid line) data are shown in Fig. 12. As can be seen from the figure, the inclusion of the 30 m data leads to a change in the gravitational torque budget, especially in the region inside 500 pc. The high positive torques derived at 400 pc in the PdBI-only maps disappear in the PdBI+30 m result. This can be explained by the larger amount of diffuse molecular gas that is recovered in the PdBI+30 m data, which lowers the high torque barrier at  ~ 200 pc found by GB05 (see their Fig. 12) and H09. The spiral arms become more prominent with the inclusion of the 30 m data, so they contribute more to the torque budget at larger radii, shifting the minimum to larger radii. This test shows that sampling the line emission on all spatial scales is important when computing gravitational torques from observations as the absolute values can vary significantly while the overall shape is preserved.

thumbnail Fig. 12

Gravitational torque curve derived from the PdBI-only data ((dashed line) reproduced from Haan et al. 2009). Gravitational torque curve derived from the PdBI+30 m data (solid line) using the same method. The addition of the 30 m data has a large influence on the torque values, especially inside 500 pc.

Open with DEXTER

Now, we turn to the comparison of the “torque map” result with the analytical computation. We are interested in the torques as we want to compute the angular momentum change in the gas (whether the gas is in- or outflowing). In the analytical case, we can simply integrate the net torque over gas orbits from the bar model potential (which is a simplification of the true underlying potential). This integration gives the angular momentum loss of a gas particle in one orbit. In the case of “torque maps” it is often argued that it is impossible to know the exact orbits from observations. So, a statistical argument is used in the “torque map” method. The assumption is that the observed gas distribution can be linked to the time spent by the gas along the orbits. The resulting torque is the gas column density weighted average of the local torques averaged per radial bin.

There seem to be significant differences between the results form both methods. In the “torque map” result we find the strongest (negative) torques at a radius of  ~650 pc. For the analytical computation the most negative torques are approximately at a radius of r ~ 1200 pc. As we saw from the intensity maps (Fig. 1), we observe very little gas in the straight gas lanes (r ≳ 1.0 kpc). We know, however, that the large-scale stellar bar has a radius of 3.0 kpc (Mulchaey et al. 1997). At 1 kpc the large-scale bar still exerts significant torques on the gas, but the “torque map” method being very sensitive to the actual observed gas distribution would not show that since there is so little gas, hence the earlier upturn. There is also a difference in the absolute values of the strongest negative torques between the “torque map” method and the analytical computation. This inequality is due to the difference between computing the net torque (model) and averaging the torques (maps).

In the inner 300 pc, another difference is seen between the “torque map” method and the analytic computation. In our analytical computation we see a second small dip in the torques, while in the “torque map” result we find positive torques. By construction, it is impossible in our model for the torques to change sign. Further, GB05 have found that there is a secondary gravitational component, a nuclear stellar oval, that plays a role in the torque budget here. No such secondary component has been included in our model.

Despite these differences, the comparison between our mass accumulation rate in Fig. 11 and the mass accumulation rate found by H09 (see their Fig. 12) shows little difference. We find values not more than a factor 2 larger than theirs. This agreement shows that our large-scale bar potential model and the subsequent analytical computation give a good representation of the true gas transportation in this region, at least from the outer disk down to the circumnuclear ring.

6.2. CO gas flow inside the ring

Another aim of this paper is to see if a large-scale bar can explain the gas distribution and kinematics over the entire 1.5 kpc radial range. For all radii down to the ring the model represents our observations very well. Inside the ring we found negative torques in our model, but positive torques based on H09’s “torque map” method. Does the influence of the large-scale stellar bar extend further in? As a first test we looked for evidence of gas spiral arms inside the ring, as predicted by the model of Englmaier & Shlosman (2000) and Maciejewski (2004). In Fig. 2 we see some indication of spiral structure in the dust. In the observed CO morphology we do not, nor do we find such evidence when we decompose the velocity field (Fathi et al. 2005; van de Ven & Fathi 2010). The latter method is sensitive to even a small arm/inter-arm density contrast. Based on this outcome, it seems that the CO gas is mainly on stable orbits in this region, which implies that it is not inflowing. However, we did find other evidence of CO gas that is potentially not on stable x2 orbits; the gas-bridge and the double CO emission peaks in the line profiles (Figs. 3 and 7, Sect. 3.3).

The primary component of the double peak emission inside the ring (Fig. 7, middle right) shows velocities indicative of stable orbits. The line-of-sight velocity of the “second peak” CO component (located within the white ellipse in the bottom right panel of Fig. 7) is 70 − 90 km s-1, which corresponds to a velocity difference from the primary CO component at these coordinates of Δv ~ 80 km s-1. It is unlikely that we have a simple superposition of two gas clouds, both on stable orbits. After separating the flux in this area into the two velocity components, we derive an H2 mass for this “second peak” CO component of 1.8 × 106 M.

To understand if this CO component might still be related to gas spiral arms induced by the large-scale bar, the gas flow direction of the “second peak” component and whether it is located in the disk of the galaxy must be ascertained. The first issue is straightforward. As Storchi-Bergmann et al. (2007) have concluded the near side of the galaxy is to the southwest and the far side to the northeast. Thus, the positive velocity of this CO component, which is south-south-west of the nucleus, indicates streaming toward the nucleus if the component is in or in front of the disk, and outflow if it is behind the disk. Outflows are usually biconical. If this component is part of an outflow, seen behind the disk, we would expect to see a second blue-shifted component from the outflow in front of the disk, which is not the case. Nor do we really expect an outflow, since this galaxy does not have a strong AGN or starburst near its nucleus. It is thus unlikely that this CO component is located behind the disk and outflowing.

If the component is inflowing, it is very unlikely that it has any significant scale height above the disk. NGC 6951 is an isolated galaxy that has been dynamically undisturbed for at least 1 Gyr. This history means that the origin of the CO component has to be internal. It could be that some intense star formation in the ring expelled CO from the ring. We do see a gap in the CO ring close to the component. However, the low velocity dispersion of the CO gas in the component makes it unlikely that the CO component was heated due to star formation, became unbound and was pushed high above the disk and is now falling back. Nor do we see evidence of a recent starforming event at that location.

Thus, it seems most likely that the CO component is inflowing and located in the disk, bringing us back to the question of whether the component’s motion can be explained by large-scale bar induced gas spiral arms, even though our initial tests were negative. Storchi-Bergmann et al. (2007) claim to see evidence of two such spiral arms in Hα line emission. The location of their Hα spiral arms does not match the location of our CO component however. Their derived velocity of  ± 20 km s-1 is also at least a factor 3 smaller than that of our CO component. If there are really gas spiral arms in NGC 6951, this discrepancy is very unexpected. CO should be a better tracer of such gas spiral arms, since it is the dominant gas tracer in this region. So it seems that the large-scale bar does not induce further gas inflow in this galaxy and that this CO component is inflowing due to some other reason.

The CO component might be connected to the stellar oval found by García-Burillo et al. (2005) inside the circumnuclear ring. The position angle of this oval is  ~66°. Unlike the Hα spiral arms, its orientation matches well with the location of the CO component. If we draw a line through the dynamical center with the oval’s PA, the CO component lies almost completely at the southern leading side of the oval (Fig. 7, bottom left). A stellar oval is in essence a “fat” bar and can drive gas inward similar to the large-scale stellar bar on larger scales (i.e. the bars-within-bars scenario). The location of the CO component with respect to the oval is correct for this scenario, although we cannot conclude inflow from our gravitational torque results. The reason why the influence of the stellar oval is not more prominent in the gas might be that its bar strength and thus the torque it exerts on the gas is low (García-Burillo et al. 2005,their Fig. 8b) and the amount of gas reaching the influence radius of the oval is limited by the circumnuclear ring and its associated star formation. The CO gas we detected in the gas bridge lies in the region where it would feel positive torques due to the nuclear oval; in the averaged torque budget, as shown in Fig. 12 this gas can dominate and mask the smaller amount of inflowing gas at these radii. Alternatively, viscosity torques may play a more significant role here, as was also dicussed by García-Burillo et al. (2005).

This CO component could be one of the last links in the chain of gas inflow down to the center. The large-scale bar drives most gas down to the ring. Based on the gravitational torques, this is not the smallest possible radius – only the mass accumulation rate becomes small. So, it is possible for gas to move further in, reach the sphere of influence of the stellar oval or experience viscosity torques, which brings it further toward the nucleus. The CO component we observe does not completely reach the dynamical center. However, spatially and kinematically, it connects to the central HCN detection of Krips et al. (2007). The HCN extends over  ~0.5″ around the dynamical center, and our CO component extends to about this radius from the center. The detected HCN also has a detected velocity range of  ± 70 km s-1, with the positive side where the CO component is, which also has a velocity of 70 km s-1. This concentration potentially completes the chain of gas inflow.

7. Summary and conclusions

We have presented high-resolution 12CO(1−0) and 12CO(2−1) maps of the central region of NGC 6951 using the IRAM PdBI and 30 m telescopes. This galaxy exhibits significant indirect evidence for recent gas inflow. It has a large-scale bar, a circumnuclear SB ring, and an AGN. Molecular gas, as traced by CO, is the dominant phase of the ISM in these regions.

Investigation of the CO data cubes shows that the gas distribution can be well explained by gas response to a large-scale stellar bar. We further quantify this scenario by reproducing the observations using a model where the gas is solely responding to a large-scale stellar bar. The model used here has the advantage that it provides direct predictions for the radial motions of the gas. Our best model is able to reproduce the main characteristics of the observed CO morphology and kinematics.

Gravity torques and gas mass accumulation rates that follow from this model were computed and compared to previous gravity torque maps by García-Burillo et al. (2005) and Haan et al. (2009). Their method is complementary to ours and we explain the differences between the two results. As expected, we find that the large-scale bar effectively transports gas inward to r ~ 400 pc, but it no longer dominates the gravitational torque budget on smaller radii.

The stellar oval reported by García-Burillo et al. (2005) is a likely candidate to take over gas transport inside the circumnuclear ring. Detailed investigation of the data cubes revealed double-peaked CO profiles at several positions inside the gas ring. Most notably, we detect a second CO complex near the dynamical center, with a velocity offset of Δv ~ 80 km s-1, which is not on a stable orbit. From simple geometric and physical arguments, we conclude that this cloud is in the disk of the galaxy, is inflowing toward the nucleus, is likely driven by the stellar oval, and is indeed the last step in bringing gas close to the nucleus.


2

Robust weighting here means the weighting function W(u,v) in a uv cell is set to a constant if the natural weight is larger than a given threshold, and W = 1 otherwise.

3

Groningen Image Processing SYstem.

Acknowledgments

T.v.d.L. acknowledges DFG-funding (grant SCHI 536/2-3).

References

Appendix A: Models’ goodness of fit

In Sect. 4.2 we constructed a bar model based on an analytic potential, and derived a corresponding model data cube. The best model was estimated by comparing the model’s channel maps and intensity map against their observed counterparts. In Fig. A.1 we present overlays of the observed and modelled spectra at different positions in the cubes to further show the goodness of fit.

Figure A.1 shows that the model (black dashed lines) follows the observed line emission (solid lines) very well. There are some local differences in the intensity. These differences occur mainly at the locations of the “twin peaks”, (–1.2, 5.8) and (–1.2, –5.5). There the observed emission is larger than we reproduce with the model. A rejected model, with slightly different values for the parameters char. length (210 pc), char. velocity (220 km s-1), bar strength (–0.030) and bar pattern speed (60 km s-1 kpc-1), is shown with grey dashed lines. This model was rejected because it overpredicts the intensity in the ring.

thumbnail Fig. A.1

CO(1−0) emission spectra at different positions in the observed data cube (solid line), the best model cube (black dashed line) and a rejected model cube (grey dashed line). The spatial positions in the galaxy for each spectrum, in arcseconds, are given in the subpanels. The model shows a good agreement with the observations. The positions (–1.2, 5.8) and (–1.2, –5.5) correspond to the “twin peaks” described in Sect. 3.1

Open with DEXTER

All Tables

Table 1

Overview of PdBI data.

Table 2

Overview of PdBI+30 m data.

Table 3

Properties of NGC 6951.

Table 4

CO line fluxes.

Table 5

MH2 masses.

Table 6

Bar potential model parameters.

Table 7

Bar potential model radial density parameters.

All Figures

thumbnail Fig. 1

Top: integrated 12CO(1−0) emission in natural a) and robust b) weighting for the SSC data. The zeroth moment maps have been primary beam corrected. The 12CO(1−0) emission has been integrated from –220 km s-1 to 250 km s-1. Contours run from 5σ in steps of 10σ. The rms value is: a) 1σ = 0.025 Jy beam-1 km s-1 and b) 1σ = 0.022 Jy beam-1 km s-1. Bottom: integrated 12CO(2−1) emission in natural c) and robust d) weighting for the SSC data. The 12CO(2−1) emission has been integrated from –200 km s-1 to 200 km s-1. Contours run from c) 5σ in steps of 20σ, d) 3σ in steps of 10σ. The rms value is: 1σ = 0.039 Jy beam-1 km s-1 for both maps. The red cross indicates the position of the dynamical center in all panels (Table 3). The beam sizes are shown in the lower left corners and correspond to the values listed in Table 2.

Open with DEXTER
In the text
thumbnail Fig. 2

12CO(2−1) robust weighted intensity map (contours from 3σ in steps of 10σ steps, 1σ = 0.039 Jy beam-1 km s-1); overlayed on a J − H color map (greyscale) from HST. White corresponds to higher extinction.

Open with DEXTER
In the text
thumbnail Fig. 3

Position-velocity diagram along a PA of 113° through the dynamical center. Contours are at 2σ, in 1σ steps, with the –2σ and –1σ contours given in grey. A gas bridge is visible between the northern part of the ring (lower left here) and the nucleus (denoted with a cross).

Open with DEXTER
In the text
thumbnail Fig. 4

Emission line ratio CO(2−1)/CO(1−0) in temperature units. This ratio was computed using the natural weighted CO(2−1) and the robust weighted CO(1−0) data cubes, smoothed to the resolution of the former. The high ratio south-east of the nucleus inside the ring is not significant.

Open with DEXTER
In the text
thumbnail Fig. 5

Channel maps of the naturally weighted 12CO(1−0) data cube. The size of each channel map is 44″ by 44″. The contours are at  − 3σ,  − 2σ, 3σ, 5σ, 10σ, 15σ and 25σ, with 1σ = 2.5 mJy/beam. The velocity relative to systemic velocity of the galaxy (vsys = 1425 km s-1) is indicated in the upper left corner. The dynamical center is indicated by a cross in each channel map. The synthesized beam is given in the lower left corner of each channel and the dirty beam is shown in the lower right panel.

Open with DEXTER
In the text
thumbnail Fig. 6

Channel maps of the naturally weighted 12CO(2−1) data cube. The size of each channel map is 25″ by 25″. The contours are at  − 3σ,  − 2σ, 3σ, 5σ, 10σ, 15σ and 25σ, with 1σ = 7.8 mJy/beam. All other notation is as in Fig. 5.

Open with DEXTER
In the text
thumbnail Fig. 7

Top left: first moment map of the CO(1−0) emission. Contours: –200 to 175 km s-1 with 50 km s-1 steps, 0 km s-1 thick gray contour. Kinematic major axis (white dashed line at PA = 138°), orientation of the large-scale bar (dashed-dotted line, PA = 84°). Top right: second moment map of the 12CO(1−0) emission. Contours: 0 to 70 km s-1 in 10 km s-1 steps. The dynamical center is marked by a black cross. Middle/Bottom panels: decomposed velocity field of the CO(2−1) natural weighted line emission. Black spots correspond to blanks due to bad fits at some spatial pixels. Middle left: the central value of the Gaussian fitted for each spatial pixel. If two Gaussians where fitted, the value of the Gaussian closest to the large-scale bar model velocity is given. The white contours represent the CO(2−1) intensity map (5σ in 40σ steps). Bottom left: for the spatial pixels where a double Gaussian was fitted, we plot the central value of that second Gaussian here, i.e. the value further away from the bar model velocity. The size and PA of the nuclear stellar oval and PA of the gasbridge, that will be discussed in Sect. 6.2, are highlighted with a white ellipse and dashed lines. Middle/Bottom right: velocity dispersion corresponding to respective velocity components. The dispersion is the Gaussian fitted σ at each position. Contours from 0 to 70 km s-1. The region with significant differences in velocity dispersion (discussed in Sect. 3.3) is highlighted with a white ellipse.

Open with DEXTER
In the text
thumbnail Fig. 8

Six example spectra from the robustly weighted CO(2−1) data cube where we fitted a double Gaussian to the velocity axis. The location of each spectrum is indicated with red lines in the intensity map (middle; same area as Fig. 1 bottom, left panel).

Open with DEXTER
In the text
thumbnail Fig. 9

Top, left: deprojected orbit pattern in the inner 3.0 kpc of the bar model. Top, right: resonance diagram. The rotation frequency Ω (red), (blue), (cyan) and the bar pattern speed (yellow dashed). The iILR at 160 pc, oILR at 1150 pc and CR at 4000 pc (outside the displayed range) are indicated with black dots. Overlayed is the rotation curve, which approaches 200 km s-1. The vertical dashed line indicates the outer edge of the models radial mass distribution. Bottom, left: integrated emission of the bar model with contours of CO(1−0) intensity map overlaid. The model emission has been computed at the CO(1−0) resolution. Contours run from 5σ in steps of 20σ (1σ = 0.025 Jy beam-1 km s-1). The red cross indicates the position of the center. Bottom, right: velocity field of the bar potential model with contours of the CO(1−0) velocity field. Contours run from –200 to 200 km s-1 in 25 km s-1 steps.

Open with DEXTER
In the text
thumbnail Fig. 10

Channel maps of our best model fit at the resolution of the CO(1−0) data. The FOV of each channel map is 44″ by 44″. The velocity indicated in the upper left corner is relative to the central channel (0 km s-1).

Open with DEXTER
In the text
thumbnail Fig. 11

Left: net gravity torque derived from the best fit model. Right: gas accumulation rate as a function of radius. The orbital characteristic radius corresponds to the radius of circular orbits that a given gas particle would maintain if unperturbed by the model’s large-scale bar.

Open with DEXTER
In the text
thumbnail Fig. 12

Gravitational torque curve derived from the PdBI-only data ((dashed line) reproduced from Haan et al. 2009). Gravitational torque curve derived from the PdBI+30 m data (solid line) using the same method. The addition of the 30 m data has a large influence on the torque values, especially inside 500 pc.

Open with DEXTER
In the text
thumbnail Fig. A.1

CO(1−0) emission spectra at different positions in the observed data cube (solid line), the best model cube (black dashed line) and a rejected model cube (grey dashed line). The spatial positions in the galaxy for each spectrum, in arcseconds, are given in the subpanels. The model shows a good agreement with the observations. The positions (–1.2, 5.8) and (–1.2, –5.5) correspond to the “twin peaks” described in Sect. 3.1

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.