Stellar cusp and warm dust at the heart of NGC1068

Establishing precisely how stars and interstellar medium distribute within the central 100 pc area around an AGN, down to the pc scale, is key for understanding how the very final transfer of matter from kpc scale to the sub-parsec size of the accretion disc is achieved. Using AO-assisted (SPHERE-VLT) near-IR images in H and Ks and narrow-band of the Seyfert 2 galaxy NGC1068 we analyse the radial distribution of brightness in the central r<100 pc area down to the pc scale. The median-averaged radial profiles are adjusted by a cusp (power-law) plus a central point-source. A simple radiative transfer model is used to interpret the data. We find that the fit of profiles beyond 10pc is done quite precisely at Ks by a cusp of exponent -2.0 plus a central point-source and by a cusp of exponent -1.2 at H. The difference between H and Ks can be explained by differential extinction, provided that the distribution of dust is itself cuspy as r^-1. However, the required stellar density follows a r^-4 cusp, much steeper than any other cusp theoretically predicted and the mass of the cluster is unreasonable, even introducing a segregation in the stellar population with an excess of giant stars inward. A much more acceptable solution is found with a K profile dominated by warm dust emission, while the H profile corresponds to a stellar cusp. NGC1068 is shown to satisfy a relationship between half-light radius, cusp luminosity and exponent which suggests that the cusp is the remnant of a recent starbust. We identify the central point-like source with the very hot dust at the internal wall of the putative torus and derive an intrinsic luminosity that requires an overall extinction AK ~ 8, a value consistent with predictions by several models.


Introduction
In an Active Galactic Nucleus (AGN hereafter), considered here as the very central region encompassing the stellar cluster and a dense molecular core, the feeding of the Central Engine (CE in the following), i.e. the accretion disk, requires a matter accretion rate from 0.3 in Seyfert (Bian & Zhao 2003) up to 10 M ⊙ yr −1 (Hopkins et al. 2012) for the most luminous QSO. How is achieved the very final transfer of matter 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 scale below 100 pc is not well established (see Goodman (2003) for a review). Neither is clear any 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. 2016a; Gratadour et al. 2015). One of the most difficult problems for hydrodynamical models is to take into account a huge range of scales, from the 10 kpc one of spiral arms down to the 0.1 pc of the accretion disc, while integrating all possible processes such as cooling, star formation, self-gravity of the gas which are all highly non-linear.
Recently Hopkins et al. (2012), using results from their multi-scale hydrodynamical model (Hopkins & Quataert 2010a,b, 2011, proposed that the formation, through instabilities, of lopsided eccentric gas+stellar discs of size of about 10 pc 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 predicts also that surface densities of gas and stars both follow cuspy distributions. Performing observations in vis-IR that could test the validity of this proposal, especially evidencing the stellar disk, is not obvious because of several limits: i) we are dealing with scales which are largely subarcsec (1pc = 20 mas at 10 Mpc) ; ii) the central region is very obscured by the huge concentration of dust ; iii) the bright dazzling AGN is totally dominating the flux, making difficult the detection, in its immediate neighbourhood, of fainter structures associated with warm dust emission or stellar population. 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 on large telescopes and imaging in the near-IR gives 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 (AO) images of NGC 1068 in K s -band and H, as well as in several narrow-band filters, obtained with SPHERE on the VLT. NGC 1068 is a Seyfert 2 galaxy which has been the target for many studies about AGN since it is one of the closest active nucleus (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 to obtain high angular resolution images in the near-Infrared (NIR), as e.g. 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) program in December 2014. The H-band 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. Even more important, 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, still 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 section 2; the discussion on the power-law radial profiles at H-band and K s -band and their interpretation in terms of a reddened cuspy stellar distributions through a simple transfer model is done in section 3, after some discussion about various classes of cusps observed in center of galaxies. In section 4, we discuss an alternative and more acceptable model where warm/hot dust heated by the central engine is mixed to stellar emission and finally, we examine the central pointlike source in section 5 as well as the constraints brought on the dust opacity by its luminosity, before concluding in section 6.

