© D. Rouan et al. 2019

1. Introduction

In an active galactic nucleus (AGN), considered here as the very central region encompassing the stellar cluster and a dense molecular core, the feeding of the central engine (hereafter CE), i.e. the accretion disc, requires a matter accretion rate that ranges from 0.3 in Seyfert galaxies (Bian & Zhao 2003) to 10 M_⊙ yr⁻¹ (Hopkins et al. 2012) for the most luminous quasars. How the very final transfer of matter is achieved from kpc scale to the sub-parsec size of the accretion disc is still a matter of debate since the mechanism of loss of angular momentum at scales below 100 pc has not been clearly established (see Goodman 2003 for a review), nor is there any clear relationship between this feeding mechanism and the likely existence of an obscuring torus, the paradigm proposed by Antonucci & Miller (1985) years ago to interpret the Seyfert 1/2 duality, and that still remains valid in the light of recent ALMA or SPHERE polaro-imaging observations (García-Burillo et al. 2016; Gratadour et al. 2015). One of the most difficult problems for hydrodynamical models is to take into account a huge range of scales, from 10 kpc in spiral arms down to 0.1 pc in the accretion disc while integrating all the possible processes such as cooling, star formation, and self-gravity of the gas, which are all highly non-linear.

Hopkins et al. (2012), using results from their multi-scale hydrodynamical model (Hopkins & Quataert 2010, 2011), recently proposed that the formation, through instabilities, of lopsided eccentric gas plus stellar discs about 10 pc in size would be the dominant cause for angular momentum transport below this scale. The asymmetry of the stellar disc would be the main cause of the strong torque exerted on the gas. The model also predicts that surface densities of gas and stars both follow cuspy distributions.

Performing observations in the visible to IR that could test the validity of this proposal, especially evidencing the stellar disc, is not obvious because there are several limits: i) we are dealing with scales that are largely subarcsec (1 pc = 20 mas at 10 Mpc); ii) the central region is very obscured by the huge concentration of dust; and iii) the bright dazzling AGN totally dominates the flux, making it difficult to detect fainter structures associated with warm dust emission or stellar population in the immediate neighbourhood. Of course millimetric interferometry, especially using ALMA, could solve the two first issues, but it is not suited to probe stellar emission. On the other hand, the combination of adaptive optics (AO) on large telescopes and imaging in the near-IR (NIR) gives us the opportunity to tackle the problem, at least on nearby AGNs. This is the goal of this paper where we present an analysis of adaptive optics images of NGC 1068 in K_s and H bands, and in several narrow-band filters, obtained with SPHERE on the VLT.

NGC 1068 is a Seyfert 2 galaxy that has been the target for many studies about AGN since it is one of the closest active nuclei (14.4 Mpc, following Bland-Hawthorn et al. 1997) and therefore one of the brightest. One consequence is that the nucleus is bright enough to be used as the guide source for the wavefront sensor of the AO system, allowing us to obtain high angular resolution images in the NIR, for example with PUEO-CFHT (Gratadour et al. 2003), NAOS-VLT (Rouan et al. 2004; Gratadour et al. 2006; Prieto et al. 2005), SINFONI-VLT (Davies et al. 2007; Müller Sánchez et al. 2009), or more recently using SPHERE-VLT (Gratadour et al. 2015).

The later observations were conducted during the SPHERE Science Verification (SV) programme in December 2014. The H- and K_s-band broad-band images, showed a polarisation angle map with a clear centro-symmetric pattern, tracing both parts of the ionisation bicone. More importantly, they featured a central non-centro-symmetric pattern approximately 60 pc × 20 pc wide, with aligned polarisation, that was interpreted as the trace of the outer envelope of the torus, revealed through a double scattering process (Gratadour et al. 2015; Grosset et al. 2018). It is this set of images, complemented by images obtained in narrow-band filters in September 2016, again with SPHERE, that we re-processed, but without using the polarisation information. In this paper we concentrate on the radial distribution of brightness.

The paper is organised as follows. The data processing and radial profile fitting is presented in Sect. 2. The discussion on the power-law radial profiles at H and K_s band and their interpretation in terms of a reddened cuspy stellar distributions through a simple transfer model is presented in Sect. 3 after some discussion about various classes of cusps observed in the centre of galaxies. In Sect. 4 we discuss an alternative and more acceptable model where warm/hot dust heated by the central engine is mixed to stellar emission. Finally, we examine the central point-like source in Sect. 5 as well as the constraints brought on the dust opacity by its luminosity, before concluding in Sect. 6.

2. Data set and data processing

The K_s- and H-band images obtained with SPHERE in December 2014 are shown in Fig. 1, using a Asinh scaling law for the intensity that enhances the low flux levels. Details about observations and data reduction have been published in Gratadour et al. (2015). The pixel size is 12.25 mas, corresponding to 0.855 pc assuming a distance of 14.4 Mpc. The field of view is 6.27 × 6.27 arcsec², or 438 × 438 pc², and the angular resolution, as measured on the point spread function (PSF), is 58 mas and 53 mas in K_s and H, respectively.

Fig. 1. Ks-band (left) and H-band (right) images of the central 438 × 438 pc of NGC 1068, obtained with SPHERE, where 1 arcsec is 70 pc. The intensity scaling corresponds to an Asinh law. For the Ks-band image the flux at peak (white area) is 2.3 × 10−3 Jy pixel−1, and in the light blue area it is ≈6.0 × 10−6 Jy pixel−1. For the H-band image the flux at peak (white) is 2.0 × 10−4 Jy pixel−1, and in the light blue area it is ≈7.5 × 10−6 Jy pixel−1.

In the H-band image can be seen the elongated disc, slightly asymmetric with respect to the CE, well defined by a strong peak. The disc size is about 135 pc × 180 pc, choosing as the limit the distance where the slope of the radial distribution flattens. Assuming that the corresponding iso-contour delineates the disc, the barycentre of this contour is offset by 6 pixels (5.1 pc) to the west and 1 pixel (0.9 pc) to the north with respect to the central peak. In the K_s-band image, the asymmetry is less pronounced and, for instance, the centre of gravity of the contour delineating about the same area is only shifted by 2.5 pixels (2.1 pc) and 1.5 pixels (1.3 pc) respectively to the west and north.

