© ESO, 2011

1. Introduction

NGC 1068, at a distance of 15 Mpc (Paturel et al. 2002), (z = 0.003786, 1″ = 73   pc with H₀ = 75 km s^-1 Mpc^-1), is a key object for the study of the inner structures of Seyfert galaxies. Using a 10 m class telescope coupled to adaptive optics, we are able to probe scales on the order of a few parsecs in the near-IR.

The standard model of AGN (Elvis 2000) assumes a very massive black hole with an accretion disk surrounded by a dusty molecular torus, whose direction with regard to the line of sight is responsible for the discrepancy between Seyfert I and Seyfert II nuclei. The size of the central source is still an open question (Ferguson et al. 1997; Nenkova et al. 2002, 2008). Its near-IR luminosity is dominated by thermal radiation from hot dust at the sublimation temperature (Thatte et al. 1997; Gratadour et al. 2003). Stars may contribute as well because Thatte et al. (1997) found that at least 7% of the nuclear luminosity is caused by a nuclear stellar core for NGC 1068.

Fig. 1 Radio jet contour map (left) from Gallimore et al. (1996b), [O iii] emission (center) from Macchetto et al. (1994) and mid-IR contour (right) from (Bock et al. 2000) over the deconvolved M′-band image (Gratadour et al. 2006).

Fig. 2 Nodule identification and slit orientation on M′-band image (left) and its deconvolved image (right). The distance between the ticks inside the slit is 0.1″ i.e. about 7.5 pc. In white, the 47° PA slit aligned with the nodules. In green, the 12° PA slit. The colors are in arbitrary units and are chosen so that the structures mentioned in this article appear clearly.

Several highly ionized forbidden emission lines, called coronal lines, can be observed in NGC 1068 (Oliva & Moorwood 1990; Geballe et al. 2009). Their ionization potentials (>100 eV) imply an excitation by hard photons from the AGN or an excitation by shocks from the NGC 1068’s radio jet (Gallimore et al. 1996b,a; Axon et al. 1998). At a lower spatial resolution than NaCo, both mechanisms have been proposed. Kraemer & Crenshaw (2000b) found that the central source is able to produce coronal lines. However, they also found evidence for an additional ionizing source, which could be the UV-X radiation generated by fast shocks (~1000 km s^-1). Geballe et al. (2009) also favored a collisional excitation for the ionized gas with near-IR data. Indeed, the jet interacts with the interstellar medium (hereafter ISM) very close to the core of NGC 1068. Bicknell et al. (1998) found that the deviation of the jet, at 0.3″ from the core, is caused by its interaction with a giant molecular cloud. Sánchez et al. (2009) confirmed that the “northern tongue” correlates with the mid-IR emission and coincides with a knot of radio continuum, providing direct evidence of interaction between the jet and the molecular cloud.

Rouan et al. (2004) and Gratadour et al. (2005, 2006) observed the nearest regions from the core of NGC 1068 (<70 pc) at the diffraction limits of the VLT with NaCo. They discovered four elongated “nodules” north of the radio jet that form a wave-like pattern with the same orientation as the jet (Fig. 1). These structures are probably caused by instabilities in the ISM because of the compression by the jet (Fig. 1). Rouan et al. (2004) and Gratadour et al. (2003) also found that the color temperature of these nodules could be explained by the stochastic heating of very small dust grains by the central source. Moreover, Thatte et al. (1997) found evidence for stellar activity in these regions (up to 50 pc from the central source).

This work is based on spectra of these nodules at the L′- and K_s-bands obtained with NaCo. For the first time, coronal lines can be associated to structures seen in broad-band images. We present in Sect. 3 results from long-slit spectra with a slit crossing the central source, the structure IR-1b, and nodules IR-3 and IR-4. In Sect. 4 we present results from long-slit spectra with a slit aligned with the four nodules. We report evidence of stellar activity at the location of the nodules and provide constraints on the contribution of stellar activity to the luminosity of the closest regions to the core. We furthermore report the detection of several coronal lines at the location of the nodules. In Sect. 5 we compare our data to Cloudy simulations with different ionizing continua to identify the possible source responsible for the coronal lines at the location of the nodules.

2. Observations and data processing

Spectra of NGC 1068 were obtained during two different observing runs, hence, in the following, we will separate the data into two different data sets, namely DS1 and DS2. All these data were acquired with NaCo (Rousset et al. 2003; Lagrange et al. 2003) at the Nasmyth focus of YEPUN (VLT-UT4) using the long-slit spectroscopic mode. The visible wavefront sensor (WFS) of NaCo was used in its most accurate mode (14  ×  14 sub apertures) because the nucleus of this galaxy is bright enough.

2.1. Data acquisition