Data set and data processing
The K s -band 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). 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 2 , or 438 × 438 pc 2 , and the angular resolution, as measured on the PSF, is 58 mas and 53 mas, respectively in K s and in H. with SPHERE. 1 arcsec is 70 pc. The intensity scaling corresponds to a Asinh law. For the K s -band image, the flux at peak (white area) is 2.3 10 −3 Jy/pixel and in the light blue area it is ≈ 6.010 −6 Jy/pixel. For the H-band image, the flux at peak (white) is 2.010 −4 Jy/pixel and in the light blue area it is ≈ 7.510 −6 Jy/pixel.
On the H-band image one can note the elongated disc, slightly asymmetric with respect to the CE, well defined by a strong peak. The disk 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 disk, the barycenter of this contour is offset by 6 pixels ( 5.1 pc) to the West and 1 pixel ( .9 pc) to the North with respect to the central peak.
As regards the K s -band image, the asymmetry is less pronounced and, for instance, the center of gravity of the contour delineating about the same area is only shifted of 2.5 pixels (2.1 pc) and 1.5 pixels (1.3 pc) respectively to West and North. Now, an important question is: what is the source of emission at H-band and K s -band? The answer is not trivial, especially concerning the K s -band. Indeed 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 near-IR spectrum when approaching the central peak is well explained by such very hot dust. In the following, our main concern will be to discriminate between the respective contribution of the two sources of radiation.
Let's first examine the medium scale. We present 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). On the same figure, two fits are represented: one as a power-law for the most internal part (r ≤ 110pc) and an exponential one for the external part, up to r = 430 pc. Both fits are fairly good and, assuming that's essentially the stellar component which is seen at H, one concludes that when going inward, a disk distribution (here of scale length R c = 410pc) is superseded by a cuspy (i.e. power-law) distribution of exponent 1.2 for r ≤ 110pc. In Fig. 3 we plot the median radial brightness profile around the central peak, measured on 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 from the one of the nucleus of NGC 1068, so that the AO performances can be considered as very similar. One can note several features.
First the H-band profile is almost perfectly fitted by a power-law of exponent -1.2; we will come back in the next section on the meaning of this behaviour. Second, we note that the K s -band profile aisles, beyond a radius of about 17 pc, is 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 a this power-law profile truncated at r < 4pc 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 well appearing 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 have 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 -1pc. : PSF at K s ; red dash: difference between K s -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 .
We will come back in section 5 on the interpretation of this important result that there is a quasiunresolved 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.
The data processing is essentially the same as the one used for the H-band and K s observations of Dec. 2014.
We plot on Fig 4 the median averaged profiles at the four wavelengths on a log-log scale. One notices 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 fitted power-laws for the range 5 to 60 pc have exponents -1.20, -1.36, -1.95, -2.11 and -2.03, for respectively H-band, cont-H, cont-K1, cont-K2 and K s -band.
There are several possible origins for the power-law emission, and we'll examine different cases in the next sections, however, the monotonic variation of the slope with respect to wavelength on a sample of 5 different filters is a strong indication that the cause of this variation is rather a continuously variable phenomenon than some spectral features such as emission lines. With this hypothesis in mind and taking into account the fact that the SNR is significantly lower in narrowband images, we have chosen in the following to consider only the H-band and the K s -band bands in our analysis, so that thanks to the broad range in wavelength coverage, any effect depending on wavelength will 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 corresponds to a situation often observed in central region of galaxies (Böker 2008), is that in both bands we sense the emission from a stellar cusp, also called 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 Ks. Of course, some mixing of this dust thermal radiation with the photospheric emission from a stellar cluster is likely, especially at H.

A reddened stellar cusp ?
In first instance, 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) have 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".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. Note that the angular resolution they reached is only 0".9, thus far below the SPHERE performance of 0".06 in K s -band (see sect. ).
Concerning 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".4 resolution speckle observations estimated that essentially all the flux at K s -band could be attributed to the point-like source. This is however not what we observe on the K s -band profile of Fig. 1, where the aisles extend far beyond those of the PSF. On the other hand, the fact that the aisle beyond 17 pc (0".24) are very well fitted too by a power-law 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: why are the exponents 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 disks 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 the K s -band cusps respectively, thus two very different scales, a fact that we should try to explain.
Compared to the structures observed at millimetric wavelengths, either in molecular lines (CO, HCN) or in the continuum (García-Burillo et al. 2014, 2016bImanishi et al. 2018), we note that the extent of the cusp at H-band reaches the inner radius of the off-centered circumnuclear molecular disk. This means that we do not expect significant extra extinction from the disk itself, since the stellar cusp is almost entirely within its inner contour.