We should now ask an important question: what is the source of emission at H band and K_s band? The answer is not trivial, especially concerning the K_s band. We expect an important contribution from stellar emission, but we also expect a significant contribution from very hot dust, close to sublimation temperature, arising from the internal radius of the likely torus. For instance, Gratadour et al. (2003) showed clearly that the steeper slope of the NIR spectrum when approaching the central peak is nicely explained by this very hot dust. In the following, our main concern is to discriminate between the respective contributions of the two sources of radiation.

Let us first examine the medium scale. We present f in Fig. 2 the median brightness profile in H band extending up to 430 pc, i.e. well beyond the molecular ring of r ≈ 200 pc detected by Müller Sánchez et al. (2009) and confirmed by ALMA (García-Burillo et al. 2014). In the same figure, two fits are represented: a power-law fit for the most internal part (r ≤ 110 pc) and an exponential fit for the external part (up to r = 430 pc). Both fits are fairly good and, assuming that it is essentially the stellar component seen at H, we conclude that when going inward, a disc distribution (here of scale length R_c = 410 pc) is superseded by a cuspy (i.e. power-law) distribution of exponent 1.2 for r ≤ 110 pc.

Fig. 2. Median radial brightness profile around the central peak, measured in the H-band. Blue solid line: Brightness profile in H band up to 430 pc, estimated by performing a median of the pixels on successive circles centred on the central peak. Red dashed line: Power-law profile r−1.2 fitted on the central distribution. Green dashed line: exponential profile fitted on the external observed distribution.

In Fig. 3 we plot the median radial brightness profile around the central peak, measured in the K_s-band and H-band images, together with the PSF profile at K_s band measured on the nearby reference star BD-00413 whose visible magnitude is close to that of the nucleus of NGC 1068, so that the AO performances can be considered very similar. Two features should be noted. First, the H-band profile is almost perfectly fitted by a power law of exponent −1.2; we come back to the meaning of this behaviour in the next section. Second, the K_s-band profile wings, beyond a radius of about 17 pc, are also well fitted by a power law, but of exponent −2.0. Interestingly, a fairly good fit to the whole data set is provided by a combination of this power-law profile truncated at r <  4 pc plus the PSF profile properly weighted. This appears as the blue curve in Fig. 3 superimposed on the K_s-band actual profile. In particular, we note that the small bumps on the K_s-band profile, close to the peak, are well fitted by the first and even the second Airy rings that clearly appear on the PSF profile. The straightforward interpretation of this fit is that the K_s-band profile results from the contribution of a central point source and of a smoother distribution that is directly related to the H-band emission. We tested other combinations by convolving the PSF with a Gaussian of various widths, but beyond a FWHM of 2 pixels (1.7 pc) the fit degrades. The central source is thus practically point-like, and in any case its radius is not significantly larger than 0.8–1 pc.

Fig. 3. Median radial brightness profile around the central peak. Solid red: measured at Ks; purple: PSF at Ks; red dash: difference between Ks band and weighted PSF; green: truncated power-law profile in r−2; solid blue: H-band profile; blue dash: power-law profile in r−1.2.

In Sect. 5 we come back to the interpretation of this important result, namely that there is a quasi-unresolved central point source in K_s band with no obvious counterpart in H. We concentrate in the next two sections on the power-law, or cuspy, component.

In Sept. 2016 we also obtained images with SPHERE in three narrow-band filters denoted cont–H (1.558 μm, Δλ = 0.0245 μm), cont–K1 (2.103 μm, Δλ = 0.033 μm), and cont–K2 (2.267 μm, Δλ = 0.034 μm). The details of the data reduction will be described in Grosset et al. (in prep.), a paper dedicated to the analysis of polarisation versus wavelength, but can also be found in Grosset (2017). The data processing is essentially the same as that used for the H- and K_s-band observations of December 2014.

We plot in Fig. 4 the median averaged profiles at the four wavelengths on a log-log scale. The following can be seen: i) the fairly linear behaviour of all profiles from 4 to 100 pixels (3.5 to 85 pc), as illustrated by the superimposed dash lines, and ii) the absolute value of the slope increasing with wavelength from H band to K_s band. The best fit power laws for the range 5–60 pc have exponents −1.20, −1.36, −1.95, −2.11, and −2.03 for H-band, cont–H, cont–K1, cont–K2, and K_s-band, respectively.

Fig. 4. Median radial brightness profile around the central peak for the cont–H (blue), H-band (green), cont–K1 (cyan), cont–K2 (magenta), and Ks-band (red) images. The log-log scale shows a power-law behaviour for all four bands. The cont–K1, cont–K2, and Ks-band fluxes have been scaled on the H-band flux at the centre. The dashed line on each curve is the best linear fit for the range 5−60 pc.

There are several possible origins for the power-law emission, and we examine the different cases in the following sections; however, the monotonic variation of the slope with respect to wavelength on a sample of five different filters is a strong indication that the cause of this variation is a continuously variable phenomenon rather than some spectral features such as emission lines. With this hypothesis in mind and taking into account that the signal-to-noise ratio is significantly lower in narrow-band images, we have chosen in the following to consider only the H-band and the K_s-band in our analysis; thanks to the broad range in wavelength coverage, any effect depending on wavelength will thus be clearly enhanced.

The first hypothesis, which could be well in line with the quite similar behaviour of the H-band and K_s-band radial brightness profiles and which corresponds to a situation often observed in the central region of galaxies (Böker 2008), is that in both bands we sense the emission from a stellar cusp, also called a nuclear stellar cluster. The second case, which is generally invoked in the case of NGC 1068 (Thatte et al. 1997; Gratadour et al. 2003) and other AGNs (Leja et al. 2018), is that there is warm to hot dust, heated by the CE and radiating in the infrared, especially at K_s. Of course, some mixing of this dust thermal radiation with the photospheric emission from a stellar cluster is likely, especially at H.

3. A reddened stellar cusp?