Data set 1 (DS1) is made of K-short spectra obtained in November 2002. Thirty-six spectra of NGC 1068 were recorded between two series of 15 spectra of HD 16835 (m_K = 6.705 from 2MASS archive), an F0 star that we used for photometric and point spread function (PSF) calibration as well as atmospheric transmission corrections. The auto-jitter mode was used so that each spectrum is randomly translated in the plane of the camera within a box of 8″  ×  8″ to ensure a good correction of any residual detector effect. The calibration source is almost as bright as the NGC 1068 core so that the AO system showed a comparable behavior in both cases. During the observations of NGC 1068, the Strehl ratio ranged from 0.4 up to 0.6, and during the reference star acquisition, the Strehl ratio ranged from 0.5 up to 0.6. The airmass ranged from 1 to 1.2 and the seeing was fairly good, about 0.8″  ±  0.2 during all observations of NGC 1068. The total integration time on NGC 1068 is 7200 s and the pixel scale was 27 mas. The 86 mas slit was used and oriented along the polar angle PA = 12°. As shown in Fig. 2, this positioning allows us to probe the infrared light of several structures to the north of the nucleus: structure IR-1b, IR-3 and IR-4. South of the central core, no identified structure is crossed by the slit and only the southern diffuse emission is probed.

Data set 2 (DS2) is made of K-short and L′-band spectra of the nodules IR-1 to IR-4 discovered in our previous imaging study (Rouan et al. 2004; Gratadour et al. 2006), obtained in September 2005. The 86 mas slit was oriented along PA = 47° and centered on IR-2. During this observing run, the atmospheric conditions were quite fairly poor, with an average seeing of 1.55″, a strong wind and a fast turbulence, leading to Strehl ratios varying from 0.1 to 0.4 at 2.2 μm and in turn to an average SR of 0.2 at the K_s-band, i.e. 0.4 at L′-band. Pixels of 27 mas and 54 mas were used at the K_s- and L′-band respectively. Eighty spectra of NGC 1068 were recorded in the K_s-band and 60 in the L′-band, representing a total integration time of 5600 s and 1200 s, respectively. In this case, BD-00 413, a G0 star, was used as a calibration star (m_K = 9.43 from the 2MASS archive).

2.2. Data reduction

The basic data reduction of the spectral images was similar to what was described in Gratadour et al. (2005) for classical imaging. After this we used dedicated spectroscopic procedures were used. First, the dispersion of the grism was not strictly linear over the whole detector and thus introduced a warp in the lines along the slit leading to a displacement of the exact position of a wavelength along the slit as large as 4 pixel between each edge of the detector. The induced warp can be evaluated by fitting a second-order polynomial on the lines of the argon calibration lamp spectrum for the K_s-band observations, and on atmospheric lines in the calibration star spectral images for the L′-band. The functions obtained for each coefficient of the second order polynomial were then fitted by a linear function. The spectral images were then over-sampled by a factor of ten so that the corrected positions became integer numbers. Finally, the images were resampled after correction to their original size.

The achieved spectral resolution has been evaluated on the corrected spectral image of the spectral lamp. All individual spectra along the slit were averaged and all identified lines of the calibration spectrum were fitted by a Gaussian distribution. The resulting spectral resolutions reached from 1720 to 1740 in the K_s-band over the entire wavelength range, which is very close to the expected spectral resolution of the instrument.

The correction of atmospheric effects was necessary as well. For that purpose, we directly compared the averaged processed spectra of the calibration star with a theoretical spectrum that consists of a continuum following a Rayleigh-Jeans distribution at these wavelengths. This way, we obtained a spectrum of the atmospheric transmission during the observations. Each spectrum extracted from the total spectral image of NGC 1068 was then corrected with this atmospheric transmission. Taking into account the airmass variation during the observations, the effect of the spatial variability of the atmosphere is lower than 8%. The effect of the temporal variability can be neglected because only 25 min separate the last frame of NGC 1068 and the first frame of the calibration star, while the exposure time for each object is one to two hours.

2.3. Photometric calibration

We used the average spectrum of the calibration star to calibrate the flux on our spectra. Only a fraction of the flux contained in the PSF passes through the slit. To estimate this fraction, we computed the PSF obtained with a Strehl ratio of 0.6 with an in-house AO end-to-end simulation code. This led to an aperture correction of 55% in a 70 pixel box around the slit in the spatial dimension. We thus summed 70 individual spectra of the star along the slit, corrected from detection (dark current, flat field, background, warp, and rotation) and atmospheric effects and used the result as a calibration spectrum containing 55% of the star luminosity. Given all sources of noise (airmass, temporal variability, and pointing error) as well as the accuracy of the estimated aperture correction, we put an upper limit of 10% on the error of the flux given below for DS1 and up to 20% for dataset 2 owing to the relatively high Strehl ratio variation during the observations.

3. Results from the 12° PA slit

Fig. 3 Black: spectra from DS1 for four different positions in the slit. Red: fits by one blackbody with temperatures summarized in Table 1.

3.1. Central continuum

K_s-band spectra extracted from DS1 by bins of 0.054″ (i.e. about 4 pc) were used to estimate the color temperature of the closest regions to the core. These spectra were smoothed by the convolution with a Gaussian of 50 pixel FWHM in the spectral dimension to reduce inaccuracies caused by potential emission lines, and then fitted by a blackbody spectrum. The temperature of the best fit gives us an estimate of the color temperature of the considered region. The precision reached on the estimated color temperature is better than 1% as shown in Fig. 3. North of the maximum emission the signal to noise ratio (SNR) on the continuum was sufficient for this fit up to about 20 pc from the central source. South of the maximum, only 12 pc from the central source could be probed. The results are outlined in Table 1.

