© The Authors 2023

1. Introduction

The central region of NGC 1068 provides a unique opportunity to study some of the smallest spatial scales (1 arcsec ↔ 72 pc) in an active galactic nucleus (AGN) due to its proximity (14.4 Mpc) (Bland-Hawthorn et al. 1997). The galaxy’s nearly face-on orientation offers a clear line of sight toward its central region, and its location near the celestial equator allows for observation from both hemispheres. As a result, NGC 1068 has been the subject of numerous publications since its original spectroscopic observation by Karl Seyfert in 1943 (Seyfert 1943).

Observations of NGC 1068 at subarcsecond resolution from UV to radio wavelengths have been made possible by the use of space-based imaging, adaptive optics, and interferometry. These high-resolution observations have provided important insights into the physical processes occurring in the central region of this galaxy.

At gigahertz frequencies, a kiloparsec-scale (∼10″) radio jet oriented along a position angle of approximately 30 degrees is clearly visible in NGC 1068 (Gallimore et al. 1996; Muxlow et al. 1996). At smaller spatial scales (Gallimore et al. 2004), it appears that the jet is initially launched at a position angle of approximately 12 degrees from a compact radio source named S1, but it is bent approximately 0.3″ north of it at the position of another compact radio source named component C (see Fig. 1). This deviation of the jet is proposed to be due to an interaction with an interstellar cloud (Gallimore et al. 2004).

Fig. 1. Top: schematic representation of the main components in the central region of NGC 1068: the CND, the outer edges of the outflow bicone, and the three regions (B, C, and D) analyzed in this work. The field of view of the SINFONI observation is indicated with a pink square. Bottom: [Si VI] emission line map in the central region of NGC 1068 obtained with the SINFONI observation (pixel scale 50 mas) and position of the SPHERE slit (red lines, roughly aligned with the jet orientation). The coordinates are centered on the position of the K band continuum maximum. Following the naming convention from Evans et al. (1991), the three regions studied in this paper are cloud B (in the center; at the same position as the radio component S1 and the compact torus), cloud C (0.25″ north), and cloud D (0.7″ north).

At longer wavelengths, the molecular content of NGC 1068 is detected through CO emission lines, notably with ALMA (García-Burillo et al. 2019). The main feature of the central region is a massive circumnuclear molecular disk (CND) that is 400 pc wide and has a mass of 10⁸ M_⊙. The CND has a central hole of approximately 200 pc in diameter, within which the AGN is located, but not at the center. A smaller molecular disk with a mass of 3 × 10⁵ M_⊙ has been observed at the position of S1, which has been identified as the mass reservoir and obscures the structure around the AGN (García-Burillo et al. 2016, 2019; Imanishi et al. 2018). A molecular clump has also been observed north of S1, roughly at the bending point of the radio jet component C. High-resolution continuum observations in the ALMA spectral range have revealed a similar geometry to the observations made at gigahertz frequencies (García-Burillo et al. 2019; Imanishi et al. 2020).

In the IR, the details of NGC 1068 become accessible with adaptive optics-fed instruments. Gratadour et al. (2006) used deconvolution techniques on NAOS-CONICA observations to produce high angular resolution images of the central arcsecond in the K, L, and M bands. The nucleus S1 is the brightest source, particularly at short wavelengths. However, several compact sources with significant flux can be observed around it, particularly in the M band, where the northern sources are almost as bright as the nucleus. The positions of the two brightest sources, that is, the nucleus and IR1B, correspond to the positions of the radio components S1 and C, respectively.

At the shortest wavelengths, the central region of NGC 1068 has been imaged in the optical and UV wavelengths with the Hubble Space Telescope Faint Object Camera (Macchetto et al. 1994). The images of the source in the optical and UV continuum, as well as with the [O III] emission line filter, are very similar to each other on a large scale, revealing an extended ∼400 pc structure oriented along a position angle of approximately 30 degrees. The [O III] image revealed the complex filamentary structure of the narrow line region (NLR), which is an outflowing ionized bicone. Several bright emission line clouds were observed in the central arcsecond (Crenshaw & Kraemer 2000; Das et al. 2006).

In this paper, we present a detailed study of the spatially resolved near-IR (NIR) spectra of three clouds (B, C, and D) from the inner NLR of NGC 1068 (named following the original convention from Evans et al. 1991). In Sect. 2.1, we present the observations, followed by a quantitative description of the spectra in Sect. 3. In Sect. 4, we present the associated modelings, while a discussion of the nature of the objects is presented in Sect. 5. Finally, in Sect. 6, we summarize our conclusions.

2. Observations and spectra extraction

2.1. Observations

This study is based on the merged analysis of two spectroscopic observations obtained with SPHERE/VLT and SINFONI/VLT. These observations are discussed individually in the following paragraphs.

2.1.1. SPHERE

We used a new observation performed with SPHERE in the IRDIS long-slit spectroscopy (LSS) mode (Vigan et al. 2008; Beuzit et al. 2019) as part of the 0104.B-0242(A) ESO observing program (PI: J.-L. Beuzit). The observation was conducted under changing atmospheric conditions. The seeing was around 0.7″ at the beginning of the observation (with a coherence time of the atmospheric turbulence between 5 and 6 ms) but quickly degraded to more than 1″ (with a coherence time below 4 ms). As a result, the adaptive optics’ performance varied a lot during the observation, providing a good resolution for the first few frames and quickly worsening. As a compromise between SNR and angular resolution, we selected the three best exposures of the object in terms of angular resolution and discarded the others.

The data reduction followed the instructions from the SPHERE manual¹ and was performed as described in detail in Vermot et al. (2019). We applied detector-level corrections, including dark subtraction, flat-field division, distortion correction, linear interpolation of the bad pixels, and realignment of the various exposures. As a second step, we performed spectroscopic calibrations, correction for the atmospheric transmission, and flux calibration from the 2MASS J and H magnitudes of BD-00413, the standard star observed during the run.

The final product is a flux-calibrated spectrum with a 100 mas angular resolution and an 8 nm spectral resolution. The slit covers 11″ × 0.09″, is centered on the maximum of continuum emission, and is oriented with PA = 11° (see Fig. 1).

2.1.2. SINFONI

We used archival SINFONI integral field spectroscopy obtained as part of the observing program 076.B-0098(A). The instrumental setup provides spectrally dispersed images in the H and K bands with R ∼ 1500 and r = 100 mas.

