TWO-DIMENSIONAL MAPPING OF YOUNG STARS IN THE INNER 180 pc OF NGC 1068: CORRELATION WITH MOLECULAR GAS RING AND STELLAR KINEMATICS

NGC 1068 is the closest, brightest, and for this reason, the most studied Seyfert 2 galaxy. Its distance of only 14.4 Mpc (Crenshaw & Kraemer 2000) corresponds to a scale of 70 pc arcsec⁻¹ at the galaxy. In spite of its proximity, there is only limited information about the stellar population in the vicinity of its active nucleus.

Lynds et al. (1991) reported a high contrast structure surrounding the nucleus in a Hubble Space Telescope (HST) Wide Field and Planetary Camera 2 (WFPC2) optical image of NGC 1068 in the continuum (centered at 5470 Å) interpreted as a stellar cluster. Thatte et al. (1997) reported the detection of this stellar cluster from stellar CO absorption features in the H- and K-band spectra and estimated its size as FWHM ≈ 50 pc and age in the range 5–16 × 10⁸ yr. Crenshaw & Kraemer (2000) have analyzed two long-slit spectra obtained with the HST Imaging Spectrograph (STIS) covering the inner ∼500 pc in the wavelength range 1150–10270 Å. They have found that the continuum in this inner region has contributions from two components: electron-scattered light from the hidden nucleus and stellar light. They conclude that the electron-scattered light dominates in the UV, while stellar light dominates at optical wavelengths. They found that the scattered light is strongest in a hot spot at ≈03 northeast of the nucleus, being enhanced also in regions of strong narrow-line emission. The stellar component was found to be more concentrated southwest of the hot spot, approximately coincident with the nucleus, being dominated by an old (age ⩾ 10⁹ yr) stellar population.

Crenshaw & Kraemer (2000) also support the presence of a nuclear stellar cluster ∼200 pc in extent and argue that it cannot be too young (in the sense of having blue ionizing stars), because the light concentration seen in the optical continuum is not seen in the UV (Macchetto et al. 1994). The UV Faint Object Camera HST images of Macchetto et al. (1994) resemble instead the narrowband [O iii] optical images, consistent with the interpretation that they are due to scattering, which is also supported by the fact that this light is highly polarized. It was the discovery of this polarized light emission, showing a Seyfert 1 spectrum, that led to the unified model of active galactic nuclei (AGNs) (Antonucci & Miller 1985; Antonucci 1993). Another argument against the presence of ionizing stars is the absence, in the UV spectra of Crenshaw & Kraemer (2000), of stellar absorption lines characteristic of the atmosphere of such stars.

More recently, Emsellem et al. (2006), using optical observations with the SAURON integral field spectrograph, have claimed that this galaxy is a candidate for a stellar velocity dispersion (σ_⋆) "drop" in the central ∼2 arcsec: σ_⋆ decreases to ∼100 km s⁻¹ in this region, being surrounded outward by higher values in the range 150–200 km s⁻¹. This σ_⋆-drop has been interpreted as the result of central gas accretion followed by an episode of star formation.

Davies et al. (2007) have also observed NGC 1068 and other nearby AGNs using the near-IR integral field spectrograph SINFONI at the Very Large Telescope (VLT). In most of the galaxies, they found signatures of circumnuclear starbursts with ages of 10–300 Myr (i.e., no longer ionizing). In the case of NGC 1068, they confirm the presence of a σ_⋆-drop, and also report an excess of stellar continuum within the inner 1'', which they attribute to the presence of a nuclear disk that has experienced recent star formation. This would also explain the decrease in velocity dispersion.

Martins et al. (2010) performed stellar population synthesis on the nuclear and extended regions of NGC 1068 by means of near-infrared spectroscopy (0.8–2.4 μm) using the NASA Infrared Telescope Facility SpeX-XD data. A 08 × 15'' slit oriented in the north–south direction was employed. Their main result is that traces of young stellar population are found at ∼100 pc south of the nucleus. They also report a ∼25% contribution of a power-law continuum at the nucleus and surrounding regions, with excess contribution (up to 60%) of this component at 100–150 pc to both sides of the nucleus, associated by them with regions where new stars are being formed. They also report that hot dust (1000 K ⩽ T ⩽ 1400 K) has an important contribution to the nuclear region while cold dust appears mostly to the south. A significant contribution of an intermediate-age stellar population to the continuum, especially in the inner 200 pc, was also found.