The equivalent temperature of the central core (950 K) is consistent with our previous results at a lower resolution (Gratadour et al. 2003), strengthening our previous interpretation: the very central source in the near-IR is the dust sublimation cavity around the central engine. Additionally, the enhanced resolution of these NaCo data reveals a peculiar evolution around the central core. While the equivalent temperature remains fairly high and is almost constant up to 20 pc to the north of the nucleus, it decreases more rapidly to the south. This disparity can be understood using Fig. 2. To the south the continuum is dominated by the scattered light from the central source because there is no structure crossed by the slit and assuming that the central core is almost point-like at 2.2 μm. To the north, because structure IR-1b is crossed by the slit, the continuum is the sum of the scattered light from the central source and the intrinsic continuum of source IR-1b. Without spatial deconvolution, it is difficult to separate both contributions, however, the fast decrease to the south seems to indicate that, in the north, at a distance farther than 15 pc, the continuum of source IR-1b becomes dominant.

It is very likely that this continuum is caused by hot dust distributed in clouds around the central engine. Measuring a constant equivalent temperature between 12 and 20 pc is quite puzzling if we assume that this dust is heated by the central source. In this dataset, we were unable to find evidence for star formation on the central spectrum (no CO band head). Thatte et al. (1997) argue for a relatively low stellar contribution at a lower spatial resolution. From our observation data, the central continuum seems to be entirely dominated by a 800 to 900 K continuum, which hardly reproduces using stellar light. Stochastic heating of very small dust grains could provide a good solution to the puzzle as already inferred for the northernmost IR structures in our previous work (Gratadour et al. 2005, 2006). However, an accurate confrontation with a numerical model is fairly difficult because of the above-mentioned contamination by the scattered light from the central source.

3.2. Emission line distribution

Fig. 4 Ks-band spectra extracted from various regions along the slit with identified IR emission lines. These spectra were renormalized to be displayed on the same plot. The line intensity variation is not absolute.

Using Fig. 2 we distinguished several regions along the slit for which we extracted the averaged spectrum:

• the southern region, spanning 42 pc to the south ofthe nucleus;

• the central region, spanning 10 pc around the central core;

• northern region A, spanning 10 pc around 15 pc to the north of the nucleus (i.e. in between the central core and IR-1b);

• northern region B, spanning 10 pc around 25 pc to the north of the nucleus (i.e. IR-1b);

• northern region C, spanning 20 pc around 40 pc to the north of the nucleus (i.e. empty/filamentary space between IR-1b and IR-3);

• the IR-3 region, spanning 10 pc around 50 pc to the north of the nucleus;

• the IR-4 region, spanning 10 pc around 60 pc to the north of the nucleus.

The K_s-band spectra of these regions are displayed in Fig. 4. In order to build this figure, we normalized all the spectra to the average intensity of the central core spectrum and shifted them in intensity by an arbitrary amount to visualize all of them in the same plot. Hence this figure is representative of the emission line intensity evolution with respect to the continuum intensity and the reader should not interpret the line intensity variation in an absolute manner but rather as a comparison to the continuum intensity. This representation also helps to evaluate the signal-to-noise ratio (SNR) obtained in each region because the noise is amplified by the same amount as the signal with this method.

As expected at this angular resolution (4 pc per bin), apart for some atmospheric residuals around 2.4 μm and beyond, the central region spectrum is totally featureless at the K_s-band, the emission being dominated by the continuum of the central source.

In the south, no emission line is clearly detected up to 42 pc from the nucleus. As previously argued, no identified structure is crossed by the slit in this region and these spectra show that the diffuse emission to the south of the nucleus contain no emission line.

In the north, along the slit, narrow emission lines of different species (H i Brackett γ at 2.18 μm,  [Ca viii]  and  [Si vii]  forbidden lines at 2.33 μm and 2.49 μm respectively) appear and their intensities evolve from 15 to 70 pc to the nucleus which agrees well with the morphology of the broad-band IR image. Indeed, emission lines appear at the location of structures identified in the deconvolved M′-band image (IR-1b, IR-3 and IR-4), whereas at 40 pc north, where no structure is identified, very little to no emission is detected at these wavelengths.

For the first time, thanks to the spatial resolution provided by NaCo, near-IR coronal emission can be precisely associated to structures observable in the near-IR broad-band images. It shows that coronal emission is not distributed homogeneously in a narrow region around the central core but rather in dense clouds spanning a region as wide as 70 pc from the nucleus.

Little molecular hydrogen emission is detected at 25 pc north from the nucleus and beyond. Only the H₂ 1-0 S(1) line at 2.12 μm is detected and it is quite weak. It follows the same distribution as the ionization lines, which is concentrated on the identified IR structures.

3.3. Emission lines properties

We measured the emission line properties (line width and intensity ratios) using Gaussian fits on the spectra of the various identified regions. A second-order polynomial was first fitted on the continuum for subtraction.

