Molecular gas in NUclei of GAlaxies (NUGA)
Key Words.:galaxies: individual: NGC 4569 – galaxies: active – Galaxy: kinematics and dynamics
This work is part of the NUGA survey of CO emission in nearby active galaxies. We present observations of NGC 4569, a member of the Virgo Cluster. We analyse the molecular gas distribution and kinematics in the central region and we investigate a possible link to the strong starburst present at the nucleus. 70 of the M of molecular gas detected in the inner 20 is found to be concentrated within the inner 800 pc and is distributed along the large scale stellar bar seen in near-infrared observations. A hole in the CO distribution coincides with the nucleus where most of the H emission and blue light are emitted. The kinematics are modelled in three different ways, ranging from the purely geometrical to the most physical. This approach allows us to constrain progressively the physical properties of the galaxy and eventually to emerge with a reasonable fit to an analytical model of orbits in a barred potential. Fitting an axisymmetric model shows that the non-circular motions must be comparable in amplitude to the circular motions (120 km s). Fitting a model based on elliptical orbits allows us to identify with confidence the single inner Lindblad resonance (ILR) of the large scale bar. Finally, a model based on analytical solutions for the gas particle orbits in a weakly barred potential constrained by the ILR radius reproduces the observations well. The mass inflow rate is then estimated and discussed based on the best fit model solution. The gravitational torques implied by this model are able to efficiently funnel the gas inside the ILR down to 300 pc, although another mechanism must take over to fuel the nuclear starburst inside 100 pc.
NUGA, for NUclei of GAlaxies (García-Burillo et al., 2003a), is a survey of nearby active galaxies aiming for a deep insight into the mechanisms responsible for feeding the nuclei of active galaxies. High resolution (down to 0.5”) and high sensitivity (rms of 2 mJy beam s for 10 km s channels) observations of 12 galaxies were obtained with the IRAM Plateau de Bure Interferometer (PdBI) and the IRAM 30 m telescope. The CO 1–0 and 2–1 transitions thus mapped allow us to study the dynamics of the molecular gas in the central kiloparsec of the galaxies with a linear resolution as high as 10 pc.
The main results achieved by the NUGA project thus far can be summarized as follows:
A great variety of morphologies and dynamics is observed: 2 and 1-arm spirals (García-Burillo et al., 2003b), bars, rings (Combes et al., 2004), circumnuclear discs, asymmetries, and warps (Krips et al., 2005). As a corollary, there is no evidence for features (at scales 10 pc) that would be typical of active galaxies and that could be uniquely linked to the activity of the nucleus.
We do not see any 2D-kinematic pattern indicative of systematic inflow motions in the circumnuclear disks of some of the NUGA galaxies examined thus far down to the scales of our spatial resolution (10 pc at best). In some cases, e.g. NGC 4826 (García-Burillo et al., 2003b) and NGC 7217 (Combes et al., 2004), the observed perturbations to the disks seem less likely to drive than to inhibit the fueling of the nucleus.
A detailed study of the gravitational torques exerted by the large scale disk asymmetries on the gas in four galaxies showed that such asymmetries may not be the dominant mechanism for angular momentum removal in the nuclear regions (García-Burillo et al., 2005).
These first results suggest either that low luminosity AGN (the NUGA sample includes Seyferts and LINERs) do not require efficient fueling to sustain their luminosity (see e.g. Ho, 2003), or that the timescales for fueling and the onset of activity are so different that both cannot be observed simultaneously (see e.g. Combes, 2004). However, more galaxies in the sample need to be analysed before any definitive conclusions about AGN fueling can be drawn.
NGC 4569 (M90) is a bright SAB(rs)ab galaxy at a distance of 17 Mpc (Bottinelli et al., 1983; Tully, 1988) in the Virgo Cluster. A large scale bar is seen in NIR images (Laurikainen & Salo, 2002) and is almost aligned with the major axis of the galaxy (PA=15 according to Jogee et al., 2005). The galaxy harbours a nucleus of the transition type (type T2 in Ho et al., 1997) which exhibits, by far, the most pronounced nuclear starburst activity among the LINERs and transition nuclei with available UV data (Maoz et al., 1998). As such it became a prime target for studies of the origin of the ionization in transition nuclei. The presence of a supergiant-dominated starburst was established by Keel (1996). Barth & Shields (2000) demonstrated that a starburst with a large population of very hot Wolf-Rayet stars could produce the observed spectrum of the nucleus and AGN activity is ruled out by X-ray (Terashima et al., 2000; Tschöke et al., 2001; Ho et al., 2001) and radio observations (Hummel et al., 1987; Neff & Hutchings, 1992). Gabel & Bruhweiler (2002) showed that for the starburst to be the only source of photoionization (no AGN) extreme conditions are required, with 5 O and B stars squeezed into the inner 30 pc.
On larger scales Tschöke et al. (2001) detected a large soft X-ray emitting region above the disk of NGC 4569 which implies a large energy input into the halo. Hensler et al. (2003) reported on a bipolar outflow seen in X-rays and a giant outflow (10 kpc in length) to the western side of the galaxy in H (see also Bomans et al 2007, submitted). Using the Effelsberg radio telescope at 4.85 GHz and 8.35 GHz, Chyży et al. (2006) discovered large symmetric lobes of polarized radio emission extending up to 24 kpc from the galactic disk on each side of the nucleus. These lobes may have been powered by the nuclear starburst.
H i line observations reveal that the galaxy is H i deficient compared to galaxies of similar type and it must have lost more than 90% of its atomic gas (Warmels, 1988; Cayatte et al., 1990). This was most likely stripped due to rapid motion of the galaxy through the intracluster medium (Vollmer et al., 2004; Boselli et al., 2006).
Interferometric CO observations of NGC 4569 were presented in Helfer et al. (2003), Jogee et al. (2005) and Nakanishi et al. (2005). In this article, higher resolution and higher sensitivity CO(1–0) and, for the first time, CO(2–1) line observations are presented, and the analysis is focused on investigating a possible link between the molecular gas kinematics in the inner kiloparsec and the nuclear starburst. The observations are presented in Section 2. The gas kinematics are modeled based on analytical expressions for the gas orbits. Three different models are compared to the observations starting with a pure geometrical model and proceeding to more realistic models in Section 3. Comparing the different models to the observations allows us to progressively constrain the kinematic properties of the galaxy and emerge with a reasonable fit. The mass inflow due to gravitational torques based on the best fit version of the most realistic model is estimated in Section 4. The goal of this estimate is to quantify the possible contribution of the gravitational torques to the fueling of the nuclear starburst. The results are discussed in Section 5 and summarized in Section 6.
2 CO Observations and global properties of the emission
The galaxy was observed with the IRAM Plateau de Bure Interferometer (PdBI) in 2003 (B and C configurations) and 2004 (A configuration) and with the IRAM 30 m telescope in 2004 in the CO(1–0) (115 GHz) and CO(2–1) (230 GHz) emission lines. The PdBI and 30 m telescope receiver characteristics, the observing setup and the observing procedures are the same as in Combes et al. (2004). The quasar 3C273 was used for PdBI gain calibration. The maps obtained from the 30 m observations are used to compute the short spacings and complete the interferometric measurements using the SHORT-SPACE task in the GILDAS software (Guilloteau & Lucas, 2000) as in Combes et al. (2004). The RMS in 10 km s wide channels is 3.9 and 5.4 mJy beam for natural weighting and the beamsizes are 2.51.6 and 1.10.7 at 115 GHz and 230 GHz, respectively.
2.2 Morphology and mass of molecular gas
Both CO(1–0) and CO(2–1) emission are distributed over a region 176 in size (Fig. 1 a & b) roughly aligned with the major axis of the galaxy (PA=24 deg). The CO(1–0) distribution shows three peaks, one peak close to the center and two symmetrical peaks at 5” from the center and aligned with the overall elongation. The south-western peak is the strongest. The central peak does not coincide with the dynamical center: it appears to be the top of a slightly curved ridge encircling a hole in the CO distribution centered 1 northwest of the dynamical center. Since our observations have a S/N 10 times higher than those presented by Nakanishi et al. (2005), a weaker continuous distribution is detected in addition to the three main peaks, even when going to the highest resolution in the CO(2–1) map.
The integrated CO(1–0) flux within a radius of 20 is 450 Jy km s with an uncertainty due to calibration errors of 15. This is compatible with the fluxes measured by Nakanishi et al. (2005) and Jogee et al. (2005). The corresponding H mass is M using a conversion factor of cm[K km s] (Solomon & Barrett, 1991). 70 of this flux is concentrated within a radius of 10 (800 pc).
2.3 Velocity distribution
The line-of-sight CO velocities span a range of 400 km s. Channel maps are shown in Figures 2 and 3. The most noticeable features are: (1) the distribution is fairly symmetric with respect to systemic velocity and dynamical center; (2) the emission is elongated in most of the channels (from 150 to 150 km s) along a direction close to that of the emission in the integrated map; (3) from one channel to the next as velocity increases, the elongated emission is slightly shifted toward the northwest; (4) the highest velocities (200 km s) are reached at 5” from the nucleus, i.e. close to the two outer peaks in the integrated map; and (5) some low-level emission, mostly seen in the CO(1–0) channel maps (Fig. 2) outside the elongated core, seems to follow the “spider” pattern of rotating disks, i.e. two “arms” on each side of the major axis shift from the north at low velocities to the south at high velocities. To guide the eyes of the reader the two arms are emphasized with thick red lines drawn by hand and overlaid to the channel maps in the bottom panel of Fig. 2. Due to the low level of the signal these arms appear discontinuous and the spider pattern is not seen in the velocity field of Fig. 4 (the first moment mixes this low signal with the noise).
The fact that in individual channels the emission is elongated in a direction close to that of the major axis (PA=24 deg) implies that significant amounts of gas must lie at “forbidden” velocities, i.e. velocities of a given sign appear on both sides of the minor axis, as was noted by Jogee et al. (2005). The velocity field is shown in Fig. 4; it is clearly distorted compared to a disk with only circular rotation. In the central region a steep velocity gradient is running in a SE to NW direction, i.e. orthogonal to the major axis. The kinematic line of nodes is therefore orthogonal to the major axis. This steep gradient results from the continuous shift of the emission in the channel maps (feature (3) described above).
2.4 Line ratio
The CO(2–1)/(1–0) line ratio was computed after tapering the 230 GHz data and restoring it with the same beam as for the 115 GHz data; the ratio map is shown Fig. 1 b. The values are in the range 0.3–1.1 in temperature units. The highest values are reached in the ridge between the central and the southern peaks near the hole in CO. A high line ratio corresponds to high density and/or high temperature gas. A possible interpretation is that the molecular gas in the ridge is being compressed and heated by the nuclear starburst. This interpretation would also be consistent with the presence of a hole in the molecular gas distribution: part of the molecular gas undergoing this pressure might have been dissociated or pushed out of the disk. A fraction of the initial molecular gas in the center must also have been consumed by star formation. The high ratio in the southern peak could be either due to a star forming region or a kinematic shock resulting from the asymmetrically perturbed gravitational potential.
2.5 Comparison to optical observations
A composite color map from B, V and H broad band images and the V-I color map from HST observations are presented in Fig. 5. Dust lanes are clearly seen along the bar, which in NIR observations extends up to a radius of 5 kpc (Laurikainen & Salo, 2002). The bar appears to be aligned with the major axis of the galaxy, like the CO emission, as noted by Jogee et al. (2005). Also striking is the difference in extinction between the eastern and western sides of the galaxy. This confirms that the near side, i.e. the side that appears as the most extinguished because of a low contribution from foreground bulge stars, is the western side. This geometry is consistent with the spiral arms being trailing as in most spiral galaxies. The two outer CO(1–0) peaks coincide with peaks in dust extinction, the southern one being the more prominent both in extinction and CO(1–0) emission. The central CO(1–0) peak does not coincide with any extinction peak, and the hole in the CO(1–0) distribution close to the location of the nucleus has blue optical colors. This agreement suggests that most of the starburst is concentrated within a region of 100 pc in radius, supporting the hypothesis that the hole in the CO(1–0) emission is a consequence of this concentrated starburst having consumed or expelled the molecular gas from the center. The latter interpretation is further supported by an HST H image (Pogge et al., 2000) which shows strong emission south of the hole and several features elongated toward the east (Fig. 1)̧. Some H emission is expected due to the ionizing stars of the starburst but part of it might also trace the impact of the starburst on the ISM and the elongated features to the east are suggestive of an outflow perpendicular to the plane. The east side being the far side, the counter-outflow (i.e. emerging from the other side of the plane) might be obscured and therefore not visible in H, although again visible further to the west as reported by Hensler et al. (2003).
3 Kinematic modeling
3.1 Motivation and strategy
How can the CO emission described in the previous section be interpreted in terms of the radial distribution and kinematics of the molecular gas in the plane of the galaxy? The two outer peaks seen in the integrated map could be due to a ring or to concentrations at the ends of a nuclear bar. The elongated shape seen in the channel maps could result from radial motions (streaming motions inside a density wave or due to shocks at the ends of a bar) or to elongated orbits. To discriminate among these possible interpretations and obtain a quantitative understanding of the deprojected gaseous disk, three different kinematic models are compared to the observations. The approach followed is progressive, starting with a pure geometrical model and proceeding to more ’physical’ models, i.e. with more physical underlying hypotheses used to compute the molecular gas motions. This approach allows us to constrain the physical properties of the galaxy step by step and thus to converge more efficiently to a reasonable solution. This is more efficient because the most physical model is very sensitive to variations of the parameters and therefore more difficult to fit without a priori constraints.
3.2 General description of the modeling
The first model is based on circular orbits with additional radial motions, while the second model is constructed from elliptical orbits. These are purely geometrical hypotheses, and reality is expected to be more complex. However, the symmetry of the CO emission in the datacube (i.e. with respect to the dynamical center and the systemic velocity) suggests the gas flow is well organized and dominated by 2 modes. Thus, if we assume a pure 2 symmetry, comparing the observations to the first model can help to quantify the amplitude of non-circular motions and estimate to what extent they can be reproduced by pure radial motions. Comparison to the second model allows us to evaluate the alternative interpretation of non-circular motions, namely the possibility that the clouds populate near-to-closed orbits that are elongated. The third model is based on analytical solutions for gas particle orbits in a barred potential.
The first two models are entirely specified by a radially varying set of properties similar to that of the widely used ’tilted ring model’ (Rogstad et al., 1974). The values of the parameters are specified at some galactocentric radii and linearly interpolated between those radii. The number of radii where the parameter values are specified and their values can be chosen freely, giving us the ability to reproduce any radial distribution to a good approximation. Simulated data are produced by randomly assigning gas particles to the disk according to the distributions specified, projecting them in the datacube and convolving by the synthesized beam. As most of the emission comes from within the PdBI primary beam full width at half maximum ( at 115 GHz), the effect of primary beam attenuation is neglected. For more details about the fitting procedure see Appendix A.1.
3.3 Axisymmetric model
|rotation vel. (km s)||40.0||30.0||40.0||50.||70.0||100.0||110.0|
|radial vel. (km s)||0.0||0.0||0.0||-120.0||0.0||0.0||0.0|
|col. dens. (no units)||6.0||1.0||0.5||4.0||0.4||0.4||0.4|
|velocity disp. (km s)||40.0||10.0||10.0||20.0||10.0||10.0||10.0|
|scale height (pc)||20.0||15.0||10.0||10.0||10.0||10.0||10.0|
The parameters of the axisymmetric model (see details in Appendix A.2) reproducing the main features of the observations with 7 radii are listed in Table 1. The corresponding channel maps for 5 velocities are shown in the second column of Fig. 6, while the corresponding integrated intensity and velocity maps are shown in Fig. 7 b. A face-on view of the model is shown in Fig. 8 a.
The three peaks in the integrated map are reproduced by a ring with mean galactocentric radius of 600 pc and a central concentration. The column density must be higher by a factor of 10 inside the ring (at 600 pc) with respect to the density at 300 pc and 700 pc. The column density must increase again inside a radius of 300 pc to reach a value higher by a factor of 15 at the center. The distribution of the CO emission in the channel maps requires an inward radial velocity of 120 km s in the ring at 600 pc. Lower radial velocities in the ring fail to reproduce the data.
This model can reproduce the shape of the CO emission in the channel maps at all velocities. The slight shift of the emission toward the northwest from channel to channel and, as a consequence, the twisted isovelocity lines are also well reproduced (compare Fig. 7a with 7b). Moreover, the ring interpretation provides a natural explanation for the fact that the elongated shape is aligned with the major axis of the galaxy: it is merely an inclination effect. The width of the emission along the minor axis, however, is larger than in the observations, and the two symmetrical peaks are not as highly contrasted as observed. Also, the radial velocity in the ring required to reproduce the observations is extremely high and cannot be interpreted as the manifestation of real radial motions.
We conclude from this first modeling that, although a ring of 600 pc radius with a central concentration can explain to some extent the spatial distribution of the CO emission, the amplitude of the non-circular motions is unrealistically high, i.e. comparable to the amplitude of the rotation motions. No physics could explain the collapse of a ring at this speed . This result implies that pure radial motions alone cannot explain the departures from circular rotation, which must be at least partly due to the shapes of the orbits themselves. Paradoxically, while the best fit solution of the circular model reproduces the observations well, it implies that the model is not a good representation of reality and that elongated orbits must be considered.
3.4 Elliptical orbit model
|charact. radius (pc)||1.0||85||170||510||620||680||850||1700||2800||5500|
|charact. rot. vel. (km s)||40||25||10||20||40||90||120||150||200||220|
|position angle (deg)||110||130||160||210||160||100||20||-180||-230||-330|
|col. dens. (no units)||3.0||2.0||0.2||0.7||0.8||0.3||0.03||0.03||0.03||0.02|
|velocity disp. (km s)||30||35||40||40||40||30||10||10||10||10|
|scale height (pc)||20||18||17||15||10||10||10||10||10||10|
Although elongated orbits are not expected to be elliptical in general, ellipses provide a first order approximation to any non-circular closed or nearly-closed shape. To estimate a possible contribution from the orbit elongation to the non-circular motions and to get a picture of the orientations of the orbits in the galactic plane (the orbit pattern) we fit a model in which the gas particles populate elliptical orbits and their angular momentum is constant over the orbit. “Nearly-closed” orbits means orbits that are stable on the time scale of a rotation period, which is possible only if the orbits do not intersect, since otherwise the cloud shocks would re-orient the gas on different orbits. Thus, for this model to be self-consistent, the variations in position angle and ellipticity must satisfy this ’no crossing’ condition. The density of clouds on a given orbit is not uniform. If it were assumed independent of the rotation sense it would be inversely proportional to the velocity. However, it is known that molecular clouds are formed in density waves or shocks, from which stars form that then destroy the parent cloud. The density is therefore expected to increase abruptly when entering the density wave or shock and to decrease smoothly after. In this simple model, density waves correspond to overdensities due to orbit crowding, i.e. are kinematic waves: they form at the aphelia of the ellipses when the position angles change with radius. To mimic the formation/destruction of the clouds we apply an exponential taper to the density. The application of this taper significantly improves the model fit. More details about this model are given in Appendix A.3.
The parameters of such a model best reproducing the observations with 10 radii are listed in Table 2. The corresponding channel maps for 5 velocities are shown in the third column of Fig. 6, while the corresponding integrated map and velocity field are shown in Fig. 7 c. The face-on view of the model solution is shown in Fig. 8 b and the orbital pattern in Fig. 9 a. Strong variations in the column density, the angular momentum, the orbit orientations and the velocity dispersion are required at a radius of 600 pc to reproduce the peaks in the integrated map and the shape of the distribution in the channel maps. As in the previous model the column density must increase by an order of magnitude at 600 pc, however inside the 600 pc radius the column density only drops by a factor of 4. A central concentration is still required within the innermost 150 pc. The angular momentum must decrease abruptly inside 600 pc. The velocity dispersion must be 40 km s at 600 pc. The most remarkable result is that to reproduce the elongated core in the channel maps, it is necessary to invert the variation of the orbit position angles at the same radius of the column density peak (in the range 500–600 pc). The orbit position angle must be maximum at this radius and equal to . This feature of the orbit pattern seems to be strongly constrained by the observations.
From Fig. 6 it can be seen that the low level emission is well reproduced by this model, giving us confidence that this emission is real and that the rotation velocity of the disk at 5.5 kpc is 220 km s. This value is consistent with the estimates of Guhathakurta et al. (1988) based on H i observations and Rubin et al. (1999) based on optical observations.
The main differences between the observed and the modeled channel maps occur at high velocities, with the model being unable to reproduce the emission closest to the center seen up to 150 km s. Moreover the model cannot reproduce the continuous shift of the emission to the northwest and therefore the central velocity gradient present in the velocity field.
In addition, to obtain emission that is spatially elongated up to high velocities, we must set a velocity dispersion of at least 40 km s in the ring. This is a high velocity dispersion, a fraction of this may come from a shock, but beam averaging of the velocities at the ends of the orbits can also be expected to contribute significantly. This consideration implies that the real orbits must have sharper ends than ellipses.
The orbit pattern and the face-on view of the emission strongly resemble gas kinematics in a barred potential as observed in many galaxies. The large gradient of the orbit position angle at 600 pc mimics a tightly wounded spiral forming a ring. The change in sign of the gradient of the position angle is a remarkable feature that is known to occur at the Inner Lindblad Resonance (ILR) when there is only one, or between ILRs when there are two of them (see e.g. Wada, 1994). The fact that this feature coincides with a column density peak, which is expected at resonances, indicates that this is actually a single ILR. The two symmetric peaks in the face-on integrated map and the large velocity dispersion at the resonance would then correspond to the shocks frequently observed at the locations of dust lanes in bars, which are due to the crowding of orbits that favors cloud collisions (see e.g. Combes, 1996). Identification of the ILR is further supported by the orientation of the bar inferred from the model, which is similar to that of the large scale stellar bar. Indeed, the shocks are expected to lead the bar (see e.g. Combes, 1996), which would therefore have a position angle somewhere between 0 and 30 (see Fig. 8 and see also the maximum value for the orbit position angle in Table 2: 210 ). This is consistent with the large scale stellar bar position angle of 15 (Jogee et al., 2005). This agreement implies that the resonance observed is unlikely to result from the action of an inner bar, but is rather the ILR of the large scale bar.
Is this interpretation supported by a model based on gas orbits in a barred potential? Why is the velocity field with the strong velocity gradient twist in the center not reproduced? Is this due to the elliptical approximation of the orbits or are pure radial motions (as opposed to radial motions intrinsic to non-circular but closed orbits) also contributing? To answer these questions, we compare the observations to a model based on analytical solutions for the gas orbits in a barred potential.
3.5 Barred potential model
|galaxy position angle (deg)||20|
|characteristic length (pc)||400|
|characteristic velocity (km s)||220|
|bar position angle (deg)||20|
|pattern speed (km s kpc)||60|
|radial dissipation rate (km s pc)||0.002|
|azimuthal dissipation rate (km s pc)||0.005|
|charact. radius (pc)||10||250||400||600||700||1500||3000|
|col. dens. (no units)||1.5||0.1||0.5||1.0||0.05||0.02||0.02|
|velocity disp. (km s)||35||10||30||30||10||10||10|
|scale height (pc)||10||10||10||10||10||10||10|
value fixed by the elliptical orbit model solution
This model was introduced and developed by Sanders & Huntley (1975), Lindblad & Lindblad (1994), Wada (1994), Sakamoto et al. (1999) and Baker (2000). It is based on analytical solutions for the orbits of gas particles obtained by a perturbation calculation of the epicyclic stellar motions. The potential is described by an axisymmetric component and a weak 2 perturbation. The perturbation has the same radial profile as the axisymmetric component and is modulated by a sine in azimuth. The amplitude of the perturbation relative to the axisymmetric component is defined as the strength of the bar. For the radial profile of the potential we use the logarithmic shape. This choice is justified by the ratio of the resonance radii. Indeed, the corotation radius of the large scale bar, assumed to be 1.2 times the bar extension, is 5 kpc according to Laurikainen & Salo (2002) and according to the elliptical orbit model the ILR radius is 600 pc, implying 8. As shown in Appendix B, for a logarithmic potential this ratio is 6.3. For comparison a Plummer potential gives a ratio of 3 (Wada, 1994).
As in Baker (2000), two dissipation terms are introduced that characterize the damping of radial and azimuthal oscillations and are assumed to reproduce the dissipative behavior of gas particles. These terms cause the major axis of the orbits to lead the bar, the offset being maximal at the ILR. The limitations of the model are (1) it is singular at corotation; (2) the epicyclic motions used to compute the perturbations are based on the axisymmetric component only so that the model is not fully self-consistent; and (3) computing closed orbits assuming dissipation rates is not self-consistent, because dissipative particles imply open orbits. As a consequence, the radii of the resonances obtained with this model (Appendix B) cannot be taken as accurate results. However, this can be seen as a first order approximation that works only when the real orbits are near to being closed. In this case the resonance radii obtained give a basic picture of the scales involved. The orbit patterns produced by this model exhibit some of the main features expected in barred potentials, e.g., the elongated orbits, the phase shift change at an ILR, and the spiral pattern. Also, they match observations and results of numerical simulations (see Lindblad & Lindblad, 1994). This analytical approach should be particularly suitable for NGC 4569 which, as noted in Section 2, is rather symmetric and dominated by an 2 component in CO and shows a weak stellar bar in the NIR (Laurikainen & Salo, 2002). The results of the elliptical orbit fit allow us to fix some of the parameters describing the potential. More details about this model are given in Appendix A.4.
The best fit parameters for 7 radii, including the parameters fixed by the results of the elliptical orbit fit, are listed in Table 3. Some of the channel maps of the simulated data are shown in the last column of Fig. 6; the integrated map and velocity field are shown in Fig. 7 d. The face-on view of the model solution is shown in Fig. 8 c and the orbital pattern in Fig. 9. The best fit bar position angle of 20 is consistent with the large scale stellar position angle of 15 derived by Jogee et al. (2005). This result confirms the conclusion drawn from the elliptical orbit model, namely that the kinematics inside the inner kiloparsec are driven by the gravitational potential of the large scale stellar bar111Note that although this conclusion drawn from the elliptical orbit model was used as an hypothesis to fix some of the parameters of the model, these fixed parameters are not related to the bar position angle, so the argument is not circular..
For this model as for the others, a column density jump of an order of magnitude is required at a radius of 600 pc and a central concentration of similar column density is required. A velocity dispersion of 30 km s is required at the radius of the column density increase. Although this is still a high value, it is less than for the elliptical orbit model because the orbits have much sharper ends (Fig. 9) that contribute significantly to the velocity dispersion through beam averaging. This suggests that the orbit shapes in this model are more realistic than the ellipses.
As can be seen in Fig. 6, this model reproduces the data better than the elliptical orbit model, especially at high velocities. The low level emission is not so well reproduced but this is expected as this model diverges at corotation, making it unreliable at large radii. The velocity field (Fig. 7) is also the closest to the data when compared with the two other models discussed above.
The fact that an analytical model that includes certain physical assumptions (about the potential shape and gas properties) reproduces the data well lends support to its underlying assumptions. In particular, the interpretation of the column density increase at 600 pc as an ILR of the large scale bar is well confirmed by this fit. This also supports the assumption that there is only one ILR or that the ILR region is small with little space for x orbits. If there were two ILRs, the gradient of the orbit position angle would change in sign at a similar radius (in the middle of the 2 ILRs as described in Wada, 1994), but the maximum position angle of the orbits would be much greater, leading to a much greater position angle for the oval ring unlikely to fit the observations equally well. The good fit also shows that even though this model is not fully self-consistent, the orbit pattern and orbit shapes favored by its best-fitting version may be reasonably close to reality. We note that the configuration of the orbit pattern implies that the velocity field is very sensitive to the orbit shapes. Indeed, the most elongated orbit are almost aligned with the major axis (the bar is almost aligned with the major axis), and the strongest velocity components are therefore orthogonal to the line of sight. As a consequence the velocity field is not only dominated by the elongation of the orbits but it is also sensitive to the orbit shapes, i.e. how the velocity vectors behave at the ends of the orbits. This can explain why the diffence between the ellipses and these more realistic orbits implies a significant difference in the velocity field.
4 Mass inflow
Two mechanisms contribute to angular momentum dissipation and mass inflow: gravitational torques, which depend on the potential and on the individual gas orbit shapes and orientations, and viscosity torques, which depend on the cloud-cloud collision rate and on the interstellar turbulence. In the previous section the observed molecular gas kinematics were successfully modeled by analytical solutions for gas particle orbits in a barred potential. We now address the question of the mass inflow based on this model solution.
As the potential of the model solution is known at each point in the galaxy, the gravitational torque can be directly computed for each particle following Eq. 12. In addition, as the orbits are known, the torques can be integrated over the orbits to estimate the angular momentum loss over a rotation (Eq. 14) as well as the corresponding radial velocity of the matter inflow (Eq. 16). As already outlined, this approach is not fully self-consistent, since the orbits should not be closed if the angular momentum changes over a rotation, but the angular momentum loss thus computed can be considered as a first order approximation to the case where orbits are nearly closed, i.e., the angular momentum loss due to the gravitational torque is a small fraction of the angular momentum. The good fit of the model to the observations suggests that, on average, the molecular gas in the inner kiloparsec lies in this regime.
The torque is proportional to the bar strength, which is only poorly constrained by the ’weak bar’ regime requirement (see discussion above), to . This implies that only upper limits can be derived for the mass inflow. Figure 10 shows the net torque (as defined Eq. 12) for 0.05 on the top left, and the corresponding radial velocity on the top right. Both curves look similar, because the radial change of the angular momentum as determined by the axisymmetric potential is almost constant above 500 pc (see Eq. 16). Due to their circular shape, the orbits with a characteristic radius less than 300 pc undergo no net torque and as a consequence the large scale bar would not be able to fuel the nuclear starburst on these scales via gravitational torques only. The maximum radial velocity is of order 10 km s and occurs for the orbit with a characteristic radius of 650 pc, i.e. at the ILR. The value of this upper limit for the radial velocities is consistent with our working hypothesis, namely, that to first order, radial motions can be neglected compared to rotational motions, consistent with the results of Wong et al. (2004). The variations of the radial velocity, shown on the bottom left of Fig. 10, imply that the radial flow is not steady over the disk and as a consequence that matter must accumulate in some orbits and be depleted in others. The actual mass flow also depends on the gas mass available at each radius (the radial velocity would directly translate into mass flow only if the gas mass did not depend on the radius, i.e., for a column density profile following ). Nevertheless, without any assumption about the radial density profile, the variation of the radial velocity shows that the gas will tend to be depleted just outside the ILR (the inward velocity increases with decreasing radius down to 600 pc) and to accumulate inside the ILR (the inward velocity decreases with decreasing radius inside 600 pc). This model is therefore consistent with the formation of rings at ILR as observed in many spiral galaxies (Buta & Combes, 1996) and often used as a signature of the ILR.
The mass accumulation rate per radial element at a given radius is the gradient of the radial momentum (product of the radial linear density with the radial velocity); it is plotted as a function of the characteristic radius for 0.05 in Fig. 10 d. To obtain the correct mass units, the column density profile of the model was normalized so that the total mass of gas equals that measured in Section 2. We see that the maximum accumulation rate inside the ILR is of order 0.2 yr pc. Outside and close to the ILR the depletion rate is of similar amplitude. The accumulation rate integrated inside the ILR (i.e. over the radius range 300–600 pc) gives an inflow rate of order 20 yr. This is significantly higher than other inflow rates estimates in the literature which are in the range 1–4 yr (Quillen et al., 1995; Laine et al., 1999; Regan et al., 1997). But this is an upper limit for this specific model and lower values are not excluded. The accumulation of molecular gas inevitably leads to a density increase and star formation. If the fraction of infalling gas transformed into stars (star formation efficiency) is constant and is of order 10, then the star formation rate (SFR) integrated over the accumulation region would be of the order of a few solar masses per year, which is realistic when compared to what is observed at ILRs in other galaxies. In NGC 4569 however, the H morphology and the CO(2–1)/(1–0) line ratio suggest that most of the star formation observed in this galaxy is taking place inside the inner 100 pc (see Section 2.5). The presence of 5 O and B stars in the inner 30 pc (Gabel & Bruhweiler, 2002) would imply a constant SFR of order 2 yr (assuming a Salpeter IMF and a life time of 5 Myr for O and B stars). This implies either that another fueling mechanism is taking over inside the inner 300 pc, or that the nuclear starburst is ending and star formation at the ILR will dominate in the future.
Viscosity torques are intrinsic to the model and are specified through the radial and azimuthal dissipative parameters. However, as outlined in Baker (2000), these parameters lead to closed orbit solutions, so they cannot be directly used to compute the dissipation of angular momentum and the mass inflow rate. Estimating the intensity of the viscosity torques based on the density and velocity dispersion at a given point in the disk would require adopting prescriptions for the visosity that are generally acknowledged as very uncertain. Thus, our model does not allow us to directly constrain the rate of mass inflow due to viscosity torques. It could well be that these torques take over the fueling of the nucleus inside the 300 pc where the gravity torques are not efficient anymore. Alternatively, an inner stellar bar could drive the gas further inward as in NGC 6946 (Schinnerer et al., 2006); however, the current observations do not show any evidence of such a small inner bar.
Models of gas dynamics in galaxies come in two main categories: ballistic models in which the gas is represented by colliding clouds on ballistic orbits (Combes & Gerin, 1985; Byrd et al., 1994), and hydrodynamic models in which the gas is assumed to be a continuous dissipative medium (Athanassoula, 1992; Piner et al., 1995; Wada & Habe, 1995; Friedli & Benz, 1993). Comparisons of the two approaches are discussed in Lindblad & Lindblad (1994), Wada (1994), Combes (1996) and Regan et al. (1999). The analytical model of orbits in a barred potential presented in the previous section effectively traces the evolution of gas particles in the gravitational potential as they are subjected to collisions. Although the equations are continuous, they are not based on hydrodynamic assumptions, and the gas is assumed to be an ensemble of independent particles rather than a continuous medium, meaning the model can be classified as ballistic. The good fit of this model to our CO observations presented here seems to validate such a ballistic approach. The results of the two approaches may actually not have so many differences that the observations can distinguish. As already suggested by Lindblad & Lindblad (1994), it may be that the orbits out of this analytical model describe a gas flow similar to the one obtained with a hydrodynamic model. This hypothesis is supported by comparing the orbit pattern as shown in Fig. 9 b to the gas flow resulting from a hydrodynamic simulation as shown in Fig. 9 of Regan et al. (1999). The shapes of the orbits resemble the gas streamlines in their model and the dust lane/shock region of their model is well reproduced by the orbit crowding in our model. The main difference is that the largest circular orbit in our model is located well inside the ILR, whereas it defines the ILR in their model. As a consequence, the ring formed at the ILR by gas accumulation is elongated in our model (Fig. 8 c) rather than circular. But such a shape has also been seen by Nakanishi et al. (2005) based on Smoothed Particle Hydrodynamics (SPH) simulations of the same galaxy.
The difference between ballistic and hydrodynamic models outlined in Regan et al. (1999) corresponds to the sharp inversion in the radial component of the velocity vector in the shock/dust lane region. This feature is argued to be reproduced by the hydrodynamic models but not by the orbits of ballistic models. However, this argument is based on the assumption that the orbits are elliptical (see their Fig. 10). As emphasized in Section 3.5 the analytical orbits of gas clouds in a barred potential used here have much sharper ends than ellipses, blurring the nominal distinction between hydrodynamic and ballistic models.
It is also interesting to note that the conclusions of both approaches in terms of mass inflow are not so different. The inflow stops at the largest circular orbits (the “nuclear ring”) in the model of Regan et al. (1999) and at the largest x orbit or at the largest nonlooping x orbit in Regan et al. (2004). In all cases, the bars seem to be incapable of driving molecular gas inside a radius of 100 pc. Similarly, in our model, the gravitational torques stops where the orbits become circular at a radius of 300 pc.
The advantage of the method followed here is that we are able to first identify the ILR based on a kinematic argument, namely the inversion of the radial variation of the ellipse orientations in the second model (Section 3.4). This allows us to fix the potential parameters of the third model, i.e. the orbit pattern, and to fit the radial profile of the column density and velocity dispersion independently. Although limited in its domain of application (weak bars and inside corotation), this method has the advantage of being insensitive to initial conditions or numerical effects that may affect evolutionary models (ballistic or hydrodynamic). In addition, it is obviously very undemanding in terms of computation resources and well adapted to iterative fitting to observations.
The main results of this article can be summarized as follows:
In NGC 4569, 70 of the M of molecular gas detected in the inner 20 is concentrated within a radius of 800 pc and distributed along the large scale bar, with two peaks at 500 pc from the center and one peak close to the center. There is a hole in the CO distribution at the location of the nuclear starburst where evidence of outflow is seen in H emission.
Fitting a model based on elliptical orbits reveals a change in sign of the radial variation of the orbit orientations at a radius of 600 pc, where molecular gas is also accumulated. This coincidence, together with the value of the orbit position angle maximum, leads us to identify this feature with confidence as a single ILR of the large scale bar.
The kinematics are well reproduced by a model based on analytical solutions for the gas particle orbits in a weakly barred potential with a single ILR at 600 pc.
The gravitational torques implied by this model are able to efficiently funnel the gas inside the ILR down to a radius of 300 pc but another mechanism must take over to fuel the nuclear starburst, i.e., there is no evidence for ongoing fueling of the starburst by the large scale bar. If such a mechanism does not exist, the nuclear starburst will stop and star formation at the ILR will dominate with a star formation rate of order one solar mass per year (assuming a 10 star formation efficiency).
This last point agrees with the result of García-Burillo et al. (2005), namely that the gravitational torque caused by kiloparsec scale perturbations of the potential seems not to be the only mechanism responsible for fueling galactic nuclei.
Appendix A Detailed description of the modeling
a.1 All models
The first two models are entirely specified by a radially varying set of properties similar to that of the widely used ’tilted ring model’ (Rogstad et al., 1974). The values of the parameters are specified at some galactocentric radii and linearly interpolated between those radii. The third model has additional global parameters. Thus, the total number of free parameters is adjustable through the number of radii, , and is given by , where is the number of properties specified at each radius including the value of the radius. The number of radii required to reproduce the data satisfactorily was found to be greater than 7. In this case, for all models the number of free parameters is too high compared to the dynamic range of the data (signal-to-noise ratio, spatial dynamic range and spectral dynamic range) to permit a rigorous minimization of an objective function like over the full parameter space. Since the goal of this modeling is to test the compatibility of the different hypotheses with the main features of the data, rather than to estimate to a high accuracy the numerous free parameters, the comparison of the simulated data to the real data has been done by eye222This comparison uses the DALIA interface to numerical simulation codes which is being developed jointly by the Observatoire de Paris and the Max Planck Institut Für Radioastronomie (Boone et al., 2006). This graphical user interface gives access to the model parameters and allows the user to directly compare the outputs of the code with real observations..
In all the models, the inclination of the galaxy is fixed to 70 (Tully, 1988), its position angle to 30 (Guhathakurta et al., 1988, except in the third model where the PA is also a free parameter) and its distance to 17 Mpc (Bottinelli et al., 1983; Tully, 1988, only used to convert angular scales into physical scales). In addition, in all models we make the effective assumption that mass column density is proportional to the CO intensity with the same constant of proportionality everywhere.
a.2 Axisymmetric model
In the first model, the properties specified are the radii, the rotational velocity, the radial velocity (positive in outward direction), the column density through the disk, the disk scale height (FWHM of the Gaussian profile) and the local velocity dispersion (square root of the Gaussian distribution variance). The number of free parameters is therefore . The vertical distribution of matter and the velocity distribution at a given position are assumed to be Gaussian.
a.3 Elliptical orbit model
To mimic the formation/destruction of the clouds we apply an exponential taper to the column density along the orbits. The taper equals 1 at the aphelion and decreases like , where is the azimuth along the ellipse. The sign of is determined by the sense of rotation of the cloud in the reference frame of the density wave for which a pattern speed needs to be specified. Here we assume that all the gas lies within the corotation radius, and that the density decreases smoothly in the sense of rotation of the whole galaxy, i.e. counter-clockwise.
The properties of the model are the characteristic radius, the characteristic rotation velocity (the angular momentum), the ellipticity, the position angle of the ellipse, the column density, the velocity dipersion and the scale height. The pattern speed of the entire model was not well constrained by the observations and since the larger the value the worse the fit, it was decided to fix this parameter at zero. The number of free parameters is therefore . As in the previous model, the vertical distribution of matter and the velocity distribution at a given position are assumed to be Gaussian.
a.4 Barred potential model
We use 2 parameters describing the logarithmic potential (the characteristic scale and velocity, see Eq. 1), 3 parameters for the bar perturbation (the position angle relative to the galaxy major axis, the strength, and the pattern speed) and 2 parameters for the gas dissipative properties (the radial and azimuthal dissipation rates). We found it difficult to obtain a good solution with a galaxy position angle of 30 as in the previous models, so we left this as an additional free parameter here. In addition, the gaseous disk is described by 3 properties (the column density, the velocity dispersion, and the scale height) at different radii as in the previous models (the radius here is the characteristic radius of the orbit, which corresponds to the radius of the circular orbit without bar perturbation). The number of free parameters is therefore .
The results of the elliptical orbit fit allow us to fix some of the parameters describing the potential. First, the characteristic velocity of the axisymmetric component of the potential corresponds to the rotation velocity at infinity. Based on previous results, we take km s, also consistent with the value of Guhathakurta et al. (1988) based on H i observations. Next, the fact that the ILR was identified at 600 pc from the center in the elliptical orbit model allows us to fix the characteristic scale of the potential (using Eq. 5) to pc. Finally, assuming there is only one ILR (as the elliptical orbit model suggests), the pattern speed is set by Eq. 6 to = km s kpc. The remaining global parameters (as opposed to those specified at each radius) are the bar strength, its position angle and the 2 dissipation rates. There is a balance between the bar strength and the dissipation rates, in the sense that the stronger the bar or the smaller the dissipation rates the closer the orbits are to pure stellar orbits. Stellar orbits are organized in families, the so-called x and x orbits, (Contopoulos & Grosbol, 1989; Sellwood & Wilkinson, 1993), which are orthogonal to each other, intersect, and can exhibit large loops. As long-lived gas orbits cannot intersect or have loops, a realistic solution for the orbits must keep the bar strength and dissipation rates balanced. The bar strength is required by self-consistency to be in the weak perturbation regime (less than 5) but below this threshold there is a degeneracy with the dissipation rates. The ratio of the two dissipation rates has little effect on the orbits, i.e. there is a degeneracy between them. So these three parameters (bar strength and the two dissipation rates) are not well constrained and only their combination matters.
Appendix B Resonance radii and pattern speed in a logarithmic potential
A logarithmic potential is defined as:
where and are the characteristic velocity and length of the potential. The rotation frequency is given by:
and the epicyclic frequency by:
The inner Lindblad resonance (ILR) occurs when equals the pattern speed of the bar. If we assume only one ILR, this happens when is at its maximum. The radius of the ILR, , is therefore obtained by solving
It scales with the scale radius of the potential, , and is given by
where is the ILR radius in units of .
The pattern speed of the bar, , is given by :
The corotation resonance occurs at the radius where the rotational frequency equals the pattern speed of the bar. So, the corotation radius, , is related to through
which leads to:
where is the corotation radius scaled to . Replacing by Eq. 5 gives:
Appendix C Mass inflow in a logarithmic potential with a bar perturbation
The non-axisymetric component of the potential (the bar component) is given by
where is given by Eq. 1 and is the strength of the bar. The resulting gravitational torque exerted on a particle at is
We define the net torque, , as the average torque over an orbit
where is the orbital period. The angular momentum variation over this orbit and per unit mass is then given by
A particle that changes angular momentum will settle on a new orbit which has an intrinsic angular momentum (as determined by the axisymmetrical component of the potential) corresponding to this new value. Hence, the radial velocity of the matter flow is given by the ratio of this angular momentum variation to the one fixed by the axisymmetric component of the potential and divided by the rotation period.
where . Using Eq. 2 this leads to
Acknowledgements.This research has made use of the GOLDMine Database. We thank the scientific and technical staffs at the IRAM 30 m and PdBI for their help in obtaining the observations presented here. We are grateful to the referee, Eric Emsellem, for helpful comments.
- Athanassoula (1992) Athanassoula, E. 1992, MNRAS, 259, 345
- Baker (2000) Baker, A. J. 2000, Ph.D. Thesis
- Barth & Shields (2000) Barth, A. J. & Shields, J. C. 2000, PASP, 112, 753
- Boone et al. (2006) Boone, F., Schilke, P., Muders, D., et al. 2006, in EAS Publications Series, ed. P. Stee, 299–305
- Boselli et al. (2006) Boselli, A., Boissier, S., Cortese, L., et al. 2006, ApJ, 651, 811
- Bottinelli et al. (1983) Bottinelli, L., Gouguenheim, L., Paturel, G., & de Vaucouleurs, G. 1983, A&A, 118, 4
- Buta & Combes (1996) Buta, R. & Combes, F. 1996, Fundamentals of Cosmic Physics, 17, 95
- Byrd et al. (1994) Byrd, G., Rautiainen, P., Salo, H., Buta, R., & Crocher, D. A. 1994, AJ, 108, 476
- Cayatte et al. (1990) Cayatte, V., van Gorkom, J. H., Balkowski, C., & Kotanyi, C. 1990, AJ, 100, 604
- Chyży et al. (2006) Chyży, K. T., Soida, M., Bomans, D. J., et al. 2006, A&A, 447, 465
- Combes (1996) Combes, F. 1996, in ASP Conf. Ser. 91: IAU Colloq. 157: Barred Galaxies, ed. R. Buta, D. A. Crocker, & B. G. Elmegreen, 286–+
- Combes (2004) Combes, F. 2004, in IAU Symposium, ed. T. Storchi-Bergmann, L. C. Ho, & H. R. Schmitt, 383–388
- Combes et al. (2004) Combes, F., García-Burillo, S., Boone, F., et al. 2004, A&A, 414, 857
- Combes & Gerin (1985) Combes, F. & Gerin, M. 1985, A&A, 150, 327
- Contopoulos & Grosbol (1989) Contopoulos, G. & Grosbol, P. 1989, A&A Rev., 1, 261
- Friedli & Benz (1993) Friedli, D. & Benz, W. 1993, A&A, 268, 65
- Gabel & Bruhweiler (2002) Gabel, J. R. & Bruhweiler, F. C. 2002, AJ, 124, 737
- García-Burillo et al. (2003a) García-Burillo, S., Combes, F., Eckart, A., et al. 2003a, in ASP Conf. Ser. 290: Active Galactic Nuclei: From Central Engine to Host Galaxy, 423–+
- García-Burillo et al. (2003b) García-Burillo, S., Combes, F., Hunt, L. K., et al. 2003b, A&A, 407, 485
- García-Burillo et al. (2005) García-Burillo, S., Combes, F., Schinnerer, E., Boone, F., & Hunt, L. K. 2005, A&A, 441, 1011
- Gavazzi et al. (2003) Gavazzi, G., Boselli, A., Donati, A., Franzetti, P., & Scodeggio, M. 2003, A&A, 400, 451
- Guhathakurta et al. (1988) Guhathakurta, P., van Gorkom, J. H., Kotanyi, C. G., & Balkowski, C. 1988, AJ, 96, 851
- Guilloteau & Lucas (2000) Guilloteau, S. & Lucas, R. 2000, in ASP Conf. Ser. 217: Imaging at Radio through Submillimeter Wavelengths, ed. J. G. Mangum & S. J. E. Radford, 299–+
- Helfer et al. (2003) Helfer, T. T., Thornley, M. D., Regan, M. W., et al. 2003, ApJS, 145, 259
- Hensler et al. (2003) Hensler, G., Bomans, D., & Boselli, A. 2003, The Cosmic Cauldron, 25th meeting of the IAU, Joint Discussion 10, 18 July 2003, Sydney, Australia, 10
- Ho (2003) Ho, L. C. 2003, in ASP Conf. Ser. 290: Active Galactic Nuclei: From Central Engine to Host Galaxy, ed. S. Collin, F. Combes, & I. Shlosman, 379–+
- Ho et al. (2001) Ho, L. C., Feigelson, E. D., Townsley, L. K., et al. 2001, ApJ, 549, L51
- Ho et al. (1997) Ho, L. C., Filippenko, A. V., & Sargent, W. L. W. 1997, ApJS, 112, 315
- Hummel et al. (1987) Hummel, E., van der Hulst, J. M., Keel, W. C., & Kennicutt, Jr., R. C. 1987, A&AS, 70, 517
- Jogee et al. (2005) Jogee, S., Scoville, N., & Kenney, J. D. P. 2005, ApJ, 630, 837
- Keel (1996) Keel, W. C. 1996, PASP, 108, 917
- Krips et al. (2005) Krips, M., Eckart, A., Neri, R., et al. 2005, A&A, 442, 479
- Laine et al. (1999) Laine, S., Kenney, J. D. P., Yun, M. S., & Gottesman, S. T. 1999, ApJ, 511, 709
- Laurikainen & Salo (2002) Laurikainen, E. & Salo, H. 2002, MNRAS, 337, 1118
- Lindblad & Lindblad (1994) Lindblad, P. O. & Lindblad, P. A. B. 1994, in ASP Conf. Ser. 66: Physics of the Gaseous and Stellar Disks of the Galaxy, ed. I. R. King, 29–+
- Maoz et al. (1998) Maoz, D., Koratkar, A., Shields, J. C., et al. 1998, AJ, 116, 55
- Nakanishi et al. (2005) Nakanishi, H., Sofue, Y., & Koda, J. 2005, PASJ, 57, 905
- Neff & Hutchings (1992) Neff, S. G. & Hutchings, J. B. 1992, AJ, 103, 1746
- Piner et al. (1995) Piner, B. G., Stone, J. M., & Teuben, P. J. 1995, ApJ, 449, 508
- Pogge et al. (2000) Pogge, R. W., Maoz, D., Ho, L. C., & Eracleous, M. 2000, ApJ, 532, 323
- Quillen et al. (1995) Quillen, A. C., Frogel, J. A., Kenney, J. D. P., Pogge, R. W., & Depoy, D. L. 1995, ApJ, 441, 549
- Regan et al. (1999) Regan, M. W., Sheth, K., & Vogel, S. N. 1999, ApJ, 526, 97
- Regan et al. (2004) Regan, M. W., Thornley, M. D., Bendo, G. J., et al. 2004, ApJS, 154, 204
- Regan et al. (1997) Regan, M. W., Vogel, S. N., & Teuben, P. J. 1997, ApJ, 482, L143+
- Rogstad et al. (1974) Rogstad, D. H., Lockart, I. A., & Wright, M. C. H. 1974, ApJ, 193, 309
- Rubin et al. (1999) Rubin, V. C., Waterman, A. H., & Kenney, J. D. P. 1999, AJ, 118, 236
- Sakamoto et al. (1999) Sakamoto, K., Okumura, S. K., Ishizuki, S., & Scoville, N. Z. 1999, ApJS, 124, 403
- Sanders & Huntley (1975) Sanders, R. H. & Huntley, J. M. 1975, BAAS, 7, 513
- Schinnerer et al. (2006) Schinnerer, E., Böker, T., Emsellem, E., & Lisenfeld, U. 2006, ApJ, 649, 181
- Sellwood & Wilkinson (1993) Sellwood, J. A. & Wilkinson, A. 1993, Reports of Progress in Physics, 56, 173
- Solomon & Barrett (1991) Solomon, P. M. & Barrett, J. W. 1991, in IAU Symp. 146: Dynamics of Galaxies and Their Molecular Cloud Distributions, ed. F. Combes & F. Casoli, 235–+
- Terashima et al. (2000) Terashima, Y., Ho, L. C., & Ptak, A. F. 2000, ApJ, 539, 161
- Tschöke et al. (2001) Tschöke, D., Bomans, D. J., Hensler, G., & Junkes, N. 2001, A&A, 380, 40
- Tully (1988) Tully, R. B. 1988, Nearby galaxies catalog (Cambridge and New York, Cambridge University Press, 1988, 221 p.)
- Vollmer et al. (2004) Vollmer, B., Balkowski, C., Cayatte, V., van Driel, W., & Huchtmeier, W. 2004, A&A, 419, 35
- Wada (1994) Wada, K. 1994, PASJ, 46, 165
- Wada & Habe (1995) Wada, K. & Habe, A. 1995, MNRAS, 277, 433
- Warmels (1988) Warmels, R. H. 1988, A&AS, 72, 57
- Wong et al. (2004) Wong, T., Blitz, L., & Bosma, A. 2004, ApJ, 605, 183