The data reduction was performed with the ESO pipeline in its standard mode, except for the flux calibration, which was performed a posteriori by matching the integrated HK band flux from the entire field of view to the one measured on 2MASS images of the galaxy. The image presented in Fig. 1 corresponds to the [Si VI] continuum-subtracted emission line map obtained from this observation.

2.2. Data interpolation and spectrum extraction

To combine the two observations, we computed a synthetic long slit on the SINFONI data cube. We used linear interpolation to rotate the cube in the spatial dimension around the maximum of continuum emission and extracted a two pixel-wide slit matching the orientation of the SPHERE one. Since this synthetic slit is 100 mas wide while the physical one from SPHERE is 90 mas, we applied a 0.9 correction factor to the SINFONI flux.

We extracted three spectra from each slit at positions corresponding to the maxima of emission lines: 0.0″ (cloud B), 0.25″ (cloud C), and 0.7″ (cloud D). Each spectrum was summed over a 100 mas bin. The Y and J bands are only covered by SPHERE, the K band only by SINFONI, and the H band by both observations. For the latter, we interpolated the SINFONI data to the spectral resolution of SPHERE and took the average between the two observations. The match was perfect between the two datasets for clouds C and D, but the SINFONI spectrum of cloud B required recalibration in order to increase its flux by a factor two, which we attribute to a combination of imprecision in the synthetic slit extraction and differences between the PSFs of the two adaptive optics systems. Lastly, the continuum of emission was subtracted from the spectrum of each cloud to extract a pure emission line spectrum.

3. Detected emission lines

The continuum-subtracted spectra of the three clouds are presented in Fig. 2, color coded by spectral domain. The green and blue parts correspond, respectively, to the SPHERE and SINFONI exclusive spectral zones, while the red part corresponds to the overlap of the two instruments. The spectra are presented, from top to bottom, with increasing distance from the nucleus. The uncertainty on the spectrum, indicated as a gray shaded area, was estimated by measuring the standard deviation of the spectra in regions free of emission lines. As can be observed, the SNR also increases with the distance from the center despite the flux of the emission lines decreasing. This is because the photon noise from the continuum of emission is at its maximum at the center.

Fig. 2. Combined YJHK spectra of the clouds with increasing distance from the nucleus. From top to bottom, the panels show a distance of 0.0″ (B/S1/nucleus), 0.25″ (C), and 0.7″ (D). The SPHERE data are displayed in green, the SINFONI data in blue, and the averaged overlap in red. The estimated uncertainty is indicated as a gray area.

A significant result can already be deduced from the qualitative analysis of these spectra: they are remarkably similar to each other. They are dominated by a very strong emission line at ∼1.10 μm, to the left of which are three medium emission lines. Redward, two medium emission lines are found at ∼1.25 μm and one at ∼1.45 μm, and lastly, another strong emission line is found at ∼1.95 μm. This is valid for the spectra of clouds C and D, which are almost indistinguishable apart from their fluxes. Part of the spectrum of cloud B is unusable, but the emission lines in the J band are also fully consistent with the description above and with clouds C and D. This indicates that the three emission line regions are excited by the same mechanism, and most probably, the properties of the interstellar medium (ISM), for example, ionizing radiation field, density, and composition, are very close in these three clouds.

Altogether, we detected a few dozen emission lines, spanning a variety of elements and ionization energies. We present in Table 1 the flux of all the emission lines detected with SNR > 3. We measured them with Gaussian fits (simple in the case of isolated lines and double in the case of overlapping lines). The uncertainties provided correspond to 1σ, as estimated from the covariance matrix of the parameters. As discussed later in this paper, the identification of the lines has been done a posteriori by analyzing the output from the best CLOUDY models. We detected permitted emission lines, notably H I and He I. The He I line is the strongest emission line of the spectrum, and we detected up to seven emission lines from different hydrogen series. Assuming the Case B recombination case (clouds opaque to Lyman UV photons), we measured no extinction (A_K ≤ 0.15 for all three clouds at 3σ). We also detected forbidden low-ionization energy emission lines, indicating the presence of a low-density environment. They are numerous and have relatively low fluxes. Lastly, we detected forbidden high-ionization energy emission lines, also known as “coronal emission lines”. Two of them are of particular importance: [Si X] at 1430 nm, which has the strongest ionization energy among our lines (401 eV), and [Si VI] at 1960 nm, which is the second strongest line in our spectral domain.

4. Physical conditions and ionization mechanism

4.1. Goal and method

In this section, we investigate the physical conditions and ionization mechanism of the clouds. To achieve this, we compare the observations with three ionization models: (1) a simple CLOUDY c17.02 (Ferland et al. 2017) setup describing a cloud with constant density and temperature (without any ionization source; models PC1, PC2, and PC3); (2) a CLOUDY model of a cloud photoionized by a central AGN (model AGN); and (3) the MAPPINGS III shock model (which itself relies on CLOUDY; model SHK).

Instead of using a specific emission line ratio as a diagnostic between the different situations, we used the models to compute synthetic spectra with the same resolution as the observation and compared them directly to the observed spectra. This approach has several advantages over the use of emission line ratios. First, it takes into account every emission line. Second, it takes into account all the non-detections, which bring important constraints on the model. Third, it takes into account the absolute flux of the emission lines. Fourth, it makes no a priori assumptions about the identification of the lines. Moreover, overlapping close lines are naturally taken into account. Finally, by studying the best model, we can make an identification of the lines a posteriori.

We removed the [Fe II] lines from the models because they should be strong and ubiquitous according to our ionization models, but they are barely detected. The explanation for this weak [Fe II] signal is discussed in Sect. 5.2.3. Otherwise, synthetic spectra were computed and directly compared to the observations. For each model, we determined the best set of parameters by comparing the reduced chi-square between the model and the observation.

4.2. Physical conditions: Constant density and temperature