To start, we can assume that the H-band power-law profile is largely due to the stellar emission. First, because of this quasi-perfect power-law behaviour, which is very similar to what is observed in many nuclei of galaxies (see next section), while it could hardly be explained by dust emission extending so far from the CE, and second, because previous spectroscopic observations by Thatte et al. (1997) clearly identified the 1.6 μm 3→6 and the 2.29 μm 0→2 band-head CO feature in a central core of FWHM of size around 45 pc (0$\stackrel{″}{.}$66). However, Thatte et al. (1997) estimate the contribution from stellar emission to be only 1/3 of the total within an aperture of 1″ radius. We note that the angular resolution they reached is only 0$\stackrel{″}{.}$9, thus far below the SPHERE performance of 0$\stackrel{″}{.}$06 in the K_s band (see Sect. 2).

For the K_s-band profile, the contribution of hot dust from a compact central source was clearly identified as early as 1987 (Chelli et al. 1987). Thatte et al. (1997), based on their 0$\stackrel{″}{.}$4 resolution speckle observations estimated that essentially all the flux at K_s-band could be attributed to the point-like source. However, this is not what we observe on the K_s-band profile in Fig. 1, where the aisles extend far beyond those of the PSF. On the other hand, the wings beyond 17 pc (0$\stackrel{″}{.}$24) are also very well fitted by a power law, which leads to the same reasoning as for the H band, i.e. that it could correspond to stellar emission. If this hypothesis is correct, then it raises the question of why the exponents are so different between K_s band and H band.

Another question raised by this profile behaviour is how the cusp characteristics compare with those of cusps actually found in early-type galaxies or in LIRGs/ULIRGs. As regards the size of the cusp, by integrating the flux in discs of increasing diameter, we have determined a characteristic half-light radius of 99 pc (1.420″) and 17 pc (0.245″) for the H-band and K_s-band cusps, respectively, thus two very different scales that need to be explained.

Compared to the structures observed at millimetre wavelengths, either in molecular lines (CO, HCN) or in the continuum (García-Burillo et al. 2014, 2016; Imanishi et al. 2018), we note that the extent of the cusp at H band reaches the inner radius of the off-centred circumnuclear molecular disc. This means that we do not expect significant extra extinction from the disc itself since the stellar cusp is almost entirely within its inner contour.

3.1. Emission lines from the ISM

We consider the possible role of differential dust extinction to explain the log-log H-band versus K_s-band slope difference in Sect. 3.3, but first let us examine another potential explanation: an extra contribution to the K_s-band flux by lines emitted by the ISM that would be significant and would feature a gradient. We can think of two kinds of line emission, either the emission of the (1-0) S(1) rotational-vibrational line (2.1218 μm) of some excited H2, or the Br_γ line (2.1655 μm) associated with the narrow- or even the broad-line region. The first case must be excluded since the high angular resolution map of the line, obtained by Müller Sánchez et al. (2009) with the spectro-imager SINFONI-VLT, clearly shows that the line emission is confined to a compact region (about 0.15″ diameter) at 0.15″ to the south-east of the central engine. Moreover, the same authors provided a spectrum of the region 0.3″ × 0.1″ centred on the nucleus that clearly shows that the contribution of the (1-0) S(1) line to the flux in the whole K_s band is only 5 × 10⁻⁴. The same spectrum indicates that the Br_γ line contribution to the same band is approximately 3 × 10⁻³ and thus should also be excluded as an explanation for the differences in profile slopes between H band and K_s band. The same conclusion can be reached considering the long slit spectrum (0.3″ slit width) in band K_s band obtained by Gratadour et al. (2003).

3.2. Variety of stellar cusps

In the literature, the notion of central cusp seems to refer to a variety of situations or contexts and we first try to make this variety more explicit.

We can identify a first category where the cusp is associated with some compact central stellar cluster. For instance, it corresponds to an extra light distribution compared to a Sérsic distribution extrapolated to the very centre in early-type galaxies (Kormendy et al. 2009; Lauer et al. 2007). The expression nuclear cluster (NC) may also be used, either for early- or late-type galaxies (Böker et al. 2004). Lauer et al. (2007) use the notion of power-law galaxies (with steep brightness increase) as opposed to core galaxies where the central brightness is nearly constant. Côté et al. (2006) use the denomination compact central nucleus. We consider in the following that all these denominations represent the same kind of object, namely a central compact stellar cluster analogous to a globular cluster, even if it differs from this class by a higher luminosity (see e.g. Böker et al. 2004), typically 3.5 mag brighter. As stated by Böker (2008), NCs are common: the fraction of galaxies with an unambiguous NC detection is 75% in late-type (Scd-Sm) spirals (Böker et al. 2002), 58% in earlier type (S0-Sbc) spirals (Balcells et al. 2007), and 70% in spheroidal (E and S0) galaxies (Graham & Guzmán 2003; Côté et al. 2006).The median characteristics that can be for instance derived from Table 2 of Lauer et al. (2007) is a half-light radius of 7 pc, within a range of 2−44 pc, and an exponent of 0.62, within a range between 0 and 1.1. Côté et al. (2006) find for Virgo early-type galaxies a similar and consistent radius of 4.2 pc for a range extending to 62 pc. These typical characteristics do not depend strongly on the type of galaxy; for instance, Böker et al. (2004) find, for a sample of late-type galaxies, a histogram of radii with a first peak at 3 pc and another one at 8 pc, the range extending to 30 pc. However, it has been clearly established that faint ellipticals (M_V ≤ 20) exhibit a more cuspy profile than massive ones (Kormendy et al. 2009; Côté et al. 2006), a feature that likely translates into a different history of merging.

The second and distinct category of stellar cusps we have identified in the literature is more clearly related to recent or ongoing starbursts. We refer in particular to the study by Haan et al. (2013) of a large sample of LIRGs and ULIRGs where HST images in the NIR were used to derive the characteristics of the nuclear stellar cusps they found in nearly all objects. We analysed the data compiled by those authors concerning ≈80 LIRG and ULIRG objects, a large fraction of which are mergers at various phases of the process. In most of the objects the authors found a cuspy profile in the nuclear region. They provide a table of the various properties of the selected objects, including radius and NIR luminosity of the cusp, as well as the exponent of the cuspy brightness distribution. Clearly the average characteristics are distinct from those of the class we discussed above: we find a median half-light radius of 475 pc within a range between 46 pc and 2.4 kpc, while the exponent of the cusp has a median value of 0.80 within a range between 0.1 and 1.7. The cusp radius is remarkably much larger than for the nuclear cluster of the class identified above, and the exponent is significantly larger and can reach high values.