Emission lines from the ISM
We'll consider in section 3.3 the possible role of differential dust extinction to explain the log-log H-band vs K s -band slope difference, but before, let's first examine another potential explanation, if there is an extra contribution to the K s -band flux of lines emitted by the ISM, that would be both significant and featuring some gradient. One can think of two kinds of line emission: either the emission of the (1-0) S(1) ro-vib line (2.1218 µm) of some excited H2 , or the Br γ line (2.1655 µm) associated to 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 spectroimager SINFONI-VLT clearly shows that the line emission is confined in 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" centered 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 band is only 5. 10 −4 . The same spectrum indicates that the Br γ line contribution to the same band is approximately 3. 10 −3 and thus should also be excluded as an explanation of the differences in profile slopes between H-band and Ks s -band. A 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).

The 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 more explicit this variety.
We can identify a first category, where the cusp is associated to some compact central stellar cluster. For instance, it corresponds to some extra-light compared to a Sérsic distribution extrapolated to the very center 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 are representing the same kind of object, namely a central compact stellar cluster analog 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 to 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. Those typical characteristics do not depend strongly on the type of galaxy, considering for instance the study by Böker et al. (2004) who find, for a sample of late-type galaxies, an histogram of radii with a first peak at 3 pc and an another one at 8 pc, the range extending to 30 pc. However, it is well 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 some 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 on-going starburst. We refer in particular to the study by Haan et al. (2013) of a large sample of LIRGs and ULIRGs where HST images in near-IR were used to derive the characteristics of the nuclear stellar cusps they found in nearly all objects. We have analysed the data compiled by those authors concerning ≈ 80 LIRGs and ULIRGs objects, a large fraction of which being mergers at various phases of the process. In most objects the authors found a cuspy profile in the nuclear region. They provide a table in which are compiled various properties for the selected objects, including radius and near-IR 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 literature corresponds to the compact stellar population which is distributed around a massive black hole in the center of a galaxy. Some theoretical works have explored the dynamics of this peculiar system (see e.g. a review by Amaro-Seoane et al. (2010)) but there are no well established observational facts, with the notable exception of the center of the Milky Way (Genzel et al. 2003;Schodel et al. 2017). Those NSC (Nuclear Stellar Clusters) 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. Especially, 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 of typically 1.3 (Freitag et al. 2006). Such exponents are thus clearly larger than the one found in the two previous classes we examined, and closer from the ones observed in NGC1068, but if we consider both the predicted small half-mass radius and the too large relaxation time T r required to reach to a cuspy distribution -due to the large mass of the BH -, we'll no longer consider this case for NGC 1068.
To make the distinction between the two classes more readable in the following, we will call Nuclear Cluster Cusp (NCC) the first one, Central Starburst Cusp (CSC) the second one. Our goal is to try to find if there are clues allowing to identify the cusp we observe with one of those classes, or a combination.
The rather large exponent we found for the cusp in NGC 1068, as well as the rather wide half-light radius are clearly ruling out the case where we would observe a mere NCC: such large exponents are never observed for this class.
As regards 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'll see in section 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 some closer examination.
To make a step further, starting from the data set compiled by Haan et al. (2013), we have searched if some general law could relate the parameters of the cusp in CSC. The purpose is to see if NGC 1068 could obey such a general law, even if it is not considered as a LIRGs. In case where the answer is yes, this would be a strong indication that the cusp is rather the result of a starburst evolution than 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, near-IR 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 near-IR cusp luminosity L cusp and the log of the cusp radius R cusp . We found that the law relating the three parameters was rather robust, with a covariance between the actual value and the predicted one of .85. The robustness of the fit can be appreciated on Fig. 5, where the luminosity of each selected object is plotted against the one predicted by the law using the two other parameters.
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, to a lower luminosity corresponds a smoother cusp. Note that on Fig. 5, we have marked with an additional cross the objects which are not belonging to a merger or 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 to be present.
In order to compare the characteristics of the NGC 1068 cusp with the derived law, we have to estimate the near-IR luminosity of the cusp. Because the obscuration is extremely high in the central region, we cannot consider visible data as it was done in the previously mentioned studies of cusps in galaxies.
To be consistent, we consider only the H-band data which are the 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 The only viable class to which the NGC 1068 cusp appears to be related to 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 active galactic nuclei live in high star formation nuclear cusps. It remains however a question: if we do 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 stars (giant) dominating towards the most central part. A combination of both effects is, as well, possible. We examine those two cases in the next subsections.