With the goal of mapping and quantifying the processes of feeding and feedback in the active nucleus of NGC 1068, as well as its stellar population, we have analyzed two-dimensional spectroscopic data obtained with the Gemini North Near-Infrared Integral Field Spectrograph (NIFS) in the J, H, and K bands. In the present paper, we report the mapping and analysis of the stellar population and near-infrared continuum within the inner 180 pc radius of NGC 1068 at a spatial resolution of ≈8 pc. We also include a brief analysis of the stellar kinematics, with the goal of relating it to the stellar population characteristics.

This paper is organized as follows. In Section 2, we discuss the observations and data reduction procedures, in Section 3 we derive the stellar kinematics, in Section 4 we obtain the properties of the stellar population, in Section 5 we discuss the results, comparing the stellar population properties with the derived kinematics and the molecular gas distribution, and in Section 6 we present our conclusions.

Spectra of NGC 1068 were obtained with NIFS (McGregor et al. 2003) on the Gemini North telescope in 2006 December under program GN-2006B-C-9. NIFS was used with the ALTAIR facility adaptive-optics system in its natural guide star mode. The compact Seyfert nucleus of the galaxy was used as the adaptive-optics reference object. A standard star observed under the same conditions as NGC 1068 has a spatial profile with an FWHM of 012 at both H and K, corresponding to a spatial resolution of ∼8 pc at the galaxy.

NIFS is an image-slicer integral-field spectrograph with a square field of view (FOV) of 3'' × 3'', divided into 29 slitlets each 0103 wide with a sampling of 0044 along each slitlet. The instrument was set to a position angle on the sky of 300° to align the slitlets approximately with the axis of the radio jet and ionization cone. This provided finer spatial sampling along the jet axis and coarser spatial sampling across the jet.

Observations at the Ks grating setting were obtained on 2006 December 12 (UT). Similar observations at the Kl grating setting were obtained on 2006 December 13 (UT), in the H band on 2006 December 14 (UT), and in the J band on 2006 December 15 (UT). The same orientation, observing procedure, exposure time, and standard stars were used on each of these nights. The spectra used in the present study cover the H band from 1.48 μm to 1.79 μm and the Kl band from 2.1 μm to 2.5 μm, both at a spectral resolving power of ∼5290. The spectra in other wavelength bands will be discussed in future papers.

Each data set consisted of multiple 90 s exposures. An offset sky position was recorded first and was followed by nine galaxy fields obtained on a 3 × 3 frame grid centered on the galaxy nucleus. The frame offset was 10, so frames were recorded at offsets of −10, 00, and +10 with respect to the nucleus along and perpendicular to the slitlets. This resulted in a full FOV 50 × 50 on the sky with maximum exposure in the central 30 × 30 and less exposure in the corners of the full field. The exposure grid was repeated nine times at each grating setting for a total exposure time on any one of the nine fields of 810 s. No pixel-level dithering of the field positions was employed. The nucleus of NGC 1068 was within each field, and so was subsequently used to spatially register each exposure. An argon/xenon arc exposure was obtained immediately after the galaxy exposures were completed. The A2 telluric standard star HIP 5886 was observed immediately prior to the NGC 1068 observation, and the A0 star HIP 18863 was observed immediately after the NGC 1068 observation. Both were selected so that they were observed through approximately the same atmospheric column as the NGC 1068 observations.

The full FOV and orientation of our NIFS observations are illustrated by the white square over an HST image through the filter F606W in the left panel of Figure 1, where we also show, for reference, within our FOV, a narrowband [O iii]λ5007A image (top right panel of Figure 1) showing the known ionization cone (Schmitt et al. 2003), as well as the radio jet (bottom right panel).

Zoom In Zoom Out Reset image size

The data reduction was performed using the GEMINI IRAF NIFS package and followed standard procedures, as described in Storchi-Bergmann et al. (2009). The nine sky exposures were first median-combined and dark subtracted to form an average sky frame. Individual galaxy exposures were then dark subtracted and the average sky exposure was subtracted. The 29 two-dimensional slitlet spectra for each exposure were cut from the raw data image and flatfielded. Known bad pixels identified from the flat field and dark exposures were corrected by linear interpolation. Each two-dimensional slitlet spectrum was then transformed to a rectilinear coordinate grid with horizontal wavelength and vertical spatial axes using transformations derived from the arc lamp and spatial calibration mask exposures. These transformed spectra were then stacked in the second spatial dimension to form a three-dimensional data cube with two spatial and one spectral axes.