Finally, the third category of cusps that is described in the literature corresponds to the compact stellar population that is distributed around a massive black hole in the centre of a galaxy. Some theoretical works have explored the dynamics of this peculiar system (see e.g. Preto 2010 for a review), but there are no well-established observational facts, with the notable exception of the centre of the Milky Way (Genzel et al. 2003; Schodel et al. 2017). These nuclear stellar cluster (NSC) cusps are characterised by half-mass or half-light radii of 1−5 pc and are in principle the result of dynamical effect, namely two-body relaxation through frequent encounters of stars or stellar BH in the strong gravitational field of a massive BH. In particular, Bahcall & Wolf (1976) established the first solid analytical expression of the stellar distribution around a massive BH and showed that it must be cuspy with an exponent 1.75 (7/4). Refinements to the theory were brought later by numerous teams (for a review, see e.g. Preto 2010), especially the fact that there must be a mass segregation with different exponents according to the mass of the stellar objects, the most massive stellar BH having the largest exponent (1.8) and main sequence stars a lower value (typically 1.3) (Freitag et al. 2006). These exponents are thus clearly larger than that found in the two previous classes we examined, and closer to the ones observed in NGC 1068, but if we consider both the predicted small half-mass radius and the too large relaxation time T_r required to reach a cuspy distribution – due to the large mass of the BH – we can no longer consider this case for NGC 1068.

To make the distinction between the two classes more readable in the following, we will call the first the nuclear cluster cusp (NCC) class, and the second the central starburst cusp (CSC) class. Our goal is to try to determine if there are clues that allow the identification of the cusp we observe with one of these classes, or a combination of the two.

The rather large exponent we found for the cusp in NGC 1068 and the rather wide half-light radius clearly rule out the case where we would observe a mere NCC; such large exponents are never observed for this class.

For the CSC class, NGC 1068 shares, at least for the H-band cusp, some characteristics, with an exponent and a radius that are indeed significantly different from the median values, but within the observed ranges for this class. Nevertheless, the K_s-band profile with its small half-light radius and large exponent could weaken this statement, but we show in Sect. 4 that the dust emission cannot be neglected. However, there is such a large dispersion of characteristics in the sample of Haan et al. (2013) that the question deserves closer examination.

To go one step further, starting from the data set compiled by Haan et al. (2013), we tried to determine whether some general law could relate the parameters of the cusp in CSC. The purpose was to see if NGC 1068 could obey such a general law even if it is not considered a LIRG. If the answer were yes, it would be a strong indication that the cusp is the result of a starburst evolution rather than the result of a dynamical effect around a massive BH. We first selected the 62 objects where the brightness distribution is indeed cuspy (γ ≠ 0) and where the information on the three parameters: radius, NIR luminosity, and exponent of the cusp, was complete (e.g. not given as an upper or lower limit). We looked for a fit of γ as a linear combination of the log of the NIR cusp luminosity L_cusp and the log of the cusp radius R_cusp. We found that the law

$$\begin{array}{c}\hfill \gamma =1.133\times [log10(L_{\mathrm{cusp}}/R_{\mathrm{cusp}})-6.63]\end{array}$$(1)

relating the three parameters was rather robust, with a covariance between the actual value and the predicted value of 0.85. The robustness of the fit can be seen in Fig. 5, where the luminosity of each selected object is plotted against the value predicted by the law using the two other parameters.

Fig. 5. Estimated cusp luminosity of a selection of the objects studied by Haan et al. (2013), as predicted by the phenomenological law we found (see text) vs the actual cusp luminosity. The blue line is the linear regression. The cross indicates the objects that are not mergers or are in a wide system.

The straightforward interpretation of the law is that at a given luminosity the smaller the radius of the cusp, the steeper the distribution and, conversely, a lower luminosity corresponds to a smoother cusp. We note that in Fig. 5, we have flagged with an additional cross the objects that do not belong to a merger or that are members of a system with widely separated galaxies (class 0 and 1 in Haan et al. 2013). We cannot see any trend indicating that those objects do not follow the general law we found. This is important since, as noted by Haan et al. (2013), the fact that there is no obvious merging situation does not prevent a cusp from being present.

In order to compare the characteristics of the NGC 1068 cusp with the derived law, we have to estimate the NIR luminosity of the cusp. Because the obscuration is extremely high in the central region, we cannot consider visible data as was done in the previously mentioned studies of cusps in galaxies.

To be consistent, we consider only the H-band data which are less affected by possible dust emission, assuming that it represents well the stellar population distribution. With a calibration done with the star BD-00413, observed the same night, we evaluate a magnitude m_H = 9.5 within an aperture of radius 100 pc. Assuming a stellar population with H-K_s-band ≈0 and an age between 4.0 × 10⁸ and 6.0 × 10⁹ yr, the bolometric luminosity would be between 3.7 × 10⁹ and 1.8 × 10¹⁰ L_⊙, according to Fig. 5 in Thatte et al. (1997). This is consistent with Davies et al. (2007) who derives a bolometric luminosity of stars in the centre of 7 × 10⁹ L_⊙. Using these estimates and the half-light radius 100 pc in Eq. (1), we find a predicted exponent between 1.06 and 1.84, a range that brackets the exponent 1.3 observed on the H-band brightness profile.

This being said, we must emphasise the somewhat large uncertainty linked to the evaluation of the bolometric stellar luminosity. This rather large exponent of 1.3 is also partly consistent with the value presented by Müller Sánchez et al. (2009), who estimated an enclosed stellar mass varying linearly with the radius (M_⋆(r)∝r), thus requiring a stellar density varying as r⁻². The agreement with the empirical law we established suggests that the cusp in NGC 1068 could be closely related to some starburst activity. For Davies et al. (2007) there are indeed strong indications, based on Brγ equivalent width, supernova rate, and mass-to-light ratio, that the stellar nucleus of NGC 1068 is the remnant of recent but no-longer-active starburst episodes. The authors estimate the age of the last episode to be 200−300 Myr.