Pure differential dust extinction
The larger dust absorption at H-band compared to K s -band must smooth more strongly the stellar brightness as we approach the central peak, provided that the dust quantity increases inward.
To test such an effect, we developed a rather simple radiative transfer model, assuming that the ISM in which is embedded the cuspy distribution of stars follows as well a cuspy, i.e. power-law, distribution, but with possibly a different exponent. Fig. 6 depicts the geometry used. At a given wavelength, the brightness at distance r from the central peak is given by 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.
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 stars density that would give the best fit to the observed H-band and K s -band radial distributions. Using the Powell IDL procedure, and leaving free both exponents, we minimised the square of the log of the ratio of observed to simulated profiles, at H-band and K s -band simultaneously. We set τ K s /τ H = 0.61, i.e. the value found in the standard ISM. We reach a solution which is rather satisfactory, as shown on Fig. 7 where is plotted the radial fluxes at K s -band and H, modeled using a power-law for the stellar radial density of exponent -4.2 and a dust density varying as r −1.05 .
On Fig. 8, we plot the deduced radial distribution of the optical depth at H-band and K s -band down to the centre of the cusp: it appears indeed cuspy, as well. From this plot, we can derive a radial optical depth τ H = 7.7 and τ K s = 4.7 between 2.5 pc and 60 pc; this corresponds to τ V = 43 in V.
This value would be close to twice the one derived by Burtscher et al. (2016) and Burtscher et al. (2015) 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 have to be even truncated at a radius of several pc or, another way to tell it, the torus could be identified with the tip of the dusty cusp. On the other hand, Grosset et al. (2018) using their Monte-carlo radiative transfer code Montagn (Grosset et al. 2016) to interpret AO-assisted near-IR polarisation measurements (Gratadour et al. 2015), estimate a minimum optical depth of 20 at K s -band to the CE (see also next section). In that 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 involved too in the polarisation seen at the centre, in which case a lower extinction would be required (Grosset, 2018, in preparation).
If indeed the exponent we derive for the dust and stellar cusps are real, they require to be explained on some theoretical basis. Especially, the exponent of -4 for the stellar density appears extremely steep and 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 center. The result is a quasi-power-law for both brightness profiles.

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-redenned ratio Flux K s /Flux H varies as a power-law vs 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)). One cannot thus expect a power-law much steeper than r −.06 and the effect on the final fit, after 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. Of course the stellar cusp exponent is now somewhat less extreme than the value -4.05 we previously determined, but remains still very large compared to any prediction done by numerical simulations or even observed in the well evidenced cusp of the Galactic Center (Genzel et al. 2003;Schodel et al. 2017) where the exponent is -2.0 for 0.4 pc < r <4 pc. We however note that the mass of the black hole SgrA* is significantly smaller than in NGC 1068.

Colors and luminosity of the stellar cusp
Up to now, 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's now look 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 solely, within a 1."47 diameter aperture (102 pc) centered on the peak, are respectively (9.51, 8.63) Assuming that we correctly evaluated the central point-source magnitude in K s -band (see below section 5), we derive then for the supposed stellar cusp a magnitude Ks cusp = 8.9. This means an observed color (H-Ks) 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 reddenning. 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-Ks) 0 cusp 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.
Concerning now the apparent luminosity of the cusp at K s -band it amounts to L K s = 3.210 8 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.810 9 L ⊙ . This is clearly too high a value, since the ratio L bol to L K s is at least 20 for a cluster of 5 Gyr (e.g. see Fig. 5 of Thatte et al. (1997)), leading to a total intrinsic stellar luminosity of the cusp of > 1. 10 11 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 which 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 11 M ⊙ according to Thatte et al. (1997) (Fig. 5). In addition, the involved quantity of dust and thus of gas would appear as well unrealistic.