The individual data cubes were then registered spatially by collapsing each data cube in the spectral direction to form a continuum image, determining the centroid of the NGC 1068 nucleus in each image, and then spatially shifting and median-combining all galaxy cubes using the IRAF imcombine task. The resulting mosaicked data cube was corrected for telluric absorption by multiplying each spectrum within the cube by a one-dimensional correction spectrum derived from the one-dimensional spectrum of the telluric standard star, HIP 5886, extracted from its data cube. A similar large-aperture one-dimensional spectrum of the total light of HIP 5886 was then used to flux calibrate the mosaicked galaxy data cube. The intrinsic spectrum of the standard star was set to the level of its Two Micron All Sky Survey (2MASS) magnitude in the wavelength band of the spectrum, and its shape was set to that of a 8600 K blackbody, which has the same 2MASS J − K color as HIP 5886.

In order to have approximately square spaxels, we have then binned together two pixels together along the slices. The resulting calibrated data cube contains 2376 spectra covering an angular region 45 × 50 (along the x-axis and the y-axis, respectively) at a sampling of 0103 × 0092, which corresponds to ∼315 × 350 pc² (or a radius of ∼180 pc about the nucleus) and ∼7.0 × 6.4 pc², respectively, at the galaxy.

2.1. Data Analysis

2.1.1. Line-of-sight Velocity Distributions

In order to map the stellar kinematics, we have obtained the line-of-sight velocity distributions (LOSVD) by fitting the absorption spectra in the H and K bands using the penalized pixel-fitting (pPXF) method of Cappellari & Emsellem (2004), following the procedure described in Riffel et al. (2008). In summary, the pPXF method finds the best fit to a galaxy spectrum by convolving stellar template spectra with a given LOSVD [L(v)], represented by Gauss–Hermite series (e.g., van der Marel & Franx 1993; Gerhard 1993; Riffel 2010):

$$\begin{equation} L(v)=\frac{e^{-(1/2)/y^2}}{\sigma _*\sqrt{2\pi }}\left[1+\sum _{m=3}^M h_{m*} H_m(y) \right], \end{equation} \tag{ 1 }$$

where y = (v − V_*)/σ, V_* is the stellar centroid radial velocity, σ_* is the velocity dispersion, v = c ln λ, and c is the speed of light. The H_m are the Hermite polynomials and h_m* are the Gauss–Hermite moments (Cappellari & Emsellem 2004). The pPXF routine outputs the stellar centroid velocity (V_*), velocity dispersion (σ_*), and the higher-order Gauss–Hermite moments h_3* and h_4*. As discussed in Riffel et al. (2008) and Winge et al. (2009), the use of a library of stellar templates, instead of a single stellar spectrum, is fundamental for a reliable derivation of the stellar kinematics. In the K band, we have used as templates the spectra of 60 stars from the Gemini library of late spectral type stars observed with the Gemini Near-Infrared Spectrograph (GNIRS) IFU and NIFS (Winge et al. 2009). In the H band, we do not have a large number of templates available, and we have thus used the spectra of five late-type stars with public NIFS data in the Gemini archive. The H-band spectra of these stars are shown in Figure 2, together with the NGC 1068 nuclear spectrum.

Zoom In Zoom Out Reset image size

2.1.2. Spectral Synthesis

In order to map the distribution of the continuum and stellar population properties, we have used the starlight spectral synthesis code (Asari et al. 2007; Cid Fernandes et al. 2005b, 2005a) adapted to spectra in the near-IR (Riffel et al. 2009a). This code fits a galaxy's spectrum by a combination of elements from a "spectral base" (described below), giving as output the weights of each base element in the synthetic spectrum (Cid Fernandes et al. 2009; Riffel et al. 2009a).