The only viable class to which the NGC 1068 cusp appears to be related is the CSC. We note that Mushotzky et al. (2014), using Herschel data on a sample of X-ray selected AGNs, reach the conclusion that most AGN live in high star formation nuclear cusps. However one question remains: if we see mainly stellar emission, why are the exponents and half-light radii so different at H-band and K_s-band? A first tentative explanation is that this difference is due to the effect of extinction by dust which may play an important role in this region. Another possibility is that there is a segregation in stellar population, with redder giant stars dominating towards the most central part. A combination of the two effects is also possible. We examine these two cases in the next subsections.

3.3. Pure differential dust extinction

The larger dust absorption at H band compared to K_s band must smooth the stellar brightness more strongly as we approach the central peak, provided that the dust quantity increases inward. To test this effect, we developed a rather simple radiative transfer model, assuming that the ISM in which the cuspy distribution of stars is embedded also follows a cuspy, i.e. power-law, distribution, but with possibly a different exponent. Figure 6 depicts the geometry used.

Fig. 6. Sketch of the simple radiative transfer model we built in order to simulate the differential H-/Ks-band extinction effect of a cuspy distribution of dust in which a cuspy distribution of stars is embedded.

At a given wavelength, the brightness at distance r from the central peak is given by

$$\begin{array}{c}\hfill F(r)\propto {{\displaystyle \int }}_{-h}^h{L}_{\star }(\lambda )\times n_{\star }(\sqrt{(}{r}^2+x^2))\times exp(-\tau (x))\mathrm{d}x,\end{array}$$(2)

where n_⋆(r) is the radial density of star, L_⋆(λ) the star luminosity at λ, and τ(x) the optical depth between the observer and the cell at location x on the line of sight:

$$\begin{array}{c}\hfill \tau (x)={{\displaystyle \int }}_{-h}^x{\alpha }_{\mathrm{dust}}(\sqrt{(}{r}^2+x^2))\mathrm{d}x,\end{array}$$(3)

with α_dust(r) the extinction coefficient of dust per unit of length at radius r. Indeed, α_dust depends on the wavelength, being proportional to Q_ext(λ).

We developed a numerical code in order to look for the exponents of the cuspy distribution of dust and of the star density that would give the best fit to the observed H-band and K_s-band radial distributions. Using the Powell minimisation algorithm and leaving free both exponents, we minimised the square of the log of the ratio of the observed to simulated profiles at H-band and K_s-band simultaneously. We set τ_Ks/τ_H = 0.61, i.e. the value found in the standard ISM. We reached a solution which is rather satisfactory, as shown in Fig. 7 where the radial fluxes at K_s and H band are plotted, modelled using a power law for the stellar radial density of exponent −4.2 and a dust density varying as r^−1.05. In Fig. 8 we plot the deduced radial distribution of the optical depth at H and K_s band down to the centre of the cusp: it appears cuspy as well. From this plot we can derive a radial optical depth τ_H = 7.7 and τ_Ks = 4.7 between 2.5 pc and 60 pc; this corresponds to τ_V = 43 in V.

Fig. 7. Red and blue solid lines: observed 2 to 50 pc radial fluxes in Ks and H band; red and blue dashed lines: Ks-band and H-band profiles resulting from our radiative transfer model, when the power law for the stellar and dust radial density have exponents of −4.20 and −1.05, respectively. The result is a quasi-power law for both brightness profiles.

Fig. 8. Radial distribution of the optical depth in H band (blue) and Ks band (red) in the model providing the best fit to the H-band and Ks-band observed brightness profiles (see text).

This value would be close to twice that derived by Burtscher et al. (2015, 2016) who estimate A_v = 24.6. If our model is correct, the need for some extra extinction arising in the so-called torus would not be required, and the dusty cusp would even have to be truncated at a radius of several pc; in other words, the torus could be identified with the tip of the dusty cusp. On the other hand, Grosset et al. (2018) used their Monte Carlo radiative transfer code Montagn (Grosset et al. 2016) to interpret AO-assisted NIR polarisation measurements (Gratadour et al. 2015) and to estimate a minimum optical depth of 20 at K_s band to the CE (see also next section). In this case, the dusty cusp we suppose here would not be sufficient to explain the most central extinction, and a torus-like extra structure is still needed. Nevertheless, it is also possible that dichroism is also involved in the polarisation seen at the centre, in which case a lower extinction would be required (Grosset et al., in prep.).

If indeed the exponent we derive for the dust and stellar cusps are real, they need to be explained on some theoretical basis. In particular, the exponent of −4 for the stellar density appears extremely steep, and is predicted in none of the theoretical studies we are aware of. We explore in the next subsection the possible role of a segregation of the stellar population with radius, where the intrinsic H-K_s colour would increase when approaching the centre.

3.4. Combination of differential dust extinction and stellar type segregation

We introduced in our numerical model the new constraint that the quantity H-K_s varies as the Log of the radius: H-K_s = C × log10(r/r_max)/log10(r_in/r_max), assuming thus that H-K_s = 0 at r_max and H-K_s = C at r_in. This means that the un-reddened ratio Flux_Ks/Flux_H varies as a power law versus r as r^{C/2.5/log10(r_in/r_max)}. For instance, if r_in = 1 pc and r_max = 100 pc, the exponent of the power law would be −C/5 On the other hand, the maximum value of C, assuming that the very central region is dominated by red giants is 0.3 (e.g. Ducati et al. 2001). We cannot thus expect a power law much steeper than r^−.06, and the effect on the final fit, after the introduction of this new ingredient in the model, cannot be to change significantly the exponent of the stellar cusp nor of the dusty cusp.

This was indeed confirmed by performing the fit that leads to the values of exponents −3.75 and −1.27 for respectively the stellar and the dusty cusps. The stellar cusp exponent is now somewhat less extreme than the value −4.05 we previously determined, but it still remains very large compared to any prediction made by numerical simulations or even observed in the well-documented cusp of the Galactic Centre (Genzel et al. 2003; Schodel et al. 2017) where the exponent is −2.0 for 0.4 pc < r <  4 pc. We note, however, that the mass of the black hole SgrA* is significantly smaller than in NGC 1068.

3.5. Colours and luminosity of the stellar cusp