To predict the emission line spectrum of a cloud with a uniform density and temperature on a large and fine grid, we first determined the physical conditions of the emitting regions by using CLOUDY c17.02. We call this model PC1. The grid covers a wide range of temperatures (log(T/K) ranging from 3.70 to 9.00 with 0.02 steps) and densities (log(n/cm⁻³) ranging from −5 to 9 with 0.05 steps. We compared the spectra in absolute units of flux, assuming an emitting volume smaller than 125 pc³.

The results for all three clouds are quite similar, with a deep minimum observed in the residual maps for temperatures ranging from 10^4.2 K to 10^5.2 K and densities of 100 cm⁻³ or greater. The corresponding 2D maps for these parameters can be found in Fig. B.1, while the 1D plots are presented in Figs. B.3 and B.2. The lower limit on density (n ≥ 10^0.5 cm⁻³) is imposed by the maximum volume of the emitting cloud, but otherwise this parameter has little effect on the spectrum.

The best models shown in the top plots of Figs. A.1–A.3 successfully reproduce the permitted emission lines (He I and hydrogen series) but do not produce any forbidden emission lines. The best-fit parameters for these models, including the corresponding emitting volumes, are provided in Table 2. The optimal temperature values for all three clouds are quite similar, ranging from 10^4.8 to 10^4.9 K.

As indicated by Figs. B.1 and B.2, a less likely secondary solution is also present in all three clouds, characterized by temperatures ranging from 10^5.5 to 10^6.5 K. These solutions predict strong forbidden emission lines but weak permitted emission lines.

In the next step, we attempted to find a two-temperature combination that can reproduce the observed spectrum. To this end, we fixed the density at n = 100 cm⁻³ and computed the spectra for all possible temperature combinations, with flux ratios between the two components ranging from 0.05 to 0.45 in increments of 0.05. The results from this approach, which we call model PC2, are detailed in Table 3, with residual maps in Fig. B.4. For all three clouds, the majority of the emission is still attributable to a medium at ∼10^4.8 K (with acceptable temperatures ranging from 10^4.5 to 10^5.2 K). Clouds B and C exhibit a secondary medium with temperature 10⁶ K (with corresponding flux ratios of 15% and 25%), while cloud D has a secondary medium with a temperature of 10^5.4 K.

A comparison of the best PC2 models for clouds B, C, and D is presented in Figs. A.1–A.3. The results show that the best models for clouds B and C are successful in reproducing the permitted emission lines as well as most of the forbidden emission lines, but they fail to generate significant flux for the strong [Si VI] coronal emission line at 1960 nm. In contrast, the best model for cloud D accurately reproduces the He, H, and [Si VI] emission lines, but it fails to generate the correct flux for several other forbidden emission lines, specifically the coronal lines [S IX] and [Si X] at 1251 and 1430 nm. The best solution for clouds B and C is a good secondary solution for cloud D, and vice versa.

Each of the three clouds exhibits a multiphase ionized medium, with temperatures ranging from 10^4.8 to 10⁶ K. We created a third model (PC3) with three emitting media; the density fixed at n = 100 cm⁻³; and the temperatures set to T₁ = 10^4.8 K, T₂ = 10^5.4 K, and T₃ = 10⁶ K. We then determined the relative contribution of each medium to the observed spectrum. The best-fit parameters are presented in Table 4, and the corresponding spectra are shown in Figs. A.1–A.3. This model produces the best results, accurately predicting the flux of both permitted emission lines and the most important forbidden lines ([S IX], [Si X], and [Si VI]) for all three clouds. The ionized phase’s mass is estimated to be approximately 1 M_⊙, 0.75 M_⊙, and 0.5 M_⊙ for clouds B, C, and D, respectively. In each case, the primary phase dominates at T = 10^4.8 K (contributing to around 65% of the flux), followed by T ∼ 10^5.4 K (contributing to approximately 20% of the flux) and T ∼ 10⁶ K (accounting for roughly 15% of the flux).

4.3. Central photoionization

In a second model, we attempted to reproduce the observed spectra with a photoionization model of a cloud around a central radiation source. For this purpose, we used CLOUDY to simulate a cloud illuminated by an AGN. The luminosity of the source was relatively well constrained and set to 10³⁸ W (Bland-Hawthorn et al. 1997; Gallimore et al. 2001; Vermot et al. 2021). The cloud’s thickness was set to 10 pc, and its distance from the nucleus was 30 pc for clouds C and D, whereas for cloud B, it was set to 1 pc. In our simulations, we varied three parameters: the density of the cloud, ranging from 10⁻² to 10⁹ cm⁻³ in 10^0.2 logarithmic steps; the temperature of the big blue bump (BBB), logarithmically ranging from 10⁴ to 10^6.5 K in logarithmic steps of 10^0.1 K (we note that the luminosity of ionizing photons is always 10³⁸ W); the optical to X-ray spectral index α_OX, as defined in Zamorani et al. (1981) and implemented in CLOUDY with an explicit negative sign. Typical AGN have α_OX = −1.4, while α_OX = 0 corresponds to the absence of X-ray emission.

With this model too, the residual maps are very similar from one cloud to another. An important degeneracy is observed between the temperature of the BBB and the α_OX parameter: the higher the temperature of the BBB, the lower the amount of required X-rays (see Fig. C.1). In all cases, the model converges toward relatively low amounts of X-rays (α_OX ≥ −1.2). For all clouds, density is reliably constrained between 10^2.5 and 10^6.5 cm⁻³. For clouds B and C, the model converges toward a density n ∼ 10⁴ cm⁻³, while for cloud C, it converges toward n ∼ 10^5.5 cm⁻³ (see Table 5 for the exact values of the best parameters).

Comparisons between the model and observation spectra for the three clouds are presented in Figs. A.1–A.3. For cloud D, the model correctly predicts the flux of the permitted emission lines as well as the strong [Si VI] coronal line, but it fails to reproduce the other weaker forbidden lines. For clouds B and C, the model predicts the presence of most of them, but it overestimates the flux of the [Ca VIII] emission line in the K band and the hydrogen series, and it underestimates the [Si VI] line at 1960 nm.

4.4. Fast radiative shock

In the final modeling attempt, the observed spectrum was modeled as emerging from a region ionized by a shock using the MAPPINGS III Library of Fast Radiative Shock Models (Allen et al. 2008). The model assumes that the shock front generates a strong radiation field of extreme ultraviolet and soft X-ray photons, leading to significant photoionization in the shock itself and in the region ahead of it. The model varies three parameters: the velocity of the shock, ranging from 100 to 1000 km s⁻¹ in steps of 25 km s⁻¹; the preshock density, ranging from 10⁻² to 10³ cm⁻³ in logarithmic steps of 10¹; the preshock transverse magnetic field B, with values ranging from 10⁻⁴ to 10³ μG, in an irregular sampling covering the extremes expected in the ISM while also sampling more finely the magnetic field values near equipartition.

The library does not explore a full regular parameter space, as the sampling of the magnetic field (number of samples and selected values of B) differs from one value of the density to another. Consequently, 2D maps of the chi-square as a function of the parameters are sparse and difficult to read, and thus we only present1D projections (minimum value of the chi-square for a given parameter). These are available in Figs. D.1–D.3. They favor fast shocks (v ≥ 750 km s⁻¹) and exclude the lowest densities (n ≥ 10⁻¹ cm⁻³) and magnetic fields (B ≥ 10⁻¹ μG).

As with the previous models, the optimal parameters of the three clouds are very similar to each other. The minimum value of the chi-square is obtained for the highest velocity (v_s ∼ 1000 km s⁻¹), low densities (n ∼ 1 − 100 cm⁻³), and intermediate magnetic fields (B ∼ 1 − 100 μG). The exact best parameters for each cloud are presented in Table 6.

For the three clouds, the best model is in excellent agreement with the observation. All the main emission lines are correctly predicted by the model (permitted and forbidden) with fluxes very close to what is observed (comparisons between the spectra are presented in the bottom plots of Figs A.1–A.3). The ratio between the strong He I and [Si VI]₁₉₆₀ emission lines is well reproduced, as well as the flux of the other coronal lines, [S IX] and [Si X]. The flux of the hydrogen emission lines is overpredicted, but the shape of the complex emission features redward of 1100 nm is remarkably well reproduced.

The best model corresponds to a shock at 1000 km s⁻¹, which is the highest velocity of the MAPPINGS model. As a consequence, it is possible that higher shock velocities could provide an even better fit to the data.

5. Discussion

5.1. Similarities between the clouds

One of the most striking results of the analysis presented in this paper is the strong similarity between the emission line spectra of the three objects, which leads to an almost identical modeling. Not only are the ratios of the various emission lines extremely similar between them, but their absolute fluxes are also of the same magnitude. We conclude that they belong to the same class of object, with similar physical conditions created by the same ionization mechanism.

5.2. Ionization mechanism

5.2.1. Comparison between the models

It is clear from the comparison to the models that the emitting medium is multiphase, with temperatures ranging from T₁ = 10^4.8 K to T₃ = 10⁶ K. The detection of strong forbidden lines with ionization energies as high as 400 eV imposes the highest temperatures, but a lower energy phase must be present to explain the strong emission from non-fully ionized elements.

Photoionization models of a central source with a simple blackbody spectral energy distribution fail to reproduce the observed spectra, as they can only reproduce one of the above-mentioned phases, depending on the temperature of the central source. More complex photoionization models with a source combining UV and X-ray contributions can produce a multiphase medium when illuminating dense clouds, reproducing the observed He I flux as well as some forbidden lines. Nevertheless, this AGN central photoionization model fails to reproduce some features of the observed spectrum, such as the complex structure blueward of He I or the ratios between the many high-energy forbidden lines.

The last model presented in this paper, which assumes that the emission lines are stimulated by a fast radiative shock, provides a remarkably good visual fit to the data and can reproduce all the essential features of the spectra (see Fig. 3 and Tables E.1–E.3). Despite significant differences between the measured emission fluxes and those predicted by our model, the reduced chi-square value we presented remains relatively low, which can appear surprising. This is because it is calculated across all pixels in the spectrum, incorporating the many non-detected emission lines in the process, which effectively lowers the value of the reduced chi-square. Furthermore, we employed a meticulous Gaussian fitting method to measure the fluxes reported in Tables E.1–E.3 that yields more precise flux values than a simple sum of the corresponding pixels.

Fig. 3. Cloud D: comparison between the observed spectrum (top, black) and the best model (bottom, red). The numerical values of the fluxes from the model’s main emission lines are compared with those from the observation in Table E.3 from the appendix.

Given the presence of an AGN in the vicinity of the emitting media, it is reasonable to assume that fast shocks could be responsible for the excitation of the clouds. Therefore, we conclude that this model is the most plausible one to explain the observed emission lines. For clouds C and D, the chi-square analysis significantly favors this shock model over the central photoionization one (see Table 7). For cloud B, the two models have a similar chi-square since the K band spectra is too noisy to discirminate between the two.

Our best fit analysis reveals that the shocks must be at high velocity (close to or exceeding 1000 km s⁻¹); the preshock density is relatively high, in the range of 1 − 10 cm⁻³; and the magnetic field is near equipartition. These parameters provide insights into the structure of the medium, which is composed of three regions: the shock itself, which heats the gas up to temperatures of 10⁷ K and generates strong ionizing radiation; a preshock region that is ionized by this radiation and forms an H II region; and a cooling recombination region located behind the shock. For a quantitative description of the ionization, temperature, and density profiles around a shock with similar parameters (v_s = 1000 km s⁻¹, n = 1 cm⁻³, B = 3.23 μG), we refer to Figs. 4 and 8 in Allen et al. (2008).

5.2.2. Improved diagnostic diagram

“Baldwin, Phillips and Terlevich” (BPT) diagrams (Baldwin et al. 1981; Veilleux & Osterbrock 1987; Kewley et al. 2001) are a powerful tool using emission line ratios to distinguish between an AGN and star formation activity from optical spectra. In a series of recent publications, D’Agostino et al. (2019a,b) proposed an improved 3D version of this type of diagram that takes into account the measured velocity dispersion derived from the lines as well in order to distinguish between shocks, AGN activity, and star formation and applied it to NGC 1068. The observation, obtained with the Wide Field Spectrograph (WiFeS; Dopita et al. 2007, 2010), is at a lower spatial resolution than ours and does not allow the resolution of individual clouds in the central arcsecond to be obtained. However, their results appear consistent with ours, highlighting that the entire central region (400 pc <  = > 6″) is dominated by shocks (Fig. 9 in D’Agostino et al. 2019b).

In our work, the conclusion is purely based on the analysis of the fluxes of various emission lines, while in their work, the decisive argument is the measured velocity dispersion of the emission lines; when high, this is a strong indicator of shock excitation (Rich et al. 2011; Ho et al. 2014). We cannot directly apply their method to our observation since we are observing other emission lines and do not have enough SPHERE spaxels to compute a proper emission line ratio function. However, despite our lower spectral resolution, we can attempt to measure the velocity dispersion from the emission lines in the SINFONI dataset for clouds C and D.

The resolution of our observation is R ∼ 1500, corresponding to an expected full width half maximum (FWHM) of 1.35 nm at 2 μm for an unresolved emission line. For our strongest emission line, [Si VI] at 1965 nm, we measured an FWHM of 4 nm for cloud C and 3 nm for cloud D, corresponding to intrinsic velocity dispersions of 245 km s⁻¹ and 150 km s⁻¹, respectively. For [Ca VIII], we found similar values: 220 km s⁻¹ and 130 km s⁻¹. These values are not as high as the highest ones reported in D’Agostino et al. (2019a,b) for the central region, but they still correspond to the beginning of the shock-dominated sequence from the 3D diagram, especially if we consider that a) a higher velocity dispersion component could be blended in our line and b) spatial mixing could be present in their lower spatial resolution IFU observation. Overall, we consider that our results are fully consistent with theirs, strengthening the interpretation that the excitation of the emission lines in the central region of NGC 1068 is dominated by shocks.

