Unveiling the inner morphology and gas kinematics of NGC 5135 with ALMA
The local Seyfert 2 galaxy NGC5135, thanks to its almost face-on appearance, a bulge overdensity of stars, the presence of a large-scale bar, an AGN and a Supernova Remnant, is an excellent target to investigate the dynamics of inflows, outflows, star formation and AGN feedback. Here we present a reconstruction of the gas morphology and kinematics in the inner regions of this galaxy, based on the analysis of Atacama Large Millimeter Array (ALMA) archival data. To our purpose, we combine the available 100 pc resolution ALMA 1.3 and 0.45 mm observations of dust continuum emission, the spectroscopic maps of two transitions of the CO molecule (tracer of molecular mass in star forming and nuclear regions), and of the CS molecule (tracer of the dense star forming regions) with the outcome of the SED decomposition. By applying the BAROLO software (3D-Based Analysis of Rotating Object via Line Observations), we have been able to fit the galaxy rotation curves reconstructing a 3D tilted-ring model of the disk. Most of the observed emitting features are described by our kinematic model. We also attempt an interpretation for the emission in few regions that the axisymmetric model fails to reproduce. The most relevant of these is a region at the northern edge of the inner bar, where multiple velocity components overlap, as a possible consequence of the expansion of a super-bubble.
keywords:galaxies: active – galaxies: ISM – galaxies: kinematics and dynamics – galaxies: nuclei – galaxies: Seyfert – submillimetre: galaxies.
Active galactic nuclei (AGN) are thought to play a major role in the formation and evolution of galaxies, providing mechanisms for feedback from the supermassive black hole (SMBH) to its host galaxy and the intergalactic medium (see, e.g., Fabian 2012; Somerville & Davé 2015 and references therein). In co-evolutionary scenarios (e.g., Hopkins et al. 2007; Lapi et al. 2014) star formation (SF) activity and SMBH properties are closely connected, both in high redshift quasars and in local Seyfert nuclei, which are fuelled by accretion of material onto the SMBH. In the absence of major merging events and companions, as it is the case for the vast majority of local galaxies, the mechanisms that link SF and accretion activity lie in the inner galactic regions (within 1 kpc from the black hole, BH; Downes & Solomon 1998; Bryant & Scoville 1999) and are thought to be responsible for the feeding of the BH and the halting of the SF through feedback mechanisms from the BH itself.
The picture is complicated by the fact that the processes that drive inflows of gas towards the center are not univocally identified and seem to depend strongly on the scale considered. At kiloparsec scale, stellar bars are the most robust and efficient instabilities, constituting preferential directions of infall of the material towards the centre, possibly triggered by galaxy interactions. At the edges of the bars, shock fronts and dense dusty star forming regions are not unusual as a result of the sudden change of direction of the spiralling infalling material (e.g., Saito et al. 2017). Similarly, it is not uncommon to identify dynamically decoupled bars on different angular scales with overlapping dynamical resonances (Combes et al. 2013 and reference therein). In fact, on smaller scales, secondary nuclear bars seem to proceed from the primary bars (e.g., Shlosman & Heller 2002), and on 10-100pc scales, nuclear disks become unstable. On even smaller scales, the clumpiness of the interstellar medium (e.g., dynamical friction and viscous torques due to turbulence) starts to play a major role (e.g., Combes 2012).
In recent years, outflows of either ionised or molecular gas from the SMBHs have been observed and considered responsible for halting the infall and the SF in the inner galactic regions (e.g., Feruglio et al. 2010; Cicone et al. 2014; Harrison et al. 2016). Fuelling of the SMBH and outflows seem to constitute the self-regulating combination of processes on the small scales responsible of the morphology and dynamics of the whole galaxy. The co-existence of such different processes, how the small scales dynamics influences the overall galaxy morphology, the timescale on which different processes happen, and if different evolutionary stages justify the different observed morphologies are still open questions.
These processes can now be investigated thanks to the unprecedented sensitivity and resolution of the Atacama Large Millimetre Array (ALMA), with the observation of CO molecule transitions in local Seyfert galaxies offering the possibility of directly witnessing the molecular gas fuelling towards the nucleus, thus allowing for the first time a thorough knowledge of the accretion process in the innermost regions of AGN (i.e., on pc scales from the nucleus; see, e.g., García-Burillo et al. 2016).
For a few nearby Seyfert galaxies, sub-mm/mm continuum and spectroscopic observations are starting to become available in the ALMA archive. Recently, Gruppioni et al. (2016) have collected the multi-wavelength photometric data for a subsample of 76 Seyfert galaxies from the IRAS 12- galaxy sample (12MGS, Rush et al. 1993) selected for having available Spitzer-IRS and literature mid- and far-IR spectroscopic data. The sample, though incomplete for statistical studies, being randomly selected from a nearly complete sample, maintains a robust statistical significance and can be considered as representative of the local Seyfert galaxy population (as confirmed by a two-sided Kolmogorov-Smirnov test; see Gruppioni et al. 2016). The photometric data available in the literature – from the radio to X-ray bands – allowed for this sample a detailed spectral energy distribution (SED) decomposition (see an example in Figure 2), that, including the contribution of stars, re-emission from dust in star forming regions and from a AGN dusty torus, provided precise measurement of the main properties of AGN and its host galaxy, like, e.g., the AGN bolometric luminosity, the star formation rate (SFR), the total IR luminosity, and the stellar mass.
For one of these Seyfert galaxies, NGC 5135, ALMA data were available in the archive in two bands (both continuum and CO lines): band 6 (B6, 1.3 mm) and band 9 (B9, 500 m). In this paper we present the ALMA data, and discuss possible interpretation of the resulting complex morphology and dynamical structure of this galaxy. To this purpose, we considered the available 100 pc resolution ALMA B6 and B9 observations of dust continuum emission, the spectroscopic maps of two CO transitions (tracer of molecular mass in star forming and nuclear regions), and of a CS transition (tracer of the dense star forming regions) and the outcome of the SED decomposition. NGC5135, thanks to its almost face-on appearance, the presence of a large scale bar, the bulge overdensity of stars, the past evidences of the presence of a Supernova Remnant (SNR) and a central AGN, is an excellent target to investigate the dynamics of inflows, outflows, SF and AGN feedback.
In we present the information on NGC 5135 available from the literature, and the SED decomposition and IR spectroscopy results. In we present the ALMA archival data for this object and the data reduction/calibration that we have performed. In we discuss the complex gas morphology and kinematics in the inner regions of NGC 5135 and the possible responsible mechanisms. In we summarise our results.
Throughout this paper, we use a Chabrier (2003) initial mass function (IMF) and we adopt a CDM cosmology with = 71 km s Mpc, = 0.27, and .
NGC 5135, a nearby galaxy at = (implying a scale of 281 pc/arcsec), was first selected by Phillips et al. (1983) as part of a sample of nearby high-excitation emission-line galaxies, from the original catalogue of 719 bright galaxies by Sandage (1978). It was then classified as a Seyfert 2 by Véron-Cetty & Véron (2006) based on the [O III]5007/H and [N II]6584/H line ratios (BPT diagram: Baldwin et al. 1981). From its total IR luminosity (i.e., [8–1000 m] integrated), log(/ L)11.17, NGC 5135 was then classified as a Luminous Infrared Galaxy (LIRG) by Sanders et al. (2003). Its DSS image is shown in Figure 1, with overplotted the ALMA B9 continuum contours: the ALMA map is zoomed in the top right corner of the figure, while the zoom at the bottom right corner shows the same ALMA B9 contours on top of the HST (F606W) image (González Delgado et al. 1998), while the broad-band SED and the result of the decomposition performed by Gruppioni et al. (2016) are shown in Figure 2. From the broad-band (from UV to mm) SED fitting and decomposition, Gruppioni et al. (2016) have estimated the main physical parameters of NGC 5135, such as the stellar mass, 5.1610 M, the total IR luminosity, 10 L, the star formation rate, (15.61.9) M/yr, and the AGN bolometric luminosity, 10 erg/s. The values derived from the SED fitting correlate well with the luminosities of the mid-IR lines tracing SF (e.g., PAHs, [Ne II], [Si III]) and AGN activity (e.g., [Ne V], [O IV]), measured by Tommasin et al. (2008); Tommasin et al. (2010) from high resolution Spitzer-IRS spectra.
Nuclear studies of NGC 5135 showed a complex morphology and the presence of numerous coexisting structures in the region within 0.7 kpc from the AGN, including a nuclear bar (e.g., Mulchaey & Regan 1997). In order to describe the central engine, Bedregal et al. (2009) observed NGC 5135 with VLT SINFONI in [Si VI] line, and H-band () and K-band () continuum. The [Si VI] emission, that requires an ionisation energy of 166.7eV, is considered a good tracer of the AGN activity in the areas surrounding the nucleus, and it is mainly produced by gas excited just outside the broad line regions (Bedregal et al. 2011).
Our target was also observed by Levenson et al. (2004) in the X-rays with Chandra (0.4-8 keV), and the AGN emission component was separated from the stellar one. Two different components were identified in the central X-ray sources: a northern component associated to the AGN, and a southern associated to a star forming region (hereafter called SFr). The SFr found in the central region seems to contribute to the high obscuration of the AGN, producing an hydrogen column density . The SF regions were previously observed also by HST (González Delgado et al. 1998; Alonso-Herrero et al. 2006) and by VLT SINFONI (in Br m; Bedregal et al. 2009). Finally, a Supernova remnant (SNR) and a massive outflow (likely supernova driven) located in the Southern part of the nuclear region, were found by using VLT SINFONI observations of [Fe II] as a tracer (Colina et al. 2012).
At radio wavelengths, Ulvestad & Wilson (1989) found an intense emission at 6 and 20 cm in VLA-observations, with total fluxes of 163.2 mJy and 58.8 mJy respectively. The complexity of the radio emission, that presented a diffuse faint North emission with a bright source in the Southern region in both bands, allowed the authors to classify NGC 5135 as an ”ambiguous” radio source.
In order to complete this scenario, we have analysed ALMA archival data, adding the cold dust and molecular gas view to the great wealth of available information about NGC 5135.
3 The ALMA (archival) data
ALMA B6 observations were obtained in Cycle 3 (project 2013.1.00243.S, PI: Colina) in two different configurations of the 12-m antennas array. The “compact” configuration was obtained with 38 antennas distributed over the baseline range between 15.0 and 348.5 m, while the “extended” one with 44 antennas over the baseline range 15 – 1600 m, reaching resolutions of 0.16 arcsec at the sky frequency. The combination of the two configurations allows to cover angular scales in the range 0.1 – 12.32 arcsec. These scales correspond to spatial scales of 28.9 pc – 3.56 kpc at the redshift of NGC 5135. Two 1875 MHz spectral windows, centred respectively at the sky frequencies of 227.4 GHz, and 241.6 GHz, cover the expected frequency of CO(2–1) and CS(5–4) molecular line emissions with a channel width of 488 kHz (corresponding to 0.42 km at the observing frequency). A third 2-GHz wide spectral window centred at 228.8 GHz, with a channel width of 15.625 MHz, was dedicated to improve the continuum characterisation. The data available in the archive were calibrated using the ALMA calibration pipeline using CASA version 4.2.2. The calibrated data-sets of the two configurations were concatenated and combined with CASA 4.7.1. The achieved continuum beam size of the concatenated image is 0.25 0.20, while the 1 RMS is 1.9 Jy beam. The maximum recoverable scale allowed by the array configuration is 12 arcsec. Finally, the CASA standard cleaning was performed on the combined image, considering briggs - weighting.
ALMA B9 observations (cycle 3; project 2013.1.00524.S, PI: N. Lu) were taken on June 2 2015 with an array composed by 39 12-m antennas included in a range of baselines between 20.5 and 886 m. Four spectral windows were centred on the sky frequencies 678.2GHz, 680.2GHz, 683.7GHz and 682.0GHz respectively. At these frequencies, the angular scales observed are in the range 0.06 – 2.68 arcsec, corresponding to 17.33 pc – 774.38 pc at the source redshift. The expected sky frequency of CO(6–5) transition falls in the fourth spectral window.
Manual calibration was performed for the archival data. The phase calibrator () is faint at these observing frequencies ( 0.4 Jy). With respect to the script available in the archive, we improved the calibration by flagging four mis-behaving antennas and updated the flux density model for the flux density calibrator by a factor 0.84 (to the value of 1.39 Jy), exploiting an observation of the calibrator coeval to those of the target (not yet available in the calibration catalogue at the time of the quality assessment analysis). The beam of the B9 continuum image is 0.30v0.29, while the reached 1 RMS is 2.11 Jy beam. The maximum recoverable scale for the B9 configuration is 2.68 arcsec.
We must stress that, given the smaller size of the maximum recoverable scale in our B9 configuration, any direct comparison between the fluxes measured in the two bands will not be possible, since B9 might miss (i.e., filter out) diffuse fluxes – if any – on scales longer than the maximum recoverable ones.
The continuum images resulting from our ALMA data analysis are shown in the top panels of Figure 3 (B6 and B9 ). The 0 moment (i.e., intensity map in Jy km s beam) of the CO(2–1) () and CO(6-5) () are shown in the panels of the same figure. The 1 RMS obtained for the moment 0 maps is 1.1 and 1.5 Jy km s beam, respectively for CO(2–1) and CO(6–5). In Figure 4 we show the velocity field () and the velocity dispersion () for the CO(2–1) and (6–5) transitions (on the and respectively), while in Figure 5 we show the intensity map of CS (5–4), for which we have obtained an 1 RMS of 1.5 Jy km s beam. These images show very clearly the complex gas and dust distribution in this source, with hints of a central bar structure connecting two potential spiral arms.
Several regions deserving particular attention and interpretation are identified by combining the information derived from the continuum and the CO moment maps. In order to select the boundaries for these regions (labelled with letters and if needed encircled by ellipses in Figures 3 and 4), and measure their fluxes, we have used the CASA-viewer. We have considered the lower resolution band (B9 in our case, due to the antennae configuration), then we have selected the extraction region as to include all the flux in that band in the continuum image, applying the same extraction region to both bands. Note that, although the flux across the three SF regions (“C”, “D”, “E”) extends for about 3 arcsec in declination, we have tried to measure the flux connected to the single knots separately by fitting the three regions singularly. The observed CO(2–1) spectra for the most relevant emitting regions, are shown in Figure 6, and are described in detail in the following paragraph.
Region “A” (centred at RA(J2000) DEC(J2000) , as found from the B6 continuum image; see Section 4) includes the AGN position, as identified in VLT-SINFONI [Si IV] and Chandra X-ray maps (RA(J2000) DEC(J2000) ; Bedregal et al. 2009; Levenson et al. 2004). It is clearly visible in both the continuum maps, implying the presence of cold and dense dust in the very central regions, but it is not clearly detected in the observed spectral line transitions. Similarly, Schinnerer et al. (2000) did not detect the nucleus of NGC 1068 in CO, although very recently the NGC 1068 torus and multiple knots of SF have been resolved in CO with ALMA at much higher resolution (4.2 pc; García-Burillo et al. 2016). Similarly to what happens in NGC 1068, where the AGN is part of a circum-nuclear disk that could be resolved only with very high resolution observation (García-Burillo et al. 2016), the flux from region “A” can be contaminated by SF and only partially due to the AGN dusty torus. However, at odds with the dust distribution in NGC 1068, the continuum emission from region “A” in NGC 5135 appears unresolved at the resolution of the B6 data, but already isolated (i.e., not embedded in a circum-nuclear ring, which is expected to extend out to larger distance). We stress that the circum-nuclear disk observed by García-Burillo et al. (2014) in NGC 1068 at a resolution of 60–70 pc was successively resolved into both AGN torus and SF knots, once observed at much higher resolution (4.2 pc; García-Burillo et al. 2016). Therefore, given fact that the B9 emission in region “A” is significantly larger than the beam (see the 4 contours in Figure 3), and given the resolution of these data, we are not able to separate out SF knots (if any) and AGN. In any case, since it is not possible to estimate the continuum emission from the AGN torus without SF contribution unless the AGN torus is clearly resolved, we can consider the flux measured from region “A” as an upper limit to the cold dust emission component from the torus. The ALMA B6 (1.3 mm) and B9 (450 m) continuum measurements corresponding to the central region (e.g., “A”) are shown as green stars in Figure 2. Our ALMA measurements seem to show much higher cold dust contribution from the nucleus than expected from torus only models. ALMA higher resolution measurements would be necessary for this galaxy, in order to be sure to properly model and interpret the emission from the AGN (i.e., without any SF contamination), so far totally unconstrained by data in the Rayleigh Jeans tail. Only by resolving the torus emission itself we can derive its real contribution, although if the flux measured in the central region will come out to be mostly due to the AGN, this result will challenge the current torus models, implying the presence of a colder dust component. However, as previously noticed, the B9 continuum measurement can just be taken as an indication and is not directly comparable with the B6 one, since the smaller recoverable scales in the B9 configuration filter out signals at larger scales, if any.
Regions “B” and “C” correspond to the bar edges. The sudden change of the in-falling direction in this region seems to produce a dense and hot environment where SF is active. The CO(2–1) spectrum of region “B” peaks at the galaxy velocity (i.e., 0 km s in Figure 6), but shows a positive velocity tail (possibly a secondary peak at about 40-50 km s), while the spectrum of the “CD” regions peaks at negative velocities (e.g., 30/40 km s). Regions “C”, “D” and “E” had been already identified as SF loci by Bedregal et al. 2009: this is confirmed now by the presence of strong CO emission. Region “E” is also the dominant region in the CS(5–4) brightness distribution (flux density 0.290.03 Jy/beam, see fig. 5), thus likely to be denser than the others. A detailed study of the main properties, physical parameters and SF rates in these regions will be presented and discussed in a further paper from our team, where a detailed chemical analysis will also be performed.
Region “G” corresponds to the SNR detected by HST and VLT SINFONI observations, and clearly visible here in all the observed emissions. Observation with VLA at 6 cm (Ulvestad & Wilson 1989) also shows a strong emission at this location, suggesting strong synchrotron emission from SNR.
Two regions of high velocity dispersion, namely regions “F” and “H”, were not detected in any previous observations:
– Region “F” is aligned along the bar, but seems to lie outside its edge. It is characterised by a double peaked spectrum (see Figure 6) with two opposite velocity components (70 km s). It is not detected in continuum emission.
– Region “H” lies along the southern arm and appears only in the CO(2–1) maps, not in continuum or CO(6–5), showing a high velocity dispersion and a secondary peak in the spectrum corresponding to large negative velocity ( 70/80 km s; see Figure 6).
4 Barolo kinematic modelling
In order to investigate the kinematics n the innermost regions of NGC 5135, we used the BAROLO software to fit our ALMA data cubes. BAROLO (3D-Based Analysis of Rotating Object via Line Observations) is a code developed by Di Teodoro & Fraternali (2015) which derives rotation curves of galaxies from emission line data cubes by fitting a 3D tilted-ring model. The tilted-ring model is based on two assumptions: the material which is responsible for the emission lines is contained into a thin disc, and the kinematics is dominated by rotation. In this model, the galaxy disc is divided into a series of concentric rings characterised by seven free parameters (three of which describe the morphology and four the kinematics): the kinematical center (), the inclination with respect to the observer () and the position angle of the major axis () define the characteristics of the rings projected on the sky plane; the velocity dispersion (), the rotational velocity (), the systemic velocity () and the radial velocity () of the gas, define the kinematic characteristics of the model. In the BAROLO version used in this work (v1.3), the three velocity components are combined to define the observed velocity along the line of sight (), mapped by moment 1. The derived velocity map can be written in terms of a generic harmonic expansion (e.g., Schoenmakers et al. 1997; Swaters et al. 1999) that we can write as (Roelofs 1995 for details):
where is the position angle of the major axis on the receding half of the galaxy, taken anticlockwise from the North direction on the sky. The code finds the best-fitting values by minimising the difference between the data-cube and the kinematic model smoothed to the same resolution of the data. BAROLO provides three methods for calculating the residuals by comparing the model and the data for every pixel in the 3D space of the data-cube: in our case, we have chosen a kind of , described as follows. If M is the flux associated with a model-pixel and D is the flux of a data-pixel (the flux of the RMS of the data-cube is associated to D in case it is a blank pixel; see Di Teodoro & Fraternali 2015 for details), then BAROLO minimises the quantity:
A further parameter to be specified is the column density of the gas. For this parameter, two different normalisations are implemented in BAROLO: pixel-by-pixel (i.e., local) and azimuthally averaged (i.e., azim). In the former case, the model assumes equal value for the integral of each spatial pixel along the spectral dimension in the model and in the observations. This normalisation allows us to consider non-axisymmetric density models and avoids that regions with anomalous gas distribution, like e.g. clumpy emission or holes, affect the global fit. The normalisation value of the second model is instead set to the azimuthally averaged flux in each ring. Although the azim normalisation cannot reproduce the details of the galaxy morphology (which are instead reproduced on a pixel-by-pixel basis by the local normalisation), it provides useful information about some of the model parameters (e.g., the inclination angle ) and is important to validate the local model solutions. In fact, since the observed gas distribution is non axisymmetric, it is not obvious that the solutions found by the two different normalisation models coincide, in case the result is not the correct one: to confirm that the kinematic model correctly reproduces the observed global motion of the gas, the local model solutions must be compatible with those found by the azim one.
Some of the model parameters have been fixed in order to reduce the computation-time and the degeneracy. For some of them the values have been fixed based on assumptions about the gas distribution (e.g., we have fixed the thickness of the disk by assuming a thin disk distribution), while other parameters were measured from the data, then fixed after checking their validity. For instance, the kinematic center - derived by setting a circular region around the AGN and assuming the emission peak of the B6 continuum coincident with the position of the AGN - was initially let free to vary. Once verified that BAROLO was finding the same kinematic position (within errors), we have fixed its value. Other two parameters were fixed in our model, in addition to the galaxy centre: the position angle () and the inclination of the galaxy disc (). By using the code kpvslice, an extension of the Karma software (Gooch 1996), we found deg as the best approximation of the major (projected) axis orientation. The inclination () of the disc is estimated by using the DSS image of NGC 5135 (see Figure 7) and matching the isophotal contours of the outer disc with an ellipsis. The inclination of the major axis of the outer ellipsis defines the inclination of the disk. This method provided an estimate of the inclination angle = deg, with the associated error derived by considering the average difference between the two red ellipses in Figure 7.
4.1 Results: Kinematic models for molecular gas of NGC 5135
In Figure 8 we show the results of our fit to the ALMA CO(2–1) data-cube (up to an angular size of 6 arcsec, corresponding to 1.7 kpc from the galaxy centre) with BAROLO, in seven representative velocity channels (out of 52) of spectral resolution 10 km s. The velocity of the galaxy is centred on the CO(2–1) line emission peak at 227.43 GHz, derived from the continuum B6-concatenated observation at the central AGN position. In the row of the figure we show the data velocity field, while in the and rows we plot the results of the local and azim model respectively, in the same velocity channels. Again, we stress how the comparison between the azim and the local model is a probe of the goodness of the local fit, since the azim model reproduces the “global” behaviour of the data and contains the same amount of gas as the local model. The region showing two opposite velocities (i.e., 70 km s, region “F”), cannot be reproduced by the models, while the other high velocity dispersion region located in the southern spiral arm (region “H”) can be only partially reproduced, with a tail at 40 km s not fitted by the model.
The comparison of our molecular gas distribution with the BAROLO model expectations gives us several important information about the gas structure and kinematics in the inner regions of NGC 5135. In Figure 9 we show the three moments (0, 1 and 2) of the CO(2–1) molecular gas distribution in the centre of NGC 5135 from the data cube ( row), and from the local ( row) and azim ( row) models obtained with BAROLO. Note that we have limited our moment analysis to the central regions of the galaxy, within a disc of radius 4 (1.167 kpc), because, as the ring radius increases, the CO(2–1) distribution strongly departs from axisymmetry, thus BAROLO cannot find acceptable solutions. Both models reproduce reasonably well the observed velocity field (moment 1) of the CO(2–1), shown in panel (b) of Figure 9 (models in panels (e), (h)), while only the local model provides a good fit to the observed gas distribution (see panel (d) versus (a)). In fact, since the azim model averages the brightness distribution ring-by-ring, it is not able (by definition) to reproduce the pixel-by-pixel asymmetry of the data (see panel (g) versus (a)). The complex structure of the velocity dispersion map (expressed as line width) observed in the data, showing zones of very large velocity dispersion, not expected for a rotation disk with super-imposed radial motions, is not reproduced by any of the models. In particular, in panel (c) we clearly note the two zones of very large velocity dispersion already identified in Section 3: one at the edge of the potential bar, showing two opposite velocities corresponding to the same spatial position (i.e., the bright yellow spot in panel (c) of Figure 9, corresponding to region “F”), and the other in the southern arm of the spiral structure (i.e., the orange region in the figure, corresponding to region “H”). Being the model axisymmetric, it cannot reproduce these local regions of high velocity dispersion in the data, but finds an averaged solution on the entire ring. This happens with both the azim and the local model. We note however that BAROLO is nonetheless able to identify a gradient, in fact the larger values of the velocity dispersion in outer rings of panels (f), (i), are due to the presence of region “F” and are not intrinsic dispersions: without the high dispersion regions, the model would reproduce reasonably well the average observed values within each ring. Finally, we note that the presence of the nuclear bar at the centre of NGC 5135 introduces motions in the innermost gas that the BAROLO software is unable to reconstruct perfectly. Therefore, we cannot consider our best-fit model as fully explanatory of the gas kinematic in NGC 5135, but only an approximation useful to understand the overall gas motion.
In Figure 10 we show the position-velocity (P-V) diagrams along the kinematic major () and minor () axes of the galaxy disk of NGC 5135 obtained with BAROLO. In the presence of a purely rotating disk, the P-V along the minor axis is expected to be flat, while it assumes an “S” shape in presence of radial motions (e.g., Fraternali et al. 2002). Although the P-V diagram along the major axis is similar to what we would expect for a regular spiral galaxy described by a rotating-disk model, the gas distribution along the minor axis shows three peculiar regions that cannot be described by the same model, even including a radial velocity component. Note that in NGC 5135 we are not only observing the P-V pattern expected for radial motions, but also velocities along the minor axis as large as those along the major axis. This is an unusual behaviour in regular spiral galaxies. Moreover, although the overall galaxy structure is well explained by the rotation disk plus radial motions model (see P-V along the major axis), there are clearly three velocity regions along the minor axis that either of the two models (i.e., local and azim) are unable to explain. This is visible at an offset of 1.5 arcsec, as a zone totally empty of data that the model would instead predict to be filled. The two “blobs” above and below the empty region (at 1.5 arcsec offset, 70 and 70 km s, corresponding to the two peaks in the left panel of Figure 10), are likely related to the presence of an expanding “superbubble” (see next Section).
In Figure 11, we show the rotation curve (; panel), the velocity dispersion (; panel) and the radial velocity (; panel) reconstructed by BAROLO as a function of the distance from the galaxy center (up to 1 kpc).
The stability of our results is confirmed by the good agreement between the solutions of the two BAROLO models, which are perfectly consistent within the error-bars.
In Figure 12, the same parameters of Figure 11 are plotted (as derived from the local model)
for different values of the position angle (i.e., varying the PA by steps of 5 within the error, around the best-fit value, i.e., 10):
this shows how robust our results are with respect to changing the PA value. Although the major axis is extracted along the PA, and the minor axis along PA90,
the P-V diagrams, as well as the derived velocities, are almost unaffected by these changes.
From the resulting rotation curve (), we observe an increasing internal rotation of the molecular gas, passing from 20 km s to 50 km s
at 1 kpc from the centre, likely to increase further in the outer parts of the galaxy. This is not an obvious results, since it shows a significant (increasing) rotation of the molecular gas
also in the very inner parts of the galaxy, expected to be “flat” and lower in typical nearby spiral galaxies (e.g., Frank et al. 2016).
Note that the presence of the bar in NGC 5135 could cause an underestimation of the rotation velocity: since the gas is streaming along the bar, observed velocities are close to rigid-body motion of the bar potential, and hence slower than circular velocity (e.g., Sofue 2016).
The velocity dispersion ( panel) shows a modest increase with increasing distance, passing from 20 to 30 km s, while the radial velocity () keeps almost constant, at around 30 km s, from the galaxy center to 1 kpc. We stress again that the presence of the bar can influence also the derived , whose significantly high values can be due to large non-circular motions of the gas inside the bar.
5 The possible origins of the complex gas kinematic in NGC 5135
As discussed in the previous section, the overall structure of the molecular gas CO(2–1) in NGC 5135, its density and kinematics can be well explained by a model considering a rotating
disk perturbed by radial motions. The direction of the radial flow of gas (i.e., in- or out-flow) can be estimated if the closer spiral arm is known, based on the inclination of the disk.
Assuming that the arms are trailing (i.e., the southern arm is receeding), we find that the molecular gas should flow towards the center of the galaxy: we are facing a gas inflow
towards the BH (for this reason the radial velocity in Figure 11 is negative).
However, the gas structure is not properly regular, as there are several regions that fall outside the rotating disk gas inflow model explanation.
In Figure 3 we have identified these regions, plus others that would deserve further explanation (marking them with letters and - if needed - encircling with ellipses), for a detailed discussion on the origin and nature of each of them.
Previous observations with different instruments (e.g., HST, VLT-SINFONI; e.g., González Delgado et al. 1998; Alonso-Herrero et al. 2006; Bedregal et al. 2009; Colina et al. 2012) have identified the central AGN (marked with the “A” letter), few SF regions at the left edge of the inner bar (“C”, “D”, “E”), and a SNR in the Southern arm of the inner spiral gas structure (“G”). The AGN region (“A”) is observable in continuum (both B6 and B9), but not in CO line, in agreement with the finding that the BH is strongly obscured (i.e., Compton-Thick, CT, from X-ray data; Levenson et al. 2004): the cold toroidal dust distribution has the size of the beam in ALMA B6, thus being unresolved at a linear scale of 56 pc. The SF regions (“C”, “D”, “E”) are clearly identifiable in the data in continuum, CO (both (2–1) and (6–5)) and CS emission. Regions “B” and “F” were not previously identified in any other bands, and seem to trace a shock front at the right edge of the bar, with an expanding super-bubble (see below) and gas likely thrown out of the galaxy at high velocity (in two opposite directions, thus showing a very large, double peaked width of the CO line). The alignment between zones “A”, “B” and “C” seems to indicate the presence of a bar, with enhanced SF at its edges. Zone “B” is likely to trace a shock front between the AGN and zone “F” (see further explanation), with the latter identifying the two opposite velocities region (clearly visible also in Figure 10 along the minor axis at 1.5 arcsec and 70 km s) at the edge of the bar edge. region “F” lies just outside the bar, close to the shock front and is the more puzzling region observed in these data (a tentative explanation is given in the next paragraph).
Both zones (“F” and “H”) were not identified by any previous observations in other bands.
Since we have chosen to be conservative and limit our moment analysis to the inner 4 arcsec (1.167 kpc from the center), where the BAROLO software provides more secure results, region “H” is out of the area considered for the fit. Indeed, the channel-by-channel plot (i.e., Figure 8, extending to 6 arcsec, 1.7 kpc) shows that the CO(2–1) distribution in this region can be only partially reproduced, with a tail at 40 km s remaining unfitted. Our interpretation of region “H” is that here there could likely be an overlap of two gas components: one following the southern spiral arm motion, the other flowing in opposite direction (e.g., the second peak at about 80 km s in the spectrum shown in Figure 6). The gas component following the arm flow is fitted by the model, the material flowing in the other direction is not.
Regarding region “F”, local emission from a SNR could be a simple explanation, but, differently from zone “G”, in “F” we do not observe significant emission either in continuum nor in spectral lines, but only a very large CO(2–1) velocity dispersion (i.e., produced by the presence of the two opposite velocity peaks, from two gas clouds spatially overlapping). Moreover, a SNR in this region cannot be confirmed by data in other bands, as in the case of region “G”. Hence the SNR origin hypothesis is likely to be discarded for no clear evidence of the presence of the remnant. The possibility of an outflow from the AGN not lying along the bar (i.e., perpendicular to the galaxy plane and just projected in zone “F”) could be considered, but in this case the origin of the shock front in “B” will be unrelated to the outflow and unclear. Moreover, an outflow of material due to the AGN is also not clearly justified by any other evidences, and it is very unlikely that, by chance, the projected position of the out-flowing material along a jet would fall exactly at the edge of the bar. Finally, the outflow interpretation cannot explain the material expelled in two opposite velocities in “F”. Our favourite interpretation is that in “F” we are facing a flow of material towards the AGN: an in-flow of gas along the bar (connecting the SF region “C”, the AGN “A” and the SF region “B”) could justify the presence of the shock front (“B”). The gas, in-falling through the inner arms by spiralling towards the central BH, could suddenly change its direction as approaching the nucleus, to start following the bar structure: such a flow is likely described by the continuum and CO(2–1) velocity pattern observed in the continuum and moment maps (see panels of Figures 3 and 4).
According to our interpretation, the observed high velocity dispersion region in zone “F” associated to the zone totally empty of gas between the two opposite velocity clouds (that the BAROLO model would instead predict to be filled) resembles an expanding super-bubble. Super-bubbles are known to be associated to very massive stars (OB-associations): strong stellar winds and subsequent SN explosions from those stars inject energy and mass into the ambient ISM, creating shock fronts that sweep-up the ISM (e.g., Castor et al. 1975). As the SN explosions start occurring within the cavity formed by the stellar wind bubbles, super-bubbles are created (e.g., McCray & Kafatos 1987) which may eventually reach kiloparsec sizes. These explosions never form a visible SNR, but instead expend their energy in the hot interior as sound waves. Both stellar winds and stellar explosions thus power the expansion of the super-bubble in the ISM. The interstellar gas swept up by super-bubbles generally cools, forming a dense shell around the cavity (observed in HI and H). The bubble interior contains hot (10 K), rarefied material, usually associated with extended diffuse X-ray emission (thus appearing as an empty cavity in cold and molecular gas). In this interpretation the two blobs of gas above and below the observed empty cavity in region “F” (at 1.5 arcsec offset, with velocities of 70 and 70 km s, corresponding to the two peaks in Figure 6), can be easily explained as gas surrounding the “super-bubble” being thrown away in different directions by the expanding spherical shell (see Kamphuis et al. 1991; Boomsma et al. 2008 for detected velocities up to 100 km s associated to “holes” in HI). In Figure 13 we show different CO(2–1) velocity channels, where the “super-bubble”, surrounded by material moving in opposite directions (e.g., peaks at -70 km s and 70 km s, not present at the systemic velocities), shows up as an empty area (encircled in the plot), with a clear shock front at its left hand side, peaking around the systemic velocities (0 km s).
In Table 1 we list the different zones identified by our data in the inner regions of NGC 5135 and summarise their main characteristics and our interpretation.
|A||central, unresolved at 56 pc, continuum only||AGN dusty torus, also studied by Bedregal et al. (2009)|
|B||high moment 0 CO(2–1), at the edge of the bar||shock front, enhanced SF|
|C||high continuum and CO mom 0 (both B6 and B9)||SF region, also studied by Bedregal et al. (2009)|
|D||high continuum and CO mom 0 (both B6 and B9)||SF region, also studied by Bedregal et al. (2009)|
|E||high continuum and CO mom 0 (both B6 and B9), but also CS(5–4)||SF region, likely denser than C, D|
|F||high mom 1 and 2 of CO(2–1), opposite velocity cloudsempty zone||expanding “super-bubble”, not observed before|
|G||high continuum and CO mom 0 (both B6 and B9)||SNR also studied by Colina et al. (2012)|
|H||high mom 1 and 2 of CO(2–1)||material flowing opposite to the arm, not observed before|
The complex kinematics observed in the central kpc of NGC 5135 is also intriguing when compared to the most recent models attempting to describe the dynamics and SF in the central molecular zone of the Milky Way. In particular, the models described by Kruijssen et al. (2014); Kruijssen et al. (2015); Krumholz & Burkhart (2016) predict that the gas is inflowing within the inner Lindblad resonance, where acoustic instabilities are stable to gravitational collapse and drive turbulence in the gas, maintaining large line-widths. As the gas streams in, the rotation curve passes from flat to solid body. In this region a substantial quantity of gas is accumulated, the shear drops and the turbulence is dissipated, making the gas prone to gravitational instability. In the nuclear regions of NGC 5135 we are likely facing radial motions, with the gas streaming from kpc scales to . At scales smaller than 200 pc, the BAROLO rotation curve is steep, then it flattens at 200R400 pc, to then steepen again towards the outer regions (at least up to 1 kpc). However, within 200 pc, both the azim and the local models of BAROLO find a minimum value for the CO(2–1) velocity dispersion (see Figure 11), which then increases at larger distances (due to the presence of the SF regions the SNR and the super-bubble, observed at the outer edge of this “ring”). This, in analogy with what happens in our Galaxy, might suggest that in NGC 5135 the gas starts accumulating and becoming gravitationally unstable at a scale of 200-300 pc from the nucleus.
Although extremely interesting and complex, the picture of the inner regions of NGC 5135 emerging from the current ALMA data is still incomplete: the dusty torus is unresolved in continuum and undetected in CO, and the CO velocity maps raise the issue of what generates the two (opposite) high velocity dispersion components at the nuclear bar edge and the possible “super-bubble”. This is a totally new result, needing further investigation. To resolve the cold dust (at 10–20 pc, possibly at clump level) of the toroidal distribution surrounding and fuelling the AGN at the centre of the nearby Seyfert NGC 5135, measuring its Rayleigh-Jeans tail and testing the current models, could be obtained only with higher resolution observations (e.g., 0.03–0.05 in continuum – 5–7 times better than current), both in continuum and line, allowing us to spatially resolve the torus and to trace its molecular structure with high density tracers. Such high resolution would also allow us to study in detail the complex kinematic of the high velocity dispersion region, in order to spatially resolve the zone with opposite velocities, investigating the nature of the material flowing in or out the nuclear region, and the possible presence of the super-bubble.
In order to investigate the mutual interplay between AGN and SF activity within galaxies, their relative role in outflow or inflow phenomena, that can be related to either negative feedback towards SF or AGN fuelling with enhanced SF, we have analised the ALMA archival data in B6 and B9 for the nearby Seyfert galaxy NGC 5135. Continuum and line maps have been obtained, showing a complex gas distribution and kinematics. We have modelled the data with a code fitting the line data cubes of rotating galaxies with a 3D tilted-rings model (BAROLO, Di Teodoro & Fraternali 2015), in order to obtain a clearer picture of the phenomena occurring in the inner regions of our target. Several possible explanations for the complex geometry and gas dynamics in this galaxy have been discussed, with the main result being the discovery of radial motions overimposed to the galaxy disk rotation, with clear evidence of them consisting of gas inflow towards the nucleus, occurring through inner spiral arms and a bar. Several regions of enhanced SF due to shocks, outflows and super-bubbles have been identified, clear signs of a turbulent activity near the galaxy centre and the AGN.
The main results of this work can be here summarised as follows:
the ALMA B6 and B9 continuum and CO line data show a complex morphology and gas kinematics in the inner regions of the nearby Seyfert galaxy NGC 5135;
at the locus of the central AGN (previously observed with VLT-SINFONI and in X-ray) we detect unresolved flux in both ALMA bands, at a galaxy scale of 59 pc in the B6 continuum, while do not detected any CO emission, in agreement with the AGN being heavily obscured (as from X-ray). The unresolved flux could contain a possible contribution from circum-nuclear SF: only higher resolution data will allow us to derive the AGN contribution only;
the fluxes measured at 0.45 and 1.3 mm in the central region, considered an upper limit to the AGN dusty torus emission due to the possible SF contamination, largely exceed the SED decomposition expectations for torus only (by about two orders of magnitude). If the flux is dominated by the AGN (with negligible contribution from circum-nuclear SF), this result will imply a significantly higher contribution from cold dust heated by the AGN than predicted by the torus models;
the inner regions of the galaxy observed with ALMA show a barred spiral-like structure with two arms and a bar, connecting the central AGN with two enhanced SF regions at the edges of the bar itself;
the analysis of the CO(2-1) moment maps (e.g., gas distribution, velocity field, velocity dispersion) reveals very complex structures and gas kinematics, part of which we tried to explain by rotating disk plus radial motions with turbulence at the edge of the bar;
the regions not explainable by the kinematic model identify potential gas outflows/inflows from/towards the AGN, associated to enhanced SF activity, a shock front and an expanding super-bubble.
Given its complex structure and the multiple components shown by the ALMA observations of the innermost regions of NGC 5135, this galaxy constitutes a crucial benchmark for understanding the feeding and the feedback mechanisms associated to SMBHs. New observations with higher resolution will be crucial to understand the physics and the geometry of the inner AGN region. In a forthcoming paper we will analise in detail the physical properties of the different structures identified in this work, such as the mass flow rate and the energy associated to the gas inflow and outflow, investigating the role of these regions and of the related processes in the feedback or the enhancing of SF, and in feeding the BH (e.g., the rate of gas flow towards the AGN).
The authors wish to thank F. Fraternali and E. Di Teodoro for their great help in the use of the BAROLO software and for very fruitful discussions about the interpretation of the results of this work. This paper makes use of the following ALMA data: 2013.1.00243.S and 2013.1.00581. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. CG acknowledges funding from the INAF PRIN-SKA 2017 program 1.05.01.88.04. MM acknowledges partial financial support by the Italian Ministero dell’Istruzione, Università e Ricerca through the grant Progetti Premiali 2012–iALMA (CUP C52I13000140001).
- Alonso-Herrero et al. (2006) Alonso-Herrero A., Rieke G. H., Rieke M. J., Colina L., Pérez-González P. G., Ryder S. D., 2006, ApJ, 650, 835
- Baldwin et al. (1981) Baldwin J. A., Phillips M. M., Terlevich R., 1981, PASP, 93, 5
- Bedregal et al. (2009) Bedregal A. G., Colina L., Alonso-Herrero A., Arribas S., 2009, ApJ, 698, 1852
- Bedregal et al. (2011) Bedregal A. G., Colina L., Azzollini R., Arribas S., Alonso-Herrero A., 2011, in Wang W., Lu J., Luo Z., Yang Z., Hua H., Chen Z., eds, Astronomical Society of the Pacific Conference Series Vol. 446, Galaxy Evolution: Infrared to Millimeter Wavelength Perspective. p. 83 (arXiv:1102.0012)
- Berta et al. (2013) Berta S., et al., 2013, A&A, 551, A100
- Boomsma et al. (2008) Boomsma R., Oosterloo T. A., Fraternali F., van der Hulst J. M., Sancisi R., 2008, A&A, 490, 555
- Bryant & Scoville (1999) Bryant P. M., Scoville N. Z., 1999, AJ, 117, 2632
- Castor et al. (1975) Castor J., McCray R., Weaver R., 1975, ApJ, 200, L107
- Cicone et al. (2014) Cicone C., et al., 2014, A&A, 562, A21
- Colina et al. (2012) Colina L., Pereira-Santaella M., Alonso-Herrero A., Bedregal A. G., Arribas S., 2012, ApJ, 749, 116
- Combes (2012) Combes F., 2012, in Journal of Physics Conference Series. p. 012041 (arXiv:1111.4770), doi:10.1088/1742-6596/372/1/012041
- Combes et al. (2013) Combes F., et al., 2013, A&A, 558, A124
- Di Teodoro & Fraternali (2015) Di Teodoro E. M., Fraternali F., 2015, MNRAS, 451, 3021
- Downes & Solomon (1998) Downes D., Solomon P. M., 1998, ApJ, 507, 615
- Fabian (2012) Fabian A. C., 2012, ARA&A, 50, 455
- Feruglio et al. (2010) Feruglio C., Maiolino R., Piconcelli E., Menci N., Aussel H., Lamastra A., Fiore F., 2010, A&A, 518, L155
- Frank et al. (2016) Frank B. S., de Blok W. J. G., Walter F., Leroy A., Carignan C., 2016, AJ, 151, 94
- Fraternali et al. (2002) Fraternali F., van Moorsel G., Sancisi R., Oosterloo T., 2002, AJ, 123, 3124
- García-Burillo et al. (2014) García-Burillo S., et al., 2014, A&A, 567, A125
- García-Burillo et al. (2016) García-Burillo S., et al., 2016, ApJ, 823, L12
- González Delgado et al. (1998) González Delgado R. M., Heckman T., Leitherer C., Meurer G., Krolik J., Wilson A. S., Kinney A., Koratkar A., 1998, ApJ, 505, 174
- Gooch (1996) Gooch R., 1996, in Jacoby G. H., Barnes J., eds, Astronomical Society of the Pacific Conference Series Vol. 101, Astronomical Data Analysis Software and Systems V. p. 80
- Gruppioni et al. (2016) Gruppioni C., et al., 2016, MNRAS, 458, 4297
- Harrison et al. (2016) Harrison C. M., et al., 2016, MNRAS, 456, 1195
- Hopkins et al. (2007) Hopkins P. F., Richards G. T., Hernquist L., 2007, ApJ, 654, 731
- Kamphuis et al. (1991) Kamphuis J., Sancisi R., van der Hulst T., 1991, A&A, 244, L29
- Kruijssen et al. (2014) Kruijssen J. M. D., Longmore S. N., Elmegreen B. G., Murray N., Bally J., Testi L., Kennicutt R. C., 2014, MNRAS, 440, 3370
- Kruijssen et al. (2015) Kruijssen J. M. D., Dale J. E., Longmore S. N., 2015, MNRAS, 447, 1059
- Krumholz & Burkhart (2016) Krumholz M. R., Burkhart B., 2016, MNRAS, 458, 1671
- Lapi et al. (2014) Lapi A., Raimundo S., Aversa R., Cai Z.-Y., Negrello M., Celotti A., De Zotti G., Danese L., 2014, ApJ, 782, 69
- Levenson et al. (2004) Levenson N. A., Weaver K. A., Heckman T. M., Awaki H., Terashima Y., 2004, ApJ, 602, 135
- McCray & Kafatos (1987) McCray R., Kafatos M., 1987, ApJ, 317, 190
- Mulchaey & Regan (1997) Mulchaey J. S., Regan M. W., 1997, ApJ, 482, L135
- Phillips et al. (1983) Phillips M. M., Charles P. A., Baldwin J. A., 1983, ApJ, 266, 485
- Roelofs (1995) Roelofs G., 1995, PhD thesis, THE UNIVERSITY OF CHICAGO.
- Rush et al. (1993) Rush B., Malkan M. A., Spinoglio L., 1993, ApJS, 89, 1
- Saito et al. (2017) Saito T., et al., 2017, ApJ, 835, 174
- Sandage (1978) Sandage A., 1978, AJ, 83, 904
- Sanders et al. (2003) Sanders D. B., Mazzarella J. M., Kim D.-C., Surace J. A., Soifer B. T., 2003, AJ, 126, 1607
- Schinnerer et al. (2000) Schinnerer E., Eckart A., Tacconi L. J., Genzel R., Downes D., 2000, ApJ, 533, 850
- Schoenmakers et al. (1997) Schoenmakers R. H. M., Franx M., de Zeeuw P. T., 1997, MNRAS, 292, 349
- Shlosman & Heller (2002) Shlosman I., Heller C. H., 2002, ApJ, 565, 921
- Sofue (2016) Sofue Y., 2016, PASJ, 68, 2
- Somerville & Davé (2015) Somerville R. S., Davé R., 2015, ARA&A, 53, 51
- Swaters et al. (1999) Swaters R. A., Schoenmakers R. H. M., Sancisi R., van Albada T. S., 1999, MNRAS, 304, 330
- Tommasin et al. (2008) Tommasin S., Spinoglio L., Malkan M. A., Smith H., González-Alfonso E., Charmandaris V., 2008, ApJ, 676, 836
- Tommasin et al. (2010) Tommasin S., Spinoglio L., Malkan M. A., Fazio G., 2010, ApJ, 709, 1257
- Ulvestad & Wilson (1989) Ulvestad J. S., Wilson A. S., 1989, ApJ, 343, 659
- Véron-Cetty & Véron (2006) Véron-Cetty M.-P., Véron P., 2006, A&A, 455, 773