The spectral base was constructed with the most recent Evolutionary Population Synthesis models of Maraston & Strömbäck (2011, see also Maraston 1998, 2005) following Riffel et al. (2009a), and comprise single stellar populations (SSPs) synthetic spectra covering 12 ages (t = 0.01, 0.03, 0.05, 0.1, 0.3, 0.5, 0.7, 1, 2, 5, 9, and 13 Gyr) and solar metallicity. We did not include a stellar population younger than 10 Myr as we verified it was not necessary because it becomes degenerated with a power-law component (Storchi-Bergmann et al. 2000) also included in the synthesis. The power law (F_ν∝ν^−1.5) was used to represent the possible contribution of a featureless continuum from the active nucleus (e.g., Cid Fernandes et al. 2004). We also included blackbody functions for temperatures in the range 700–1400 K (one for each 100 K interval) in order to account for possible contributions from dust emission.

Following Cid Fernandes et al. (2004), we have binned the contribution of the SSP x_j elements into stellar population components of four age ranges t_j: young (x_y: t_y ⩽ 50 Myr); young-intermediate (x_yi: 0.1 Gyr ⩽ t_yi ⩽ 0.7 Gyr), intermediate-old (x_io: 1 Gyr ⩽ t_io ⩽ 2 Gyr), and old (x_o: 5 Gyr ⩽ t_o ⩽ 13 Gyr). We present the synthesis results in two forms: percent flux contribution of each SPC at λ = 2.1 μm (x_j) and percent mass contribution to the total mass (m_j).

As described in the previous section, we have used the pPXF method of Cappellari & Emsellem (2004) to fit the absorption spectra in the H and K bands in order to obtain the LOSVD. In our previous studies (e.g., Riffel et al. 2008, 2010, 2011), we have used only the K-band spectra for the derivation of the stellar kinematics. But in the case of NGC 1068, the K-band spectra within ≈03 from the nucleus are strongly diluted by a red continuum that can be attributed to emission by hot dust from the circumnuclear torus, as in the case of NGC 4151 (Riffel et al. 2009c). In addition, the stellar K-band absorption features (CO bands) toward the northeast, where the ionization cone is located (Crenshaw & Kraemer 2000), are contaminated by strong emission from the [Ca viii]λ2.322 μm coronal line. In order to avoid these contaminations present in the K-band spectra, we have used the H-band spectra to derive the stellar kinematics for angular distances from the nucleus smaller than 03. For the other regions, we used the K-band spectra whenever the contamination by the [Ca viii]λ2.322 μm coronal line was negligible, as the stellar absorption features are stronger in the K band than in the H band, resulting in a better constrained stellar population, particularly for the outermost regions. In order to "automatically" select which band to use for the derivation of the stellar kinematics, we used the flux of the [Ca viii]λ2.322 μm coronal line: in regions with fluxes of this line higher than 3σ (σ here meaning the root mean square deviation of the continuum flux in the neighboring continuum), we used the H-band spectra, while in other regions we used the K-band spectra.

In the top and bottom panels of Figure 3, we show representative fits (in red) compared to the observed spectra (in black), extracted at locations identified in the resulting velocity dispersion map: the nucleus, 05 northeast, 10 southeast, and 15 northwest of the nucleus. The region where we have used the H-band spectra for the derivation of the kinematics is shown in the center right panel of Figure 3 (shaded region in the figure).

Zoom In Zoom Out Reset image size

The left panel of Figure 4 shows the centroid velocity map, from which we subtracted the systemic velocity of V_LSR = 1168 ± 15 km s⁻¹ (corrected for the motion relative to the local standard of rest (LSR)), adopted as the average stellar velocity within an aperture of 05 × 046 centered on the nucleus. The centroid velocity map shows a rotation pattern, with the kinematic major axis approximately along the east–west direction, at P.A. ≈ 80°, with blueshifts to the east and redshifts to the west, reaching a projected velocity amplitude of 45 km s⁻¹. One-dimensional cuts showing the centroid velocities along two "pseudo-slits" with width 03 oriented along the galaxy major and minor axes (P.A. ≈ 80° and ≈ − 10°, respectively) are shown in Figure 5. The pseudo slit along the major axis is illustrated in the central left panel of Figure 3 by the dashed lines.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The right panel of Figure 4 shows the velocity dispersion σ_⋆ map. The highest values reach σ_⋆ ≈ 190 km s⁻¹, observed right at the nucleus and in patchy regions farther than 1'' from the nucleus toward the borders of the field. Within a region between ≈02 and ≈1'' from the nucleus (corresponding to 14 pc and 70 pc, respectively, although more extended along east–west than along north–south), σ_⋆ decreases to 125 km s⁻¹. One-dimensional cuts showing the velocity dispersions along the two "pseudo-slits" passing through the galaxy major and minor axes are shown in Figure 5.