5.2.3. [Fe II] emission lines

We noted an unexpected lack of [Fe II] line emission in our spectra. According to the photoionization and shock models used to describe the other emission lines, these lines should be prominent. This discrepancy led us to remove [Fe II] lines from our analysis. In this section, we explore why [Fe II] might be so weak.

In the ISM, UV and X-ray observations often find less iron than expected. Most of it, about 90%, is thought to be hidden in dust particles (Dwek 2016; Cia 2018). One might initially guess that iron is similarly depleted in our case, thereby accounting for the weakness of [Fe II] emission. However, our best explanation for the emission lines we do see involves a fast-moving shockwave traveling at 1000 km s⁻¹. Such a shockwave would destroy dust, freeing the iron. Moreover, if iron were depleted, we should also miss signatures of other elements found in dust, such as silicon, but we observed strong silicon lines. A closer look revealed that the [Si X] and [Si VI] lines we observed come from ionization states requiring extremely high temperatures (at least 10⁵ K and 10⁶ K, respectively), whereas [Fe II] lines are typically found at temperatures of ∼10⁴ K or less, as shown in Figs. F.1–F.3.

In the MAPPINGS shock model, there are two places where temperatures are below 10⁴ K: the gas that has not been hit by the shockwave yet and the cooling zone behind the shock (see Fig. 5 in Allen et al. 2008). The model predicts that the vast majority of the [Fe II] emission should arise in the cooling zone, due to the higher n_H density. So we need to explain why there is no [Fe II] there. We think there are two possibilities.