Because of the atmospheric residuals beyond 3.4 μm, the properties of the  [Si vii]  line are difficult to measure accurately, especially in regions A and B. However, by comparing the central and southern spectra to spectra of regions A and B, it is obvious that this line appears in the latter regions. Because the integral of a Gaussian is proportional to the product of its amplitude by its FWHM, we were able to estimate the  [Si vii]/[Ca viii]  ratio by comparing their amplitude and FWHM. The  [Ca viii]  line exhibits an almost constant FWHM of 450 km s^-1 comparable to the  [Si vii]  line width, as measured on IR-3 and IR-4. Hence the amplitude ratio in regions A and B gives us the intensity ratio for these regions. The emission lines properties for each regions are outlined in Tables 2–4.

4. Results from the 47° PA slit

Fig. 5 Spectrum from DS2 at the Ks-band. The spectrum is the sum of the four nodule spectra.

Fig. 6 Idem Fig. 5 at the L′-band.

4.1. Nodule identification

Plotting in Fig. 7 the maximum flux value of the [Si vii] line as a function of the distance from the central source, we see three peaks at about 30, 40, and 50   pc corresponding to the nodules IR-1, IR-2 and IR-3. The fourth nodule seems to be blurred because of the poor seeing, but obviously a source is present at about 60   pc from the central source because of the slow decrease up to 65 pc. This identification is impossible at the L′-band because of the poor SNR. As a consequence, to increase the SRN of our data, measurements in this work at the K_s- and L′-band have been made on the average spectra of the four nodules, hereafter called “averaged nodule spectrum”.

Fig. 7 [Si vii]/Brγ flux ratio as a function of the distance from the central source in blue, superimposed on the maximum intensity along the spatial dimension of the slit as a function of the distance from the central source in black (arbitrary units). The ratio is almost constant whatever the distance. In green, the position of the nodules identified.

Fig. 8 Spectra of the nodules at the Ks-band in arbitrary units so we can display them in the same plot.

4.2. Nodules temperature

The averaged nodule continuum has a peculiar shape. At the K_s-band the continuum is flat up to 2.3 μm, then has a positive slope. At the L′-band the positive slope is reproduced well by a blackbody at  ~500   K corresponding to hot dust emission. As shown in Fig. 9, such a blackbody alone is not able to reproduce the fairly flat shape in the K_s-band. Thus, we added a graybody (F_λ ∝ λ^βB_λ) component corresponding to stars seen through gas and dust.

Assuming a plane-parallel geometry and stars mixed with gas and dust in the nodules, the flux of the graybody component is , which can be approximated by . Assuming τ ∝ Q_abs ∝ λ^-1, the equations lead to F_λ ∝ λF₀, with F₀ the flux of stars, Q_abs the absorption coefficient of efficiency, and τ the optical depth of the medium.

Fitting the nodules spectrum with the sum of a graybody and a blackbody, we derived a dust temperature T_dust = 519 ± 3   K (in agreement with Gratadour et al. 2006) and for stars T_star = 7500 ± 1800   K (Fig. 9). We derived a stellar contribution of 18% from the fit to the total luminosity up to  ~50 pc around the nucleus. This is fairly consistent with the stellar contribution calculated by Thatte et al. (1997) i.e. at least 7% of the total nuclear luminosity at these distances.

Fig. 9 Fit of the spectra with a unique component at 500 K (green), and with a sum of a graybody and a blackbody (red).

4.3. Emission lines

4.3.1. H i lines

The Brγ emission line (2.166  μm at rest) is a ionizing field/stellar activity indicator. Brγ has been detected by many authors in the central parts of NGC 1068 (Kawara et al. 1989; Gratadour et al. 2003; Martins et al. 2010).

The Brα (4.08   μm at rest) is another stellar activity indicator and is already identified in the central parts of NGC 1068 (Kawara et al. 1989).

Considering a temperature of 7000 K for the gas (see Sect. 4.2) and a hydrogen density of 10³ cm^-3 (Lutz et al. 2000), we have an intrinsic ratio Brα/Brγ = 2.96 (Hummer & Storey 1987). Following Ho et al. (1990), the differential extinction between the Brγ and Brα lines in our spectrum is A_γ − A_α = 1.32 ± 0.10.

The total extinction at 4.05 μm is then A_α = 0.58 ± 0.04 leading to A_γ = 1.90 ± 0.14. The visible extinction is calculated using A_2.2   μm = A_V × 0.108 (Cox 2000). Assuming A_γ ≃ A_2.2   μm, then A_V = 17.6 ± 1.29.

Finally, we estimated the hydrogen density of the nodules: A_V = 5.3 × 10^-22 × N_H (Cox 2000) leading to N_H = (3.32 ± 0.24) × 10²²   cm^-2 considering N_H = N(H + H ⁺  + 2H₂) the column density of hydrogen atoms. Assuming a thickness of  ~5 pc for one nodule, the hydrogen density n_H ≃ (2.0 ± 0.1) × 10³ cm^-3, which is consistent with the assumed value from Lutz et al. (2000).

4.3.2. Ro-vibrationnal H₂ 1-0 S(1) line