So far our aim has been essentially to reproduce the H-band and K_s-band slopes of the radial profiles without considering the measured fluxes themselves. Let us now look to see if there is some consistency between the estimated fluxes of the stellar cusp and the model we developed. The measured H-band and K_s-band magnitudes of the cusp alone, within a 1$\stackrel{″}{.}$47 diameter aperture (102 pc) centred on the peak, are respectively 9.51 and 8.63. Assuming that we correctly evaluated the central point-source magnitude in K_s band (see Sect. 5), we derive then for the supposed stellar cusp a magnitude Ks_cusp = 8.9. This means an observed colour (H-K_s)_cusp = 0.6. Now, if we compute in our model, for each filter, the ratio of the intrinsic cusp luminosity to the one escaping from the dusty cusp, we can derive the corresponding global (H-K_s) due to reddening. It amounts to 1.6, assuming the second version of our model, i.e. combining differential dust extinction and stellar type segregation. The intrinsic (H-K_s)${}_{\mathrm{cusp}}^0$ would then be −1.0. Even considering the uncertainties and the crudeness of the model, this value appears inconsistent with a stellar population. For instance, for a B1 star (H-K_s) = −0.14 (Koornneef 1983) and a very young stellar cluster (H-K_s) = −0.1–0.2 (Santos et al. 2013). This is a first clue that the hypothesis of pure stellar cusp is not robust enough.

Now the apparent luminosity of the cusp at K_s band amounts to L_Ks = 3.210⁸ L_⊙. On the other hand, our model gives a fraction of escaping stellar light at K_s band of 5.5 %, leading to an intrinsic K_s-band luminosity of 4.8 × 10⁹ L_⊙. This value is clearly too high since the ratio L_bol to L_Ks is at least 20 for a cluster of 5 Gyr (e.g. see Fig. 5 in Thatte et al. 1997), leading to a total intrinsic stellar luminosity of the cusp of > 1 × 10¹¹ L_⊙, i.e. between 1/4 (Raban et al. 2009) and 2/3 (Bock et al. 2000) of the nuclear bolometric luminosity. This is probably the point that invalidates the interpretation in terms of pure stellar cusp: this fraction is much larger than usually accepted, especially if we consider that a cluster 5 Gyr old of this luminosity would have a mass of 2.8 × 10¹¹ M_⊙ according to Thatte et al. (1997; Fig. 5). In addition, the involved quantity of dust and thus of gas would appear unrealistic as well.

4. Warm/hot dust extended emission

We examine now whether there is an alternative interpretation of the power-law behaviour of at least the K_s-band brightness radial profile in terms of thermal emission by warm to hot dust heated by the energetic radiation from the central engine. Of course this interpretation needs to be somewhat consistent with several known properties of the nuclear region.

One first requirement is that the region where the thermal radiation dominates over other sources, especially stellar radiation, must extend at least to 50−60 pc, i.e. where the cuspy distribution is clearly observed. This requirement means that the dust temperature cannot drop rapidly along a radius, or in other words, that the heating radiation from the CE cannot be attenuated on a small scale-length. A very rough estimate of what this requirement means can be made as follows. Let us consider that the power law in r⁻² is verified in the range 6−60 pc, thus that the ratio of thermal emission integrated on the line of sight between 6 and 60 pc should be 1/100. Assuming that the medium is optically thin and homogeneous, the dust luminosity per unit volume between 6 and 60 pc must be in the ratio 1/1000 because the effective path on which significant emission is integrated towards the observer is 10 times longer at 60 pc. Assuming a radius of sublimation R_sub = 2 pc, a temperature of dust sublimation T_sub = 1500 K, and a dust coefficient of absorption Q_abs(λ)∝λ⁻², the dust temperature would vary as T(r) = 1500(r/R_sub)^−2/6 if the heating flux varies as r⁻², i.e. if it is not attenuated. In that case, the ratio of dust luminosity per unit volume between 6 and 60 pc would then be approximatively 1000 at 2.2 μm, as required. The hypothesis that the CE light is not attenuated may seem an unrealistic constraint because we know that, at least on the direct line of sight to the CE, there must be an extinction A_v ≈ 20 attributed to the torus. However, this does not prevent, in the dust-free bicone perpendicular to the torus, the extinction from being low enough, so that a significant fraction of the volume around the CE can be illuminated directly by the energetic photons. An alternative case could be that the medium is sufficiently clumpy so that a fraction of direct light can indeed reach regions distant from the CE and heat the local dust at a high enough temperature. Scattering on the surface of clumps should also help energetic photons to travel those rather large distances. This is, in fact, the situation that we will favour in the following. To assess how realistic this hypothesis is, we built a model of a population of spherical absorbing/scattering clouds of various radii (uniformly distributed between 0 and 1.6 pc) within a sphere of radius 50 pc and launched 10⁶ photons. For a volume filling factor of 2.5% and 300 clouds, the probability for a photon to reach the external radius of the sphere without being absorbed is 42%. Figure 9 illustrates the histogram of radii reached by photons in such a simulation, where an albedo of 0.5 was taken. We note the very flat distribution of radii.

Fig. 9. Histogram of the radius reached by photons launched from the centre of a sphere populated by 300 absorbing clouds totalling a filling factor of 2.5%. The last and highest bin corresponds to escaping photons. There is a very flat distribution of radii.

Assuming thus that dust can be heated significantly at least 50 pc away from the CE, we adapted the numerical code developed for the stellar cusp simulation (see Sect. 3.3) to model the situation described above. The luminosity of the dust per unit volume is simply taken proportional to Planck’s law at the local temperature. We also assume that there could be some extinction on the line of sight. We note that this last ingredient is not in contradiction with the idea that direct light from the CE can reach the dust at any distance because this last assumption stems from the clumpiness of the medium, while the extended dust thermal emission could suffer extinction by clumps along the line of sight.

By adjusting very few parameters, namely a sublimation radius of 1 pc, a constant dust density and τ_Ks = 6.5 from the centre to 50 pc, we are able to fit fairly well the measured K_s-band radial profile, as illustrated in Fig. 10 by the red dashed line. However, despite the quality of the fit, the required K_s-band extinction is fairly large and hardly compatible with our hypothesis of heating by light directly from the CE. Moreover, when examining the case of the H band, assuming that here too dust thermal emission is dominant, we find no solution by our model that is consistent with the parameters found for the K_s-band fit. The most acceptable solution indicated by the dotted blue line in Fig. 10, requires τ_H = 14.5 a value inconsistent with the classical grain extinction curve where τ_Ks/τ_H = 0.61 (against 0.44 here) and, moreover, no fit is found. In particular, the temperature at large distances is too low to produce a significant flux in H.