One possibility is that dust forms again after the shock. If dust particles form quickly in the cooling region, they could capture a lot of the iron, explaining why we do not see [Fe II]. The same would happen to silicon, but the [Si X] and [Si VI] lines would still be visible, as they come from the hot, shocked region. Dust formation might indeed be possible here, as there are enough atomic elements to form dust (e.g., Fe, Si, C, which are released from the destroyed dust); the shock would make the medium dense; and the dust wind from the central black hole could provide “seeds” for dust formation. We note that [C I], which is a low-ionization energy state coming from an element constitutive of dust (similar to [Fe II]; see Fig. F.4), is also overpredicted by the shock model (see Tables E.1–E.3), which supports this hypothesis. As a possible follow-up observation, the nondetection of the [Si II] line at 34.85 μm could confirm it.

The second possibility is that no low-temperature area exists after the shock. Due to intense radiation from the AGN central engine, it is possible that the cooling region behind the shock never cools to below 10⁴ K. Additionally, the central region of NGC 1068 is known to have recently experienced star formation (Storchi-Bergmann et al. 2012; Vermot et al. 2019; Rouan et al. 2019), so hot stars in the nuclear star cluster (or directly inside the clouds) could contribute to the heating of the ISM in the post-shock region, for instance, in the Orion H II region, which exhibits low [Fe II] emission (Walmsley et al. 2000). In that case, we would not see any [Fe II] because iron cannot emit it in those conditions.

To confidently choose between these two options, we would need more complex simulations. However, this is not within the scope of our current work.

5.3. Nature of the objects

5.3.1. Other properties

In addition to the emission line spectra, the objects have been detected in the mid-IR (Gratadour et al. 2006), forming a structure very similar to the one observed in the [Si VI] image presented in Fig. 1. Cloud B corresponds to the nucleus, cloud C to IR-1B, and cloud D to the superposition of IR-3 and IR-4, as named in Gratadour et al. (2006). The other clouds detected in the mid-IR also have counterparts in the [Si VI] image. While the secondary clouds are much less luminous than the nucleus in the K band, this difference in flux becomes negligible at longer wavelengths, such as the M band, where most of the energy is radiated: M_nucleus = 6.6, M_{IR − 1B} = 6.7, M_{IR − 3 + IR − 4} = 7.6. Therefore, all three objects can be considered strong sources of IR radiation.

The three clouds also have molecular counterparts. In Müller Sánchez et al. (2009), the authors identified several molecular structures in SINFONI observations through an H₂ rovibrational line. Clouds B and C correspond to the structures named southern and northern tongues, respectively, identified in the H₂ observations. The southern tongue is on the line of sight of the IR peak, and cloud D is located within a larger CO molecular structure, which is an overdensity of the CND. The dynamical modeling of the southern and northern tongues indicates that they could be molecular streamers fueling the nucleus. The mass of these streamers ranges from a few 10⁶ to several 10⁷ M_⊙, estimated by various methods. The CO counterpart to component B is the molecular disk observed at the position of S1, with a mass estimated at 3 × 10⁵ M_⊙. A small CO clump is observed north of it at the position of cloud C, and cloud D lies again within the larger CND. In our [Si VI] observation of the objects, a velocity gradient can be measured along clouds C and D (see Appendix G; cloud B is too noisy). The amplitude of the velocity gradient is similar between the two clouds, and if interpreted as gravitational rotation, it would correspond to masses greater than a few 10⁶ M_⊙.