In order to calculate the errors in the centroid velocities and velocity dispersions, we have performed Monte Carlo simulations in which Gaussian noise was added to the spectra. The corresponding error bars are illustrated in the one-dimensional cuts of Figure 5, showing that the errors both in the centroid velocities and in the velocity dispersions range from 10 to 20 km s⁻¹.

As described in Section 2, we have used the starlight spectral synthesis code adapted to spectra in the near-IR (Riffel et al. 2009a) to fit the galaxy's spectra by a combination of templates from a base including blackbody and power-law continua and single age stellar populations, giving as output the contribution of each base element to the synthetic spectrum (Cid Fernandes et al. 2009; Riffel et al. 2009a).

We show in Figure 6 the results of the synthesis for two spectra (extracted within our aperture of 01×009), one at the nucleus (top panel) and the other at 1'' north of it (bottom panel). The observed (continuous line) and synthetic (dashed line) spectra were normalized to unity at 2.1 μm.

Zoom In Zoom Out Reset image size

The spatial distribution of the stellar population components is shown in Figure 7. The top panels show the percent contributions to the light at λ2.1 μm of the young, young-intermediate, intermediate-old, and old age components, from left to right: x_y, x_yi, x_io, and x_o, while the bottom panels show the corresponding mass contributions m_y, m_yi, m_io, and m_o. We have masked out from the figures the results of the synthesis in the regions for which adev (see below) is smaller than 2% in order to show only the results based on robust fits.

Zoom In Zoom Out Reset image size

The flux contribution to the continuum within the inner 180 pc of NGC 1068 is dominated by the young and intermediate age stellar populations. The flux contribution of the young population x_y, as well as its contribution to the total mass m_y, shows a "hole" of radius ≈1'' around the nucleus covering the east–north–west directions, but extending almost to ≈2'' toward the southwest–south. The x_y and m_y are distributed around this "hole," with enhanced contributions in a partial ring around the nucleus from the northeast–east toward the south and to the west–southwest. The young-intermediate age contribution x_yi and m_yi also shows a "hole" around the nucleus but with approximately half the radius of the one observed in x_y and m_y, and with similar extent in all directions. The southwest–south extension of the "hole" seen in x_y and m_y is not present in x_yi and m_yi, which show enhanced contributions in this region—a clump between ≈1 and 2'' to the right (southwest) of the nucleus in Figure 7—where m_yi ≈ 100%. The remaining regions of enhanced contributions of x_yi and m_yi seem to be anti-correlated with those of enhanced x_y and m_y, being observed mostly to the northeast through the north and west.

The contribution of the old stellar population, x_o and m_o, is restricted to within ≈05 (35 pc) from the nucleus, where the contribution to the total mass m_o is dominated by this component. The contributions of the intermediate-old stellar population x_io and m_io are small and avoid the nuclear region, being somewhat co-spatial with the regions of enhanced contribution of the young stellar population.

Figure 8 shows the contribution of the power-law and blackbody components, together with the reddening derived from the stellar population synthesis. The rightmost panel shows the distribution of the adev, which is the percent mean deviation of the fit – |O_λ − M_λ|/O_λ × 100, where O_λ is the observed spectrum and M_λ is the fitted model (Cid Fernandes et al. 2004, 2005a)—and gives the robustness of the fit. Values of adev ≲ 2% correspond to good fits, while higher values, observed toward the corners of the field, indicate less robust fits.

Zoom In Zoom Out Reset image size

Figure 8 shows that the contribution to the flux within 05 of the nucleus is dominated by the blackbody components with 700 K ⩽ T ⩽ 800 K, which also show a significant contribution to a radius ≈1''. The blackbody components with 900 K ⩽ T ⩽ 1000 K also have some contribution closer to the nucleus, while for T ⩾ 1100 K the contribution occurs even closer to the nucleus but is of negligible strength. Figure 8 also shows that the reddening obtained from the synthesis is only significant within the inner ≈05, where it reaches E(B − V) = 1.5 mag. The power-law component is extended toward the south–southeast (right in Figure 8), while the vertical structure is an artifact probably due to scattering of the strong nuclear source along the NIFS mirror slices.