Only one H₂ line is detected in the K_s-band spectrum: the H₂ 1-0 S(1) emission line. Molecular gas is a known signature of star-formation regions. We calculated the mass of warm molecular hydrogen, , using the following equation (Dale et al. 2005): where F_{ν = 1 → 0S(1)} is the integrated flux in the H₂ 1-0 S(1) line. We found M(H)   ~  80 M_⊙, thus, assuming a population fraction f_{ν = 1, J = 3} = 0.0122 (Scoville et al. 1982), a total H₂ mass 6500 M_⊙. This mass leads to a density in the nodules n(H₂) ≃ 200 cm^-3 assuming a cylindrical geometry. Thus 10% of H atoms are in a molecular state, which is compatible with the high extinction found Sect. 4.3.1.

Because only one H₂ emission line is detected, we cannot conclude on the excitation mechanism of molecular hydrogen. However, the jet of NGC 1068 with a velocity greater than 700 km s^-1 (Axon et al. 1998) is dissociative for H₂. The excitation may be caused by UV pumping from stars or by radiative UV shocks caused by the jet, or both.

4.3.3. Coronal lines

• Si vii (E_ion = 246  eV), already detected inNGC 1068 by Marconiet al. (1996) and Prietoet al. (2005) and[Si ix] (E_ion = 341  eV), also detected by Marconiet al. (1996) and Geballeet al. (2009) are the most intenselines in our averaged spectrum.

• Ca viii (E_ion = 147  eV) was also detected in NGC 1068 (see Marconi et al. 1996; Martins et al. 2010; Sánchez et al. 2009) and appears in our averaged spectrum.

• Al vi (E_ion = 190  eV) has been detected by Geballe et al. (2009). A small feature is visible at this wavelength in our averaged spectrum, but the poor signal-to-noise ratio does not allow us to conclude on its detection.

• Coronal lines indicate the presence of a strong ionizing field. Mechanisms could explain the presence and the intensity of these lines: photoionization and shocks ionization, fully described in Geballe et al. (2009) and Kraemer & Crenshaw (2000a).

• The source IR-7, located at 120 pc from the central source, is an extended cloud detected in the K_s-band averaged spectrum seen on 12.8 μm VLT/VISIR images (Galliano et al. 2005; Poncelet et al. 2007). We detected the [Si vii] and [Ca viii] emission lines, which appear blueshifted corresponding to an equivalent velocity of +600 km s^-1 with regard to the systemic of NGC 1068. If this structure is caused by interaction of the jet with the ISM and moves because it is swept along by the jet, then this velocity may be representative of the jet velocity.

4.4. CO 2-0 band head

We detected the CO 2-0 band head at the K_s-band spectrum with an equivalent width EW_nod = 1.14 × 10^-3 μm while in a region around the nodule EW_out = 4.22 × 10^-4 μm. The CO 3-1, CO 4-2 and CO 5-3 band head are not detected in the averaged nodule spectrum because of the fairly poor SNR. These absorption lines indicate the presence of stellar activity (late K- or early M-type stars). This feature has already been observed by Martins et al. (2010) at the central core of NGC 1068.

5. Origin of the coronal lines

We used Cloudy v08.00 (Ferland et al. 1998) to compare our data to a model of radiative transfer in the nodules and determine the ionizing mechanism responsible for the observed coronal lines.

Cloudy requires a geometrical and physical description of the nodules (chemical abundances, grains, distance, geometry, etc.) and the incident ionizing continuum.

We ran our simulations for a single average nodule at the location of IR-1. The adopted chemical abundances were taken from Shields & Oke (1975) and Ogle et al. (2003). Concerning dust, we used the silicate-graphite-PAH dust model proposed by Draine & Li (2001, 2007). The grain size distribution is based on the work of Laor & Draine (1993) and modified to include very small grains because stochastic heating of nano-grains has been proposed to explain the observed color temperature of the nodules at infrared wavelengths (Rouan et al. 2004; Gratadour et al. 2005). The distribution function of the grains size d is a simple power law: with β = 3.3. The minimal size is 1 nm and the maximal size is 1 μm. We assume that silicates and graphites have the same distribution function. Finally, PAHs are composed of a maximum of 100 atoms and their abundance is equivalent to that of the interstellar medium in the solar vicinity. We also assume a constant ionized gas density across the nodule.

We present below the results obtained with Cloudy for each ionizing continuum and compare them to our data to explain the presence of coronal lines in the nodule spectra.

5.1. Central source

The ionizing continuum of the central source of an AGN is fitted well by the sum of two components:  (Mathews & Ferland 1987). T_BB represents the temperature of the Big Blue Bump, α_UV is the slope of the ultraviolet continuum, α_X the one of the X continuum and a is set to fit the ratio α_OX between X rays and UV. kT_IR is taken as 0.01 Ryd. For our Cloudy simulations, kT_BB = 4 Ryd (T_BB = 6.3 × 10⁵  K) according to Spinoglio et al. (2005). Finally, α_OX =  −1.55, α_UV =  −0.93 and α_X =  −1 (Pier et al. 1994).