Fig. 10. Ks-band (red solid) and H-band (blue solid) radial profiles. Red dashed line: fit to the Ks-band profile assuming a maximum τK = 6.5 on the line of sight. Blue dashed line: H-band profile of dust emission with τH = 15 (see text).

The inescapable conclusion is that we do not actually detect dust thermal emission in H band, but only light of a stellar cusp. This leads us to consider finally, rather than the pure K_s-band profile, the same one from which the stellar cusp contribution is subtracted, assuming that the contribution is just proportional to the H-band profile. We used a weighting factor of 0.8 on the subtracted H-band profile. This corresponds to H-K_s = 0.25, a value usually found in globular clusters. We obtain the green profile in Fig. 11. A fairly good fit by our model is found (green dashed line in Fig. 11) now with τ_Ks <  1.5, a range quite consistent with the hypothesis of clumpiness. We actually assume τ_Ks ≪ 1 in the following.

Fig. 11. Ks-band (red solid) and H-band (blue solid) radial profiles; green solid line: Ks-band radial profile subtracted with a stellar cusp contribution proportional to the H-band profile; blue dashed line: estimated stellar cusp at Ks; green dashed line: best fit by the model of thermal dust emission (see text).

Given the assumed low extinction, the K_s-band magnitude of the hot dust can be estimated at 9.55. The spectral luminosity at K_s band is then 5.0 × 10³⁴ W, which can be translated into L_bol = 3.3 × 10³⁷ W, using our model. With the same model, we can relate this luminosity, derived from observation, to the content in warm dust, using some classical assumptions on the dust: MRN (Mathis Rumpl Nordsiek) size distribution as a^−3.5, a_min = 0.5 nm, a_max = 500 nm, Q_abs from Draine (1985). The mass of dust accounting for the observed luminosity, using the temperature distribution of our model, is M_dust = 0.07 M_⊙, a value which definitely corresponds to a very diluted medium. For instance, we derive a typical hydrogen density of about 2000 m⁻³, so that assuming τ_Ks ≪ 1 is thus fully justified.

We think that with this final configuration we have reached a fully consistent picture where the nature of the two cuspy distributions in H band and K_s band are in fact largely different: it is totally dominated by a stellar cusp at H band and, at K_s band, it is progressively dominated by warm to hot dust thermal radiation when going inward, this hot dust belonging to a very diluted medium. We now can understand why the half-light radii at the H-band and K_s-band cusps, correspond to very different scales, as noted previously: the dust thermal emission dominating the central region in K_s is indeed much more compact than the stellar cusp seen in H.

5. Nature of the central point-like source

As shown in Sect. 2, the K_s-band radial profile reveals that there is a quasi-point-like source superimposed on the power-law aisles discussed above. Using the SPHERE calibration provided by ESO and estimating the best PSF fitting to the point source, we find mKs_point = 10.2 ± 0.2. The point-source provides 31% of the K_s-band luminosity measured in the 50 pc diameter central part. This number is significantly smaller than that estimated by Thatte et al. (1997) from their speckle and spectroscopic data, albeit at a lower spatial resolution.

We estimated that the size of the source cannot be larger that 1.7 pc. As proposed by several authors (Thatte et al. 1997; Gratadour et al. 2003), a straightforward interpretation is that the source corresponds to the emission by dust at high temperature, close to sublimation, at the internal walls of the torus. We considered various estimates of this internal radius. Nenkova et al. (2002) propose the relation R_in = 1.2 $L_{12}^{1/2}$ pc, with L₁₂ = 10⁻¹² L_bol/L_⊙, which leads to R_in = 0.9 pc, using L_bol  =  2.1 × 10⁴⁵ erg s⁻¹ (Marconi et al. 2004); with L_bol = 4.2 ×10⁴⁴ erg s⁻¹ as estimated by García-Burillo et al. (2014) we instead find R_in = 0.4 pc. This is to be compared to the maximum radius of 0.8 pc we mentioned in Sect. 2 for the central point source. There is thus a good consistency with the picture suggested by our SPHERE data.

As an additional step, we can compare the flux we measure for the point source to an estimate of the power emitted by the assumed hot dust. At the distance of NGC 1068, K_s band = 10.2 means a K_s-band spectral luminosity of 2.9 × 10³⁴ W μm⁻¹. Assuming a sublimation temperature of 1400 K and a central cavity radius of 0.6 pc with walls optically thick at 2.2 μm, the expected intrinsic K_s-band spectral luminosity from the cavity would be $\approx \pi B_{\lambda }2\pi R_{\mathrm{in}}^2=4.3\times {10}^{37}$ W μm⁻¹ (assuming an opening angle of the bicone of 120°). This means that a 7.9 mag of extinction in K_s band is required to reconcile the two numbers. This is significantly lower than the 20 mag in K_s band that Grosset et al. (2018) have proposed to explain their NIR polarisation data, but somewhat higher than the range 2.5–4.8 of extinction proposed by several authors (Thatte et al. 1997; Efstathiou et al. 1995; Young et al. 1995). Given our rough estimate of the K_s-band luminosity of the very host dust, which is based on an approximate R_in and geometry, neglects any radiative transfer effect, including scattering, and does not consider a realistic dust composition, the uncertainty is probably larger than 2.5 mag. In any case, the essential of the obscuration must occur in the very central part, thus supporting the hypothesis of a torus with a dense internal shell at a pc scale. This does not prevent a more extended molecular extension of the torus from existing on a 7−10 pc scale, like the one identified by García-Burillo et al. (2016), and even a more diluted envelope on a 30 pc scale, as proposed by Gratadour et al. (2015) in order to explain the peculiar NIR polarisation pattern.

6. Conclusion and prospectives