5.1. Stellar Kinematics

5.2. Stellar Population

5.2.1. Young Stars and Correlation with H₂ Ring

Figure 7 shows that the main stellar population components contributing to the flux within the inner 180 pc of NGC 1068 are the young and the young-intermediate ones. Within these components, we can verify from the output of the synthesis that the age bins which contribute with most of the flux are those of 30 Myr and 300 Myr, respectively. Thus, it can be concluded that two main episodes of recent star formation have occurred in the nuclear region. The first burst occurred about 300 Myr ago, reaching as close as 05 (35 pc) from the nucleus to the southwest, west, and north of the nucleus. The location of this burst may extend even closer to the nucleus, but the dominance of the BB component in this region may preclude the detection of this component there. Particularly conspicuous is a region between 08 (60 pc) and 16 (100 pc) from the nucleus to the southwest, where it is the only stellar population component. The second burst occurred about 30 Myr ago, mainly in an approximately ring-like structure with an average radius of ≈100 pc.

In Vale et al. (2012), we have obtained emission line flux distributions from the same data we are using here. The flux in the molecular hydrogen emission line H₂λ2.1218 μm, shown in the left panel of Figure 9, is distributed in a ring—first shown by Müller-Sánchez et al. (2009)—similar to what is observed in NGC 4151 (Storchi-Bergmann et al. 2010). In Figure 9, the H₂ ring (left panel) can be compared with the young stellar population distributions x_y and m_y (central and right panels, respectively). This figure shows that x_y and m_y are enhanced at locations covered by the H₂ ring, and, in particular, the region of strongest H₂ emission is the region with the largest contribution of the young stellar component.

Zoom In Zoom Out Reset image size

The H₂ ring can be identified with the circumnuclear disk of molecular gas first reported by Schinnerer et al. (2000). In F. K. B. Barbosa et al. (in preparation), we show that the H₂ emitting gas in this ring seems to be expanding in the plane of the galaxy, confirming previous suggestions by Galliano & Alloin (2002) and results from observations with the instrument SINFONI on the VLT (Müller-Sánchez et al. 2009). Recent observations of the molecular gas emission in the millimeter spectral region (e.g., lines of CO, HCN) by Krips et al. (2011) also show this ring as well as its expansion, which, combined with the measured line ratios, have led these authors to propose a scenario in which part of the molecular gas is radially blown outward as a result of shocks. Krips et al. (2011) suggest that the shocks originate in the interaction of the AGN jet with the molecular gas.

Krips et al. (2011) also point out that an expanding shell of shocked gas can originate in stellar winds from the evolution of a super stellar cluster (although they do not favor this interpretation). This is the hypothesis put forth by Sakamoto et al. (2006) to explain the presence of "superbubbles" observed in molecular gas emission in starburst galaxies such as NGC 253 and M 82 (Matsushita et al. 2000). In the case of AGNs, the presence of superbubbles has been reported by Lipari and collaborators in a number of galaxies, examples being Mrk 231 (e.g., Lipari et al. 2005 and references therein) and IRAS 04505-2958 (Lipari et al. 2009). These authors suggest that broad absorption line (BAL) systems seen in the spectra of AGNs could be due to these superbubbles—giant expanding gas shells which originate in circumnuclear starbursts. Lipari & Terlevich (2006) propose an evolutionary scenario in which massive circumnuclear starbursts in the vicinity of AGNs produce supernovae and hypernovae which lead to these large-scale expanding shells, which are often obscured by dust. Later in the evolution, as the starburst ages and the shells expand, the AGN becomes the main source of ionization of the circumnuclear gas. This scenario is shown to reproduce many of the observed correlations among spectral properties of infrared-luminous AGNs. For lower luminosity Seyfert galaxies, Storchi-Bergmann et al. (2001) also propose such an evolutionary scenario to explain the excess of young- and intermediate-age stellar populations within the inner kiloparsec: the gas which flows in first triggers star formation in the circumnuclear region, and the nuclear activity is then subsequently triggered.