We used the ionization parameter , which is the ratio between the density of ionizing photons over the hydrogen density, to specify the intensity of the spectrum of the ionizing source. is the ionization rate of the central source. Spinoglio et al. (2005) found U = 0.1 at the edge of the NGC 1068 narrow line region (NLR) i.e. 2 pc from the central source. The density of photons decreases as r^-2, with r the distance to the source, then at the location of the cloud IR-1 (~30 pc), the ionization parameter is U ~ 4.4 × 10^-4 (log U =  −3.35). Using this value of U, the ionizing luminosity of the incident continuum generated by Cloudy is L_ion = 10^43.7   erg s^-1, a typical value for Seyfert galaxies (Pier et al. 1994; Geballe et al. 2009).

Fig. 10 [Si vii]/[Ca viii] flux ratio with respect to U for Cloudy simulations for the central source as the ionizing source. In green we depict the measured ratio for the averaged nodule, in red the measured ratio for IR-1b with the 3-σ confidence interval. For nodule IR-1 (i.e. 30 pc), U =  −3.35 and for IR-1b (i.e. 14 pc), U =  −2.69. For log U <  −2.75, the Cloudy precision is insufficient to provide confident values: the [Si vii] and [Ca viii] fluxes quickly tend to zero.

At IR-1, results from Cloudy simulation show that no coronal line can be produced by the central source of NGC 1068: observed line ratios are not reproduced by simulations (Fig. 10). Another argument is that IR-4 is twice as far away from the central source as IR-1. Hence, the ionizing flux from the central source is then four times weaker than at IR-1. However, the flux in highly ionized IR lines is almost constant with the distance to the central source (cf. Fig.  7). Our first conclusion is that a local ionizing source is necessary to explain this behavior.

It has been shown that dense matter-bounded clouds (Binette et al. 1996, 1997) can successfully produce high-ionization lines with nuclear photoionization and the extrapolation from the inner NLR of the drop in the density of photons as r^-2 is fraught with risks. However, in our case there is no obvious evidence of these clouds because up to 10 pc the spectra are dominated by the continuum of the central source (Fig. 4). Thus, a photon density decreasing as r^-2 is a reasonable assumption.

Close to the central source, at the location of IR-1b (log U =  −2.69), the influence of the central source is not negligible. Cloudy simulations produce [Si vii] and [Ca viii] with a line ratio equal to 2.36, while the observed ratio is equal to 7.40. Hence, 30% of the total ionizing field at this distance might be caused by the central source (assuming an equivalent composition for IR-1 to 4 and IR-1b). The rest is probably caused by radiative shocks caused in turn by the jet-ISM interaction. Indeed, Sánchez et al. (2009) provided evidence that a giant molecular cloud caused the jet direction to change. In Rouan et al. (2004), we also see a strong correlation between IR-1b and the the jet direction change.

Cloudy simulations show that coronal lines are produced by the central source with a high [Si vii]/[Ca viii] ratio up to 10 pc. Farther from the central source, additional ionizing sources are necessary to explain the observed spectra.

5.2. Super star clusters

Stellar activity is a known source of ionization. We demonstrated that a stellar continuum allows us to explain the shape of the observed nodule spectrum (cf. Sect. 4.2). We detect at the K_s-band a CO band head (2.3 μm) that is a K- or M-type star indicator. We also detect the Brγ and Brα lines, a young and hot OB star indicators. Let us explore the hypothesis of SSCs at the core of the nodules.

5.2.1. Characteristics of SSCs

(1)The expression of the star-formation rate is given by Eq. (1) (Kennicutt 1998), where L_Brγ is the luminosity of the Brγ line. Using the Brγ line flux, we found a SFR = (5.2 ± 0.19) × 10^-1   M_⊙ yr^-1, a typical value for galactic nuclei: in the Hα sample of Ho et al. (1997) the nuclear SFRs in most galaxies are quite modest, in average  ~0.2  M_⊙  yr^-1 (median 0.02  M_⊙ yr^-1). We found then, with Eq. (1), and thanks to this parameter, we are able to count stars in the super star cluster and calculate the total mass with the initial mass function.

Sternberg et al. (2003) gives the ionization rate for B0.5 to O4 stars for masses greater than 20  M_⊙ (most of ionizing photons comes from those stars). Assuming a Kroupa initial mass function (Kroupa 2001), we calculated the number of stars necessary to produce the observed Brγ line flux in the averaged nodule spectrum, following the method of Kornei & McCrady (2009). The best fit gives 2 821 ± 101 OB stars. However, this is an upper limit because we consider that all Brγ line flux comes from stars.

Blue giants are not the only high-mass stars in a cluster because evolved red giants are likely to contribute as well. A way to calculate the number of red stars is to use the observed CO 2-0 band head. Our aim was to reproduce the observed flux in the averaged nodule spectrum considering that only stars following the Kroupa IMF with mass greater than or equal to 20   M_⊙ are present (lower masses did not have time to evolve). No mass-luminosity relationship exists for red giants (Malkov 2007). Hence, we assumed the mass-luminosity relation L ∝   M^2.76 given by Vitrichenko et al. (2007) for massive stars and Betelgeuse as a reference for mass and luminosity (M = 15   M_⊙ and M_K = −4.38 mag). Concerning the CO band head, we assumed that the ratio is almost constant whatever the type of star. Using the ESO stellar library (Pickles 1998), we found . Finally, we calculated the number of stars necessary to produce the observed CO band head flux in the averaged nodule spectrum, and we found 237 ± 12 red giants.