The three clouds may also correspond to the radio sources detected in Gallimore et al. (2004) between 1 and 10 GHz. However, aligning the radio image with the NIR emission is uncertain due to the high uncertainty in absolute astrometry (Capetti et al. 1997) and the dissimilarities in morphology between the two images that prevent cross-registration. Despite this, if S1 and the IR emission peak are associated with the nucleus (Gallimore et al. 2004; Gratadour et al. 2006; Gámez Rosas et al. 2022), cloud B may correspond to the radio component S1 and cloud C to the radio component C. Cloud D could be associated with component NE, but with an offset of about 0.05″. At these wavelengths, the similarities between the objects are striking, particularly between S1 and C, which have comparable fluxes, spectral slopes, and geometries (see Figs. 3–5 in Gallimore et al. 2004). This comparison also holds for component NE at 5 and 8.4 GHz, but not at 1.4 GHz, where it appears more extended and luminous.

Thus, the similarities observed in the NIR emission line spectra of the clouds extend to other wavelengths. All three clouds are strong sources of IR radiation, indicating the presence of hot dust. They also have strong molecular counterparts and are massive, likely with masses greater than or equal to 5 × 10⁶ M_⊙. Finally, energetic phenomena occur in their vicinity, as evidenced by the strong radio continuum.

5.3.2. Shocked molecular clouds

Our new results indicate that clouds B, C, and D are giant molecular clouds shocked by the AGN jet, as evidenced by the high-energy forbidden lines observed in our NIR spectra and the nearby strong radio sources. The presence of hot dust and molecular content indicates that these clouds have not been completely destroyed by the jet yet, and they could be transient structures.

If the IR peak indeed traces the position of the nucleus (Pfuhl 2020; Vermot et al. 2021; Gámez Rosas et al. 2022), then cloud B is in the line of sight of the nucleus. If it is interacting with the jet, as evidenced by this work, it means that it has reached the actual position of the nucleus. As such, cloud B could simultaneously be the mass reservoir for the current accretion episode and the dense obscurer giving the nucleus its type 2 properties. Both the hot phase of the molecular content observed with SINFONI through the H₂ rovibrational line (Müller Sánchez et al. 2009) and the colder phase traced by various ALMA molecular lines (García-Burillo et al. 2019; Impellizzeri et al. 2019; Imanishi et al. 2020) point toward complex kinematics that are incompatible with a simple rotation profile and probably result from the superposition of rotation, counter rotation, and outflow kinematics. Thus, if cloud B is indeed currently fueling and obscuring the nucleus, it has not reached a steady orbital state around the central mass.

Cloud D is part of the CND, and clouds E and F (not discussed in this paper) are also part of the CND, providing evidence that the CND shocked region is extended. The presence of strong IR emission and molecular tracers at the position of the clouds reveal that the interaction with the jet is not sufficient to destroy the densest molecular clouds of the CND.

Cloud C is located in the empty region between the CND and the nucleus, and it appears to be an extended structure as shown in the high angular resolution [O III] HST/FOC image (Macchetto et al. 1994) and the H₂ rovibrational line image (Müller Sánchez et al. 2009). A kinematics analysis in he latter reference indicates that the structure is in a highly elliptical orbit around the nucleus and is streaming toward it. The associated radio component C is located on its eastern edge, which is the leading side of its orbit, as per Müller Sánchez et al. (2009). Our results indicating that the emission lines come from a shocked region are consistent with the scenario of cloud C being a tidally disrupted molecular cloud that is shocked and ionized as its orbit brings it into contact with the outflow.

In this scenario, clouds D, C, and B would chronologically trace the evolutionary stages of the gas reservoir for the AGN activity. The CND (cloud D) would constitute the bulk of the molecular reservoir. Its interaction with the outflow would result in a loss of angular momentum, causing some clouds to detach from it and fall into the gravitational potential (cloud C). Ultimately, these clouds would reach the nucleus (cloud B) and trigger the accretion onto the supermassive black hole (SMBH), potentially obscuring it in the process.

5.3.3. Cloud C: A potential secondary AGN

Cloud B is located at the position of the IR peak, which is reliably associated with S1 (Pfuhl 2020; Gámez Rosas et al. 2022), the AGN nucleus (Gallimore et al. 2004). Given the similarities observed between clouds B and C in their NIR emission lines, radio flux and spectrum, and IR flux, it is reasonable to consider that cloud C could also host an active SMBH.

In the IR images presented in Gratadour et al. (2006), the authors discuss the challenge of explaining the high temperatures observed in IR-1B with radiation from the central source and suggest that the interaction of the jet with its environment or the presence of very small grains could account for it. However, the presence of an active SMBH in cloud C would provide a simpler explanation for this temperature excess.

While component S1 is identified as the nucleus in Gallimore et al. (2004), due to an inverted spectrum below 5 Ghz and a flat spectrum above as well as the presence of bright H₂O masers in a rotating pattern, component C is not identified as a nucleus despite having an almost identical spectrum and flux as S1 and the presence of H₂O masers (without a rotating pattern). Component C is instead associated with the interaction between the jet and a molecular cloud, due to the non-exact alignment with optical emission lines and the apparent deviation of the jet at this position. However, recent observations by Morishima et al. (2023) have revealed a ring-like structure of H₂O masers around component C, which could be indicative of rotation. While this does not rule out the possibility that component C results from the interaction of the jet with a molecular cloud, the authors suggest that it could also indicate the presence of a rotating disk around a ∼10⁶ M_⊙ SMBH. This mass estimate is consistent with the velocity gradient measured on the [Si VI] emission line in our SINFONI data (see Fig. 4).

Fig. 4. Cloud C: velocity map based on the measured Doppler shift of the [Si VI] emission line. The contours represent flux levels.

At radio wavelengths, the properties of components S1 and C in NGC 1068 are comparable to those of components N1 and S in NGC 6240 in terms of flux, geometry, and (to a lesser extent) spectral slope (Gallimore & Beswick 2004). In the case of NGC 6240, these two radio sources are identified as a double AGN. The SINFONI observation of NGC 6240 (Ilha 2016) also shows the peak of the coronal emission lines at the position of the two AGN, closely but not perfectly lined up, similar to NGC 1068. The two nuclei of NGC 6240 are also strong sources of NIR to mid-IR radiation (Max et al. 2005; Mori et al. 2014).