Thus, one possibility to be considered is that the H₂ ring is actually a superbubble, for which the strongest emission comes from the galaxy plane, where there is more gas to be swept by the bubble. This is supported by the fact that NGC 1068 sits next to Mrk 231 in the diagram of Davies et al. (2007) which relates the AGN luminosity with the circumnuclear starburst age. Although the expansion velocity of the H₂ ring—under 100 km s⁻¹ (F. K. B. Barbosa et al., in preparation)—is much lower than in Mrk 321 superbubbles, in nearby starbursts, observations in the millimetric spectral region show similar expansion velocities of a few tens of km s⁻¹. The morphology and sizes (100–200 pc) of superbubbles in nearby starbursts are also similar to those of the molecular ring of NGC 1068. One difference is the mass; while in the nearby starbursts the gas masses are of the order of 10⁶ solar masses (e.g., in NGC 253), the molecular gas mass in the inner ≈100 pc of NGC 1068 is estimated to be an order of magnitude larger (Müller-Sánchez et al. 2009 and references therein).

In the case of the nearby starburst galaxy NGC 253, Sakamoto et al. (2006) have calculated the necessary power to drive the outflow and concluded that a 10⁶ M_☉ super star cluster could be the origin of such bubbles, created by successive supernova explosions, which can have a duration of ≈30 Myr, arguing that these bubbles are the end products of star formation. Having a mass an order of magnitude larger than the super bubbles, it is not clear that the H₂ ring in NGC 1068 has a similar origin.

It is interesting to note that the apparent center of the H₂ ring does not coincide with the nucleus, but is located in the region ≈1'' to the southwest where the 300 Myr stellar population is particularly enhanced. If the H₂ ring were actually a superbubble originated from stellar winds and supernova explosions, then a cluster located there—in the center of the ring—could have originated it. A problem with this scenario is the fact that the ring is expanding and, even at a lower limit of the expansion velocity of 10 km s⁻¹, it would reach its present radius of about 100 pc in only 10 Myr. As supernovae appear no later than ≈50 Myr, a cluster of age 300 Myr seems too old to have originated the H₂ ring in NGC 1068.

And why would x_y be enhanced in the H₂ ring? Why are the regions with strongest H₂ emission (e.g., at ≈15 east of the nucleus) the ones with highest concentration of young stars? Could an expanding bubble have triggered the 30 Myr episode of star formation? This seems unlikely due to the fact pointed out above that in 10 Myr the ring moves at least ≈100 pc, and this is not consistent with the present radius of the ring, of ≈100 pc, and the age of the stars in the ring, of 30 Myr.

On the other hand, the 30 Myr age is consistent with the presence of supernova explosions originated in this same population, which could be the origin of the shocks inferred by Krips et al. (2011) from the molecular line ratios. The presence of such supernova explosions could also explain the high velocity dispersion observed in the region of strongest H₂ emission at 15 east of the nucleus (F. K. B. Barbosa et al., in preparation), and could be the origin of the apparent expansion of the H₂ ring. A spatial coincidence between a 10 and 50 Myr stellar population and strong H₂ emission was also found by Martins et al. (2010) at 100 pc south of the nucleus in a recent long-slit study of the stellar population of NGC 1068.

5.2.2. Older Age Components and Continua

A comparison between the stellar velocity dispersion (σ_⋆) distribution (Figure 4) and the distribution of x_yi and x_io in Figure 7 reveals that, in most of the region where there is a drop in σ_⋆, the dominant stellar population is the intermediate-age one (yi and io). We tentatively conclude that this intermediate-age stellar population is related to the σ_⋆ drop, as we have found already for Mrk 1066 (Riffel et al. 2010) and Mrk 1157 (Riffel et al. 2011). Similar rings seem to be common in the central region of active galaxies (e.g., Barbosa 2006; Riffel et al. 2008, 2009b) and can be interpreted as being due to colder regions with more recent star formation than the underlying bulge. This interpretation is in agreement with the previous analysis of SINFONI data by Davies et al. (2007), who have estimated as 200–300 Myr the age of the stars in the nuclear compact disk of NGC 1068. In a previous work, Thatte et al. (1997) claimed to have found a stellar core with intrinsic size scale of 45 pc. Using the estimated mass-to-light ratio and assuming some contamination from old stars, they propose an age of 500 Myr for a young age component in this region, in agreement with our results.