The high angular resolution NIR images of the central 200 pc of of NGC 1068 reveal a cuspy distribution of the brightness at all wavelengths between H band and K_s band with, however, a continuously increasing absolute value of the exponent with wavelength. Using a simple numerical model, we tried first to interpret this behaviour as resulting from a very cuspy (r⁻⁴) distribution of the stellar density in a central compact cluster down to 1 pc, while the H-/K_s-band cusp exponent difference would be accounted for by a distribution of the ISM, also cuspy but with a less steep exponent ≈ − 1. Introducing some segregation in the stellar population, i.e. assuming that giant stars are dominating towards the very centre, only slightly changes the exponent of the stellar cusp that becomes −3.75, i.e. a value never observed nor predicted on theoretical grounds. The required dust opacity is in all cases very high, so that when computing the intrinsic luminosity of the putative stellar cluster, we reach a value much too high to be realistic.

The alternative interpretation that we explored is that the K_s-band brightness profile is a combination of a stellar cusp and warm to hot dust emission and it appears clearly more acceptable. Considering the residual profile when a stellar cusp component, derived from the H-band profile, is subtracted, a good fit can be obtained by a simple model of dust heated by the CE in a homogeneous and diluted medium.

In addition to the cusp and thermal warm dust radiation, the K_s-band profile features a distinct central point-like source interpreted as the very hot dust close to the central engine. The estimated extinction to the very centre is consistent with previous estimates. The overall picture revealed by this study appears closely consistent with the scheme of a clumpy medium around the CE and with the requirement of a torus on a pc scale. It also brings constraints on the structure of the ISM within the 100 central pc, which must be very diluted.

Finally, we note that the three components we identified within the central 100 pc, namely the stellar cusp, the very hot dust core, and the extended warm dust, each contribute to the K_s-band luminosity for a similar fraction, with magnitudes in K_s-band respectively of 9.8, 10.2, and 9.6.

Several questions are opened by this study, and here we mention a few of them. Can a similar stellar cusp and extended dust emission be observable in other AGNs thanks to HAR imaging or even interferometry? Can radiative transfer models produce a consistent view of all observations, including recent ALMA molecular detection of the torus, mid-IR imaging, mid-IR interferometry, and NIR polarisation? Can the rather steep stellar cusp be predicted by theoretical models combining multi-scale hydrodynamics, stellar dynamics, and star formation? Is the relationship that we established from the data of Haan et al. (2013) between cusp radius, cusp luminosity, and exponent verified for other cases of nearby AGNs? Is there a solid explanation for this phenomenological law, especially that the exponent depends only on the ratio L_cusp/R_cusp? In any case, it is in the views of our team to tackle the two first questions, especially the second by using our numerical model Montagn.

Acknowledgments

The authors would like to acknowledge financial support from the Programme National Hautes Energies (PNHE) and from “Programme National de Cosmologie and Galaxies” (PNCG) funded by CNRS/INSU-IN2P3-INP, CEA, and CNES, France.

References

All Figures

Fig. 1. Ks-band (left) and H-band (right) images of the central 438 × 438 pc of NGC 1068, obtained with SPHERE, where 1 arcsec is 70 pc. The intensity scaling corresponds to an Asinh law. For the Ks-band image the flux at peak (white area) is 2.3 × 10−3 Jy pixel−1, and in the light blue area it is ≈6.0 × 10−6 Jy pixel−1. For the H-band image the flux at peak (white) is 2.0 × 10−4 Jy pixel−1, and in the light blue area it is ≈7.5 × 10−6 Jy pixel−1.

In the text

Fig. 2. Median radial brightness profile around the central peak, measured in the H-band. Blue solid line: Brightness profile in H band up to 430 pc, estimated by performing a median of the pixels on successive circles centred on the central peak. Red dashed line: Power-law profile r−1.2 fitted on the central distribution. Green dashed line: exponential profile fitted on the external observed distribution.

In the text

	Fig. 3. Median radial brightness profile around the central peak. Solid red: measured at Ks; purple: PSF at Ks; red dash: difference between Ks band and weighted PSF; green: truncated power-law profile in r−2; solid blue: H-band profile; blue dash: power-law profile in r−1.2.
In the text

	Fig. 4. Median radial brightness profile around the central peak for the cont–H (blue), H-band (green), cont–K1 (cyan), cont–K2 (magenta), and Ks-band (red) images. The log-log scale shows a power-law behaviour for all four bands. The cont–K1, cont–K2, and Ks-band fluxes have been scaled on the H-band flux at the centre. The dashed line on each curve is the best linear fit for the range 5−60 pc.
In the text

	Fig. 5. Estimated cusp luminosity of a selection of the objects studied by Haan et al. (2013), as predicted by the phenomenological law we found (see text) vs the actual cusp luminosity. The blue line is the linear regression. The cross indicates the objects that are not mergers or are in a wide system.
In the text

	Fig. 6. Sketch of the simple radiative transfer model we built in order to simulate the differential H-/Ks-band extinction effect of a cuspy distribution of dust in which a cuspy distribution of stars is embedded.
In the text

	Fig. 7. Red and blue solid lines: observed 2 to 50 pc radial fluxes in Ks and H band; red and blue dashed lines: Ks-band and H-band profiles resulting from our radiative transfer model, when the power law for the stellar and dust radial density have exponents of −4.20 and −1.05, respectively. The result is a quasi-power law for both brightness profiles.
In the text

	Fig. 8. Radial distribution of the optical depth in H band (blue) and Ks band (red) in the model providing the best fit to the H-band and Ks-band observed brightness profiles (see text).
In the text

	Fig. 9. Histogram of the radius reached by photons launched from the centre of a sphere populated by 300 absorbing clouds totalling a filling factor of 2.5%. The last and highest bin corresponds to escaping photons. There is a very flat distribution of radii.
In the text

	Fig. 10. Ks-band (red solid) and H-band (blue solid) radial profiles. Red dashed line: fit to the Ks-band profile assuming a maximum τK = 6.5 on the line of sight. Blue dashed line: H-band profile of dust emission with τH = 15 (see text).
In the text

	Fig. 11. Ks-band (red solid) and H-band (blue solid) radial profiles; green solid line: Ks-band radial profile subtracted with a stellar cusp contribution proportional to the H-band profile; blue dashed line: estimated stellar cusp at Ks; green dashed line: best fit by the model of thermal dust emission (see text).
In the text