The decisive argument to classify NGC 6240 as a double AGN was the association of the two strong IR and radio sources to two hard, luminous X-ray nuclei discovered with Chandra (Komossa et al. 2003). However, in the case of NGC 1068, clouds B and C are too close together for X-ray telescopes to resolve them. Nevertheless, the Chandra image presented in Young et al. (2001) reveals an extension of the nuclear component in the direction of cloud C, and the spectral modeling done in the same paper highlights the difficulty in fitting the X-ray spectrum with a single hot plasma model. These elements are consistent with the presence of a secondary source.

In essence, we contend that even without S1, the central region of NGC 1068 would still be classified as an AGN based on the properties of cloud C, such as its strong inverted radio continuum, its ring-like distribution of masers, and its strong IR emission. Our results corroborate this interpretation, as they reveal strong similarities between the emission line spectra of clouds B and C, with the former being identified as an AGN. However, this argument is not conclusive since the shocked atomic content detected in our study could also arise from a shock with the jet. Given that cloud C is in the outflow region, this hypothesis is favored until further evidence for AGN activity in cloud C (such as resolved X-ray emission) is found. We plan to investigate this question further in the near future.

6. Conclusions

The combination of SPHERE and SINFONI data enabled us to produce a comprehensive YJHK spectrum for three emission line regions named B, C, and D in the central arcsecond of NGC 1068. Our analysis reveals that these regions exhibit similar emission line spectra, suggesting that they belong to the same class of objects. The various modeling techniques utilized in this study produced consistent results for all three objects.

Our investigation involving the use of CLOUDY to produce synthetic emission line spectra for several models and the comparison of these spectra with the full observations revealed the presence of multiple phases in the clouds. The temperatures of these phases range from 10^4.8 K to 10⁶ K. Our findings indicate that a central AGN photoionization model is not able to accurately reproduce the variety of emission lines observed in the spectra. However, the MAPPINGS III fast radiative shock model is able to reproduce the observed emission line spectra well, with shock velocities of v_s ∼ 1000 km s⁻¹ and near equipartition magnetic fields.

Our study suggests that the three objects are likely to be giant molecular clouds that are either orbiting or streaming toward the nucleus and are being shocked by the same radio jet as they enter the outflowing bicone. It is possible that most of the properties associated with the nucleus and its torus, such as mid-IR and radio continuum emissions as well as atomic and molecular lines, originate from cloud B, which is coincident with the nucleus and may be currently feeding the central SMBH. However, given the similarities between clouds B and C, together with their shared AGN properties, such as high IR flux, radio continuum, H₂O masers, and M_dyn ≫ 10⁶ M_⊙, it is also possible that both objects host an AGN, with cloud C being a secondary AGN orbiting cloud B.

Acknowledgments

We thank Éric Lagadec and the SPHERE consortium for carrying out this observation as part of the GTO Other Science program. This work was made possible by the support of the international collaboration in astronomy (ASU mobility) with the number CZ.02.2.69/0.0/0.0/18_053/0016972 and the institutional project RVO:67985815. ASU mobility is co-financed by the European Union.

References

Appendix A: Comparison between models and observations

The following appendices contain additional information supporting the conclusions drawn in the main body of the text.

For each ionized region, Appendix A presents a comparison between the observed spectrum and the various models discussed in the analysis: PC1, PC2, PC3, AGN, and SHK. Figure A.1 presents the results for cloud B, while Figs. A.2 and A.3 present the results for clouds C and D, respectively. Appendix B presents the reduced chi-square maps for models PC1 and PC2. For model PC1, Fig. B.1 shows the 2D maps as a function of density and temperature for clouds B, C, and D, while Figs. B.2 and B.3 show the 1D plots for the same parameters. For model PC2, Fig. B.4 shows the 2D maps of the reduced chi-square as a function of T₁ and T₂ for the three regions. Appendix C presents the reduced chi-square for the central photoionization model AGN. Figure C.1 shows the 2D maps for the parameters of the central source (α_OX and T) for the three regions, while Fig. C.2 shows the chi-square as a function of the density of the cloud. Appendix D presents the reduced chi-square for the fast radiative shock model SHK. Figures D.1, D.2, and D.3 show the 1D chi-square as a function of shock velocity, preshock density, and transverse magnetic field, respectively, for all three clouds. In Appendix E, we compare the measured flux of each line in clouds B, C, and D with the prediction of the SHK model. Appendix F presents the emitting conditions (temperature, density) for a few selected emission lines of interest: [Fe II], [Si VI], [Si X], and [C I]. Finally, Appendix G presents the [Si VI] Doppler shift maps for clouds C and D in Figs. G.1 and G.2, respectively.

Fig. A.1. Cloud B: Comparison between the observed spectrum and the various models. The top panel shows the observed spectrum in black with the estimated uncertainty overlaid in transparent red. The five following plots show the prediction from models PC1, PC2, PC3, AGN, and SHK, in that order, with the observation overlaid in transparent blue.

Fig. A.2. Cloud C: Comparison between the observed spectrum and the various models. The top panel shows the observed spectrum in black with the estimated uncertainty overlaid in transparent red. The five following plots show the prediction from models PC1, PC2, PC3, AGN, and SHK, in that order, with the observation overlaid in transparent blue.

Fig. A.3. Cloud D: Comparison between the observed spectrum and the various models. The top panel shows the observed spectrum in black with the estimated uncertainty overlaid in transparent red. The five following panels show the prediction from models PC1, PC2, PC3, AGN, and SHK, in that order, with the observation overlaid in transparent blue.

Appendix B: Physical conditions

B.1. Single temperature and density

Fig. B.1. Model PC1: Reduced chi-square between the model and the observation spectra as a function of temperature and density for clouds B (top), C (middle), and D (bottom).

Fig. B.2. Model PC1: Reduced chi-square between the model and the observation spectra as a function of temperature for clouds B (top), C (middle), and D (bottom).

Fig. B.3. Model PC1: Reduced chi-square between the model and the observation spectra as a function of density for clouds B (top), C (middle), and D (bottom).

B.2. Double temperature, n = 2 cm⁻³

Fig. B.4. Model PC2: Reduced chi-square between the model and the observation spectra as a function of T1 and T2 for clouds B (top), C (middle), and D (bottom).

Appendix C: Central photoionization

Fig. C.1. Model AGN: Reduced chi-square between the model and the observation spectra as a function of temperature and αOX for clouds B (top), C (middle), and D (bottom).

Fig. C.2. Model AGN: Reduced chi-square between the model and the observation spectra as a function of density for clouds B (top), C (middle), and D (bottom).