Finally, from this estimated number of stars with a mass greater than 20  M_⊙, we constrained the IMF and we found an averaged total mass of stars of  ~8.07 × 10⁵  M_⊙ for masses varying between 0.08  M_⊙ and 90  M_⊙.

We deduced an estimate of the age of the stars using Eq. (1): τ ~ 1.55   Myr. This age corresponds to a 30 M_⊙ star period on the main sequence. Assuming a cylindrical geometry for nodules, we obtain a density of stars greater than 5000   stars pc^-3, a number compatible with typical SSCs which contain 10⁴  to  10⁶ stars for a diameter of 1  to  5   pc, so a density higher than 5 × 10²   stars  pc^-3 (Johnson 2005).

5.2.2. Cloudy simulation spectrum

Using the Tlusty stellar atmosphere model of Lanz & Hubeny (2003), we ran a Cloudy simulation to determine if an SSC is able to produce coronal lines. Apparently, a typical SSC cannot contribute significantly to the ionizing field producing the observed coronal because the extreme UV or X luminosity of OB stars is weak. Indeed, the ionization rate derived from equation 1 leads to U = (8.42 ± 1.24) × 10^-4 and thus log U = −3.07 ± 0.06 for the stellar ionizing field and this value of the ionization rate is insufficient to produce coronal lines (Fig. 10).

5.3. Jet

A jet is a source of shocks causing UV-X rays emission. The latter is a possible excitation source for high ionization potential. Allen et al. (2008) studied the continuum of high velocity radiative shocks for velocities up to 1000 km s^-1.

In this work we used different continua in our Cloudy simulations corresponding to radiative shocks for different jet velocities because the flux of the radiative shocks at high energies increases with the jet velocity (Fig. 2).

Fig. 11 Spectra from the MAPPINGS III library of fast radiative shocks (Allen et al. 2008). In green we depict the continuum produced by radiative shocks for a jet velocity of 700   km s-1 and in red for a jet velocity of 1000   km s-1.

For NGC 1068, Kraemer & Crenshaw (2000b) used the shock-induced radiation model generated by Wilson & Raymond (1999) for a shock velocity of 1000 km s^-1 and reproduced the NGC 1068 ultraviolet and optical spectra with a fairly good precision.

We specified the intensity of the radiative shocks using Eq. (2) from Allen et al. (2008), where I_UV is the ionizing intensity, v_shock is the shock velocity and n is the density of the medium. Beyond 1000 km s^-1, continua are not available and we assume that the shape of the continuum does not vary. This approximation will lead to a slightly over-estimated jet velocity because of the under-estimated flux at high energy. The model takes the magnetic field into account. But this has a negligible impact at energies higher than 100 eV (assuming a jet velocity of  ~1000 km s^-1 and a density of 1000 cm^-3) and we chose its value as close as possible to the measured value for NGC 1068 i.e. 1200 μG (Gallimore et al. 1996b). This model contains two emission components: the prompt emission from the shock and the emission from radiative shocks. We verified that the first component is several order of magnitude fainter than the second component. Hence, the prompt emission from the shock will be neglected. (2)Results from simulations are presented in Figs. 12 and 13.

Fig. 12 Results from Cloudy simulations considering the jet as the source of the ionizing field. The solid line represents the results of the Cloudy simulations. Dashed lines represent the observed line ratio with its 3-σ confidence interval.

Fig. 13 Idem Fig. 12 for the lines ratio [Si vii]/[Si ix].

According to Axon et al. (1998), results for the jet velocity lower than 700 km s^-1 are not valid for NGC 1068 and are not taken into account. Considering the shock models of Allen et al. (2008), the observed ratio [Si vii]/[Ca viii] is reproduced for a jet velocity higher than 1200 km s^-1. The ratio [Si vii]/[Si ix] is reproduced for a jet velocity of  ~1700 km s^-1.

Both results are consistent with Axon et al. (1998), who claimed that the jet velocity should be higher than 700 km s^-1.

The radiative shocks models is able to produce coronal lines. Hence, the jet is our best candidate for the ionization source responsible for coronal lines at the location of the nodules.

6. Conclusion

For the first time, thanks to the spatial resolution provided by NaCo, near-IR coronal emission can be precisely associated to structures observable in the broad-band images. Coronal emission is not distributed homogeneously in a narrow region around the central core but rather in dense clouds spanning a region as wide as 70 pc from the nucleus at high angular resolution.

The analysis of our long-slit spectroscopic observations of these regions allowed us to highlight the effect of the jet on the coronal lines and a significant star-forming activity in the closest regions of the core of NGC 1068 and to determine the nature/properties of the nodules and the probable excitation mechanism for coronal lines:

• The Brγ and Brα emission lines were used to determine the hydrogendensity in the nodules (~2 × 10³ cm^-3) and the extinction onthe line of sight (A_V ≃  17).

• The H₂ ro-vibrationnal 1-0 S(1) is detected in our spectra. We estimated the molecular hydrogen density to be 10% the atomic hydrogen density, which is compatible with the extinction above.

• The slope of the nodules continuum at the L′-band indicates a dust temperature of  ~500 K, which is consistent with previous studies. The peculiar slope at the K_s-band is only reproduced by introducing a second continuum with a temperature of 7500 K, consistent with a temperature of OB stars reddened by dust and gas.

• The H i Brγ emission line and the CO 2-0 band head also indicate the presence of stars. Inside the nodule, a density of about 5000 pc^-3 was found, compatible with that of super star clusters.

• The ionizing field of these stars is unable to produce coronal lines but is a possible excitation mechanism for H₂ 1-0 S(1).

• Beyond 10 pc from the central source, the ionization field of the central source is not sufficient to produce coronal lines: its effect on the nodules (at 30–60 pc from the central source) is negligible.

• Ratios between coronal lines intensities are almost constant with the distance from the central source, which is another clue for a local ionizing source of excitation.

• Results from the Cloudy photoionization simulations allow us to conclude that the jet of NGC 1068 is the probable source of excitation for coronal lines (via radiative shocks).

• Using the shock models of Allen et al. (2008), we show that the line intensity ratios can be reproduced for a jet velocity of about 1700 km s^-1.

• The emission lines at the location of source IR-7 exhibit a blueshift. Considering this structure as jet-driven, we estimated the jet velocity to be higher than 600 km s^-1 in agreement with previous studies and this work.

New observations with better atmospheric conditions will be able to considerably improve the present work. Moreover, observations with a Fabry-Perot imager at a resolution similar to the one achieved in this work around the [Si vii] coronal line could provide more information on the jet-cloud interaction process in the core of this archetypal AGN.

Acknowledgments

The authors thank the anonymous referee for the comments, which have improved the quality of this paper. J. Exposito is very grateful to Jacques Le Bourlot, Didier Pelat and Guillaume Drouart for their help in the molecular hydrogen line interpretation, for the estimation of errors on measurements, and for the help on the stellar contribution in our spectra.

References

All Tables

All Figures

	Fig. 1 Radio jet contour map (left) from Gallimore et al. (1996b), [O iii] emission (center) from Macchetto et al. (1994) and mid-IR contour (right) from (Bock et al. 2000) over the deconvolved M′-band image (Gratadour et al. 2006).
In the text

	Fig. 2 Nodule identification and slit orientation on M′-band image (left) and its deconvolved image (right). The distance between the ticks inside the slit is 0.1″ i.e. about 7.5 pc. In white, the 47° PA slit aligned with the nodules. In green, the 12° PA slit. The colors are in arbitrary units and are chosen so that the structures mentioned in this article appear clearly.
In the text

	Fig. 3 Black: spectra from DS1 for four different positions in the slit. Red: fits by one blackbody with temperatures summarized in Table 1.
In the text

	Fig. 4 Ks-band spectra extracted from various regions along the slit with identified IR emission lines. These spectra were renormalized to be displayed on the same plot. The line intensity variation is not absolute.
In the text

	Fig. 5 Spectrum from DS2 at the Ks-band. The spectrum is the sum of the four nodule spectra.
In the text

	Fig. 6 Idem Fig. 5 at the L′-band.
In the text

	Fig. 7 [Si vii]/Brγ flux ratio as a function of the distance from the central source in blue, superimposed on the maximum intensity along the spatial dimension of the slit as a function of the distance from the central source in black (arbitrary units). The ratio is almost constant whatever the distance. In green, the position of the nodules identified.
In the text

	Fig. 8 Spectra of the nodules at the Ks-band in arbitrary units so we can display them in the same plot.
In the text

	Fig. 9 Fit of the spectra with a unique component at 500 K (green), and with a sum of a graybody and a blackbody (red).
In the text

Fig. 10 [Si vii]/[Ca viii] flux ratio with respect to U for Cloudy simulations for the central source as the ionizing source. In green we depict the measured ratio for the averaged nodule, in red the measured ratio for IR-1b with the 3-σ confidence interval. For nodule IR-1 (i.e. 30 pc), U =  −3.35 and for IR-1b (i.e. 14 pc), U =  −2.69. For log U <  −2.75, the Cloudy precision is insufficient to provide confident values: the [Si vii] and [Ca viii] fluxes quickly tend to zero.

In the text

	Fig. 11 Spectra from the MAPPINGS III library of fast radiative shocks (Allen et al. 2008). In green we depict the continuum produced by radiative shocks for a jet velocity of 700   km s-1 and in red for a jet velocity of 1000   km s-1.
In the text

	Fig. 12 Results from Cloudy simulations considering the jet as the source of the ionizing field. The solid line represents the results of the Cloudy simulations. Dashed lines represent the observed line ratio with its 3-σ confidence interval.
In the text

	Fig. 13 Idem Fig. 12 for the lines ratio [Si vii]/[Si ix].
In the text