We have found a small contribution from the old component, x_o, to the flux in the inner ≈05 (35 pc), but even though the flux contribution in the near-IR is not larger than ≈25%, the mass contribution is close to ≈100%. The presence of an old stellar component around the nucleus, with age ⩾ 10⁹ yr, has already been suggested by Crenshaw & Kraemer (2000), on the basis of long-slit STIS spectra, as the dominant light component in the optical part of the spectrum.

Figure 8 shows the distribution of the power-law and blackbody components. The power-law component does not contribute much at the nucleus, being spread around it. The vertical line is an artifact, which also appears in the lowest temperature blackbody component; it is an instrumental effect, due to scattered nuclear light along the vertical mirror slices of NIFS. We attribute the extra-nuclear contribution of the power-law component to scattered light from the nucleus, which has been found to be a strong component of the continuum in the ultraviolet (Crenshaw & Kraemer 2000).

As pointed out above, the blackbody components dominate completely the light at the nucleus, in particular the ones with temperatures in the range 700–800 K (Figure 8), comprising almost 100% of the flux. We interpret them as originating in the dusty torus surrounding the AGN, which is unresolved by our observations (Marco & Alloin 2000; Jaffe et al. 2004; Raban et al. 2009; Martins et al. 2010). In Figure 10, we compare the continuum emission in the K band (left) with the sum of all blackbody components (center) and with the H₂ flux distribution (right). This comparison—which shows that the continuum emission and BB flux distributions are similarly concentrated—gives support to the identification of the peak of the K-band continuum with the blackbody components contribution. Some contribution of the blackbody components is also seen beyond the inner arcsecond, which we attribute to the presence of faint dust emission or to further contribution from scattered nuclear light spread in the nuclear region. The comparison between the central and right panels of Figure 10 shows that some BB emission is seen also along the west part of the ring (to the top and right in Figure 10).

Zoom In Zoom Out Reset image size

We have mapped the stellar kinematics and age components, as well as the contribution of power-law and blackbody emission components in the inner 180 pc radius of NGC 1068 at a spatial resolution of 8 pc using integral field spectroscopy in the near-infrared. The main results of this paper are as follows.

• 1.  
  The stellar kinematics are well described by rotation with a projected amplitude of ≈45 km s⁻¹ and major axis at the position angle of ≈80°.
• 2.  
  The stellar velocity dispersion σ_⋆ shows a drop from 165 km s⁻¹ to 125 km s⁻¹ within the inner 1''–70 pc at the galaxy.
• 3.  
  The dominant stellar populations over the inner 180 pc radius are the young (mainly 30 Myr) and intermediate-age (mainly 300 Myr) ones.
• 4.  
  Within the inner 05 (35 pc), the blackbody components with 700 K ⩽ T ⩽ 800 K dominate the flux and are attributed to the dusty torus (unresolved by our observations) surrounding the active nucleus.
• 5.  
  Within the inner 35 pc, the old component (age ⩾ 5 Gyr) dominates the contribution of the stellar population in terms of mass fraction.
• 6.  
  The distribution of the youngest (30 Myr) stellar population is co-spatial with an expanding ring of H₂ emission at ≈100 pc from the nucleus. We speculate that supernova explosions from this young stellar population could be the origin of the apparent expansion of the ring.
• 7.  
  In the region of the σ_⋆-drop, the dominant stellar population is the intermediate-age (mainly 300 Myr) one.

NGC 1068—the best studied Seyfert 2 galaxy—is thus another case of an active nucleus surrounded by recent star formation in a circumnuclear ring, supporting the evolutionary scenario for AGNs (Storchi-Bergmann et al. 2001; Lipari & Terlevich 2006; Davies et al. 2007). In such a scenario, gas from the external regions of the galaxy (or even acquired from outside) flows to the nuclear region, first triggering star formation in the inner few hundred parsecs, and subsequently reaches the nucleus to feed the supermassive black hole, triggering the nuclear activity.

We thank the anonymous referee for valuable suggestions which helped to improve the paper. Based on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the Science and Technology Facilities Council (United Kingdom), the National Research Council (Canada), CONICYT (Chile), the Australian Research Council (Australia), Ministério da Ciência e Tecnologia (Brazil), and Ministerio de Ciencia, Tecnologia e Innovación (Argentina). This work has been partially supported by the Brazilian institution CNPq. R.R. thanks the Brazilian funding agency FAPERGS (process 11/1758-5) for financial support.