Appendix D: Fast radiative shock

Fig. D.1. Model SHK: Reduced chi-square between the model and the observation spectra as a function of the shock velocity for clouds B (top), C (middle), and D (bottom).

Fig. D.2. Model SHK: Reduced chi-square between the model and the observation spectra as a function of the preshock density for clouds B (top), C (middle), and D (bottom).

Fig. D.3. Model SHK: Reduced chi-square between the model and the observation spectra as a function of the preshock transverse magnetic field for clouds B (top), C (middle), and D (bottom).

Appendix E: Line fluxes from model SHK

Appendix F: Physical conditions for selected emission lines

Fig. F.1. Emitting conditions for [Fe II] at 1.644 μm: $\frac{[FeII]}{H\beta }$ normalized to one for its maximum value.

Fig. F.2. Emitting conditions for [Si VI] at 1.965 μm: $\frac{[SiVI]}{H\beta }$ normalized to one for its maximum value.

Fig. F.3. Emitting conditions for [Si X] at 1.430 μm: $\frac{[SiX]}{H\beta }$ normalized to one for its maximum value.

Fig. F.4. Emitting conditions for [C I] at 0.985 μm: $\frac{[CI]}{H\beta }$ normalized to one for its maximum value.

Appendix G: [Si VI] Doppler shift

Fig. G.1. Cloud C: Velocity map based on the measured Doppler shift of the [Si VI] emission line.

Fig. G.2. Cloud D: Velocity map based on the measured Doppler shift of the [Si VI] emission line.

All Tables

All Figures

Fig. 1. Top: schematic representation of the main components in the central region of NGC 1068: the CND, the outer edges of the outflow bicone, and the three regions (B, C, and D) analyzed in this work. The field of view of the SINFONI observation is indicated with a pink square. Bottom: [Si VI] emission line map in the central region of NGC 1068 obtained with the SINFONI observation (pixel scale 50 mas) and position of the SPHERE slit (red lines, roughly aligned with the jet orientation). The coordinates are centered on the position of the K band continuum maximum. Following the naming convention from Evans et al. (1991), the three regions studied in this paper are cloud B (in the center; at the same position as the radio component S1 and the compact torus), cloud C (0.25″ north), and cloud D (0.7″ north).

In the text

	Fig. 2. Combined YJHK spectra of the clouds with increasing distance from the nucleus. From top to bottom, the panels show a distance of 0.0″ (B/S1/nucleus), 0.25″ (C), and 0.7″ (D). The SPHERE data are displayed in green, the SINFONI data in blue, and the averaged overlap in red. The estimated uncertainty is indicated as a gray area.
In the text

	Fig. 3. Cloud D: comparison between the observed spectrum (top, black) and the best model (bottom, red). The numerical values of the fluxes from the model’s main emission lines are compared with those from the observation in Table E.3 from the appendix.
In the text

	Fig. 4. Cloud C: velocity map based on the measured Doppler shift of the [Si VI] emission line. The contours represent flux levels.
In the text

	Fig. A.1. Cloud B: Comparison between the observed spectrum and the various models. The top panel shows the observed spectrum in black with the estimated uncertainty overlaid in transparent red. The five following plots show the prediction from models PC1, PC2, PC3, AGN, and SHK, in that order, with the observation overlaid in transparent blue.
In the text

	Fig. A.2. Cloud C: Comparison between the observed spectrum and the various models. The top panel shows the observed spectrum in black with the estimated uncertainty overlaid in transparent red. The five following plots show the prediction from models PC1, PC2, PC3, AGN, and SHK, in that order, with the observation overlaid in transparent blue.
In the text

	Fig. A.3. Cloud D: Comparison between the observed spectrum and the various models. The top panel shows the observed spectrum in black with the estimated uncertainty overlaid in transparent red. The five following panels show the prediction from models PC1, PC2, PC3, AGN, and SHK, in that order, with the observation overlaid in transparent blue.
In the text

	Fig. B.1. Model PC1: Reduced chi-square between the model and the observation spectra as a function of temperature and density for clouds B (top), C (middle), and D (bottom).
In the text

	Fig. B.2. Model PC1: Reduced chi-square between the model and the observation spectra as a function of temperature for clouds B (top), C (middle), and D (bottom).
In the text

	Fig. B.3. Model PC1: Reduced chi-square between the model and the observation spectra as a function of density for clouds B (top), C (middle), and D (bottom).
In the text

	Fig. B.4. Model PC2: Reduced chi-square between the model and the observation spectra as a function of T1 and T2 for clouds B (top), C (middle), and D (bottom).
In the text

	Fig. C.1. Model AGN: Reduced chi-square between the model and the observation spectra as a function of temperature and αOX for clouds B (top), C (middle), and D (bottom).
In the text

	Fig. C.2. Model AGN: Reduced chi-square between the model and the observation spectra as a function of density for clouds B (top), C (middle), and D (bottom).
In the text

	Fig. D.1. Model SHK: Reduced chi-square between the model and the observation spectra as a function of the shock velocity for clouds B (top), C (middle), and D (bottom).
In the text

	Fig. D.2. Model SHK: Reduced chi-square between the model and the observation spectra as a function of the preshock density for clouds B (top), C (middle), and D (bottom).
In the text

	Fig. D.3. Model SHK: Reduced chi-square between the model and the observation spectra as a function of the preshock transverse magnetic field for clouds B (top), C (middle), and D (bottom).
In the text

	Fig. F.1. Emitting conditions for [Fe II] at 1.644 μm: $\frac{[FeII]}{H\beta }$ normalized to one for its maximum value.
In the text

	Fig. F.2. Emitting conditions for [Si VI] at 1.965 μm: $\frac{[SiVI]}{H\beta }$ normalized to one for its maximum value.
In the text

	Fig. F.3. Emitting conditions for [Si X] at 1.430 μm: $\frac{[SiX]}{H\beta }$ normalized to one for its maximum value.
In the text

	Fig. F.4. Emitting conditions for [C I] at 0.985 μm: $\frac{[CI]}{H\beta }$ normalized to one for its maximum value.
In the text

	Fig. G.1. Cloud C: Velocity map based on the measured Doppler shift of the [Si VI] emission line.
In the text

	Fig. G.2. Cloud D: Velocity map based on the measured Doppler shift of the [Si VI] emission line.
In the text