Warm/hot dust extended emission
We examine now if 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 requires 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 terms, 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 done as follows.
Let's consider that the power-law in r −2 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 (λ) ∝ λ −2 , the dust temperature would vary as T (r) = 1500(r/R sub ) −2/6 if the heating flux varies as r −2 , i.e. 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 look 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 that in the dust-free bicone, perpendicular to the torus, the extinction can be low enough in a significant fraction of the volume around the CE which can thus 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 for energetic photons to travel up to 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 populations 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 6 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 %, Fig. 9 illustrates the histogram of radii reached by photons in such a simulation, where an albedo of 0.5 was taken. One notes the 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 section 3.3) to model the situation described above. The luminosity of the dust per unit volume is simply taken proportional to the Planck's law at the local temperature. We also assume that there could be some extinction on the line of sight. 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 By adjusting very few parameters, namely a sublimation radius of 1pc, a constant dust density and τ K s = 6.5 from the center to 50 pc, we are able to fit fairly well the measured K s -band radial profile, as illustrated on 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. More over, when examining the case of the H-band, assuming that, here too, dust thermal emission is dominant, we find in fact no solution by our model consistent with the parameters found for the K s -band fit. The least worse, indicated by the dotted blue line on Fig. 10 requires τ H = 14.5 a value unconsistent with the classical grain extinction curve where τ K s /τ H = 0.61 (against 0.44 here) and moreover, no fit is found. Especially, the temperature at large distance is too low to produce a significant flux in H. The unescapable 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 later is just proportional to the H-band profile. We used a weighting factor 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 on Fig. 11. A pretty good fit by our model is found (green dashed line on Fig. 11) with now τ K s < 1.5, a range quite consistent with the hypothesis of clumpiness. We will actually assume τ K s << 1 in the following. Given the assumed low extinction, the K s -band magnitude of the hot dust can be estimated to be 9.55. The spectral luminosity at K s -band is then 5.0 10 34 W, that can be translated into L bol = 3.3 10 37 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 size distribution as a −3.5 , a min = 0.5nm, a max = 500nm, 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 −3 , so that assuming τ K s << 1 is thus fully justified.
We think that with this final configuration we 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 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.

The nature of the central point like source
As shown in section3, the K s -band radial profile reveals that there is a quasi-pointlike 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 mK s 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 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 1/2 12 pc, with L 12 = 10 −12 L bol /L ⊙ which leads to R in = 0.9 pc, using L bol = 2.1 10 45 erg/s (Marconi et al. 2004); with L bol = 4.2 10 44 erg/s as estimated by García-Burillo et al. (2014) we rather find R in = 0.4 pc. This is to be compared to the maximum radius of 0.8 pc we mentioned in section 3 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 some 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 34 W µm −1 . 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 ≈ π B λ 2π R 2 in = 4.3 10 37 W µm −1 (assuming an opening angle of the bicone of 120 • ). This means that 7.9 magnitude 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) propose to explain their near-IR polarisation data, but somewhat larger 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, that is based on an approximate R in and geometry, neglects any radiative transfer effect, including scattering, nor considers a realistic dust composition, the uncertainty is probably larger than 2.5 magnitudes. In any cases, 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, that a more extended molecular extension of the torus exists on a 7-10 pc scale, as the one identified by García-Burillo et al. (2016a), and even a more diluted envelope on a 30 pc scale, as proposed by Gratadour et al. (2015) to explain the peculiar near-IR polarisation pattern.

Conclusion and Prospectives
The high angular resolution near-IR 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 −4 ) distribution of the stellar density in a central compact cluster down to 1 pc, while the H/K s -band cusp's exponent difference would be accounted for by a distribution of the ISM as well 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 center, changes only slightly 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 center is consistent with previous estimates. The overall picture revealed by this study appears well consistent with the scheme of a clumpy medium around the CE and with the requirement of a torus at a pc scale. It also brings constraints on the structure of the ISM within the 100 central pc, which must be very diluted.
As a final remark, 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, contribute each 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, let us mention at least 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 near-IR 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.