© The Authors 2023

1 Introduction

For most of observational astronomy history, protoplanetary disks had been spatially unresolved, but this has changed drastically – particularly in the last decade. Nowadays, different facilities provide a wealth of detail by tracing different (vertical) layers of these objects at various wavelengths (e.g., Menu et al. 2015; Ginski et al. 2016; Andrews et al. 2018; Varga et al. 2018; Bertrang & Avenhaus 2018; Perraut et al. 2019; Keppler et al. 2020; Bae et al. 2021). The contribution of different but complementary observations is demonstrated in Fig. 1. From these observations, we have learned that protoplanetary disks are very rich in structure, such as radial gaps, rings, and spiral arms (i.e., substructure; e.g., Benisty et al. 2022; Boccaletti et al. 2020). The origination and diversity of substructure in proto-planetary disks may be attributable to different mechanisms, such as the gravitational interaction of the disk's material – especially forming or embedded planets (e.g., Mordasini et al. 2015), for which compelling evidence exists in some cases (Keppler et al. 2018).

Thermal infrared emission of protoplanetary disks in most cases remains spatially unresolved to 10 m class telescopes, even when they operate at the diffraction limit, because the equilibrium temperature of the disk material is too low to yield significant radiation at, for example, 10 μm on scales that can be resolved (e.g., Menu et al. 2015). For spatially resolving the thermal emission, long-baseline interferometry is therefore typically required (cf. GRAVITY; Gillessen et al. 2006, MATISSE; Lopez et al. 2022). The situation is different for the emission from polycyclic aromatic hydrocarbon (PAH) particles – large molecules consisting of multiple benzene rings with hydrogen atoms around the edges. These particles are so small that the absorption of a single stellar ultraviolet (UV) photon provides sufficient heating for 10 μm emission to occur far into the disk, allowing them to be observed with single-dish telescopes (Geers et al. 2007; Visser et al. 2007). Such nanoparticles experience stochastic heating followed by a quick de-excitation through infrared emission and then remain cold (and invisible) until they absorb another UV photon (Natta & Krügel 1995; Li & Lunine 2003; Siebenmorgen & Krügel 2010). Thus, for PAH emission, the expected emission profile does not follow the (nonlinear) profile resulting from a thermal equilibrium calculation but instead scales directly with the number of UV photons absorbed in a given disk location. It thus approximately follows an inverse square law with respect to the distance to the central star, modulated by the spatial distribution of the PAH molecules as well as the "flaring profile" of the disk (effective surface height vs. radius) and substructure, such as radial gaps (e.g., Doucet et al. 2007).

Contrary to the macroscopic-sized (millimeter-scale) grains that dominate the emission seen in Atacama Large Milimeter Array (ALMA) continuum images (e.g., Andrews et al. 2018), PAH molecules are strongly coupled to the disk gas (Seok & Li 2017). They are presumably present throughout the disk but are excited and therefore observable only in regions directly exposed to stellar UV radiation and thus probe the highest disk layers (Siebenmorgen & Heymann 2012).

Additionally, it is well established that the PAH emission bands show clear variation in wavelength, shape, and both absolute and relative intensities - not only among different sources but also across spatially resolved sources (e.g., Lagage et al. 2006). In the N band, it has been shown that the 7.9 and 8.6 μm bands are tightly correlated with each other, but not with the 11.3 μm band (Maaskant et al. 2014). Theoretical and laboratory work indicates that the PAH charge state predominantly determines these correlations. While neutral PAHs emit strongly in the 11.3 μm band, ionized PAHs emit strongly in the 7.9 and 8.6 μm bands (e.g., Hudgins & Allamandola 2005; Draine et al. 2021). Since the ionization is sensitive to UV fluxes (Maaskant et al. 2014), spatial tracking of the distinctly ionized and neutral PAH features can provide insight into the optical thickness profile of the disk in which they are embedded.

In early observations performed with 4-8 meter-range telescopes operating in seeing-limited mode, some of the largest disks could be spatially resolved (e.g., Fig. 1). Long-slit spectroscopy aided in this, allowing robust relative measurements between the continuum and the PAH features. However, without adaptive optics (AO) correction, the spatial resolution was significantly lower than the diffraction limit. The most prominent example is the disk of HD 97048 (Doucet et al. 2007), which was already highly resolved with Thermal Infrared MultiMode Instrument (TIMMI2) on the European Southern Observatory (ESO) 3.6 m telescope (van Boekel et al. 2004). Furthermore, the seeing-limited VLT Imager and Spectrometer for mid-Infrared (VISIR) was able to secure a high-quality image of it, which showed the disk surface and allowed the vertical disk height versus radius to be directly (geometrically) measured for the first time (Lagage et al. 2006). Follow-up studies (e.g., Pinte et al. 2019) showed that, in some objects, the PAH emission is dominated by the dense disk regions. In contrast, in other cases, the emission arises primarily in much less dense "gap" regions also prevalent in ALMA/Spectro-Polarimetric High-contrast Exo-planet Research (SPHERE) observations. In these observations, the emission is dominated by ionized PAH molecules, whereas the PAHs in the dense regions are mostly neutral (Maaskant et al. 2014).

We report on observations performed on eight Herbig Ae disks with the Very Large Telescope (VLT) spectrometer and imager for the VISIR mid-infrared New Earths in the Alpha Cen Region (NEAR) as a long-slit spectrometer with a diffraction-limited point spread function (PSF). We obtained spatially resolved N-band spectra (8-13 μm) with a spectral resolution of ≈300 and spatial sampling of ≈ 1.1λ/D, which corresponds to several tens of astronomical units for our targets. Using a nonlinear Markov chain Monte Carlo (MCMC) optimization technique, we fit the reduced spectra with a parametric model consisting of a decoupled blackbody component, an additional blackbody component that is coupled to opacity curves of ten silicate species, and four prominent PAH emission features in the N band, roughly centered at 8.6, 11.3, and 12.6 μm, following, for example, van Boekel et al. (2005). We also include the "red" tail of the 7.9 μm emission feature as a power law. Finally, we produce radial profiles of the PAH emission, with which we explore the physical properties of the disks in our sample, such as the projected geometry of their photospheres, the distribution of PAHs of different charge states across the disk, and the dependence thereof on the intrinsic structure. In Sect. 2 we present our target sample; in Sect. 3 we provide details about the performed observations and subsequent data processing; and in Sect. 4 we present our modeling methodology. In Sects. 5 and 6, we present our results and discuss them, respectively, and in Sect. 7 we summarize and conclude.

Fig. 1 Variations in the effective emission height of different tracers for a 3D disk demonstrated for HD 97048: (1) scattered-light as observed with SPHERE (Ginski et al. 2016), (2) PAH emission as observed with VISIR (Lagage et al. 2006), and (3) 13/12CO as observed with ALMA (Walsh et al. 2016). The inset on the bottom right is Fig. 4 from Lagage et al. (2006), demonstrating the 2D spatial distribution of the 8.6 μm PAH emission intensity of HD 97048.

2 Target sample

We selected several nearby disks based on the presence of PAH features and indications of extended emission in previous VISIR observations (Maaskant et al. 2014). These are: AB Aurigae, HD 95881, HD 97048, HD 100453, HD 100546, HD 163296, and HD 169142. In some sources, the continuum emission may be extended (e.g., van Boekel et al. 2004; Lagage et al. 2006). Furthermore, sufficient optical brightness for the NEAR wavefront sensor was required. Of the resulting sample, eight targets were observed during the VISIR-NEAR science verification time runs; all these targets are known to have disk gaps or ring-like disk structures in ALMA continuum data, scattered light observations, or infrared interferometry.

Table 1 lists our target sample and some fundamental properties. In Appendix A, we provide background information on the individual targets. In Fig. 2, we present reduced SPHERE polarized scattered-light images for all targets in our sample, except HD 95881, which was not observed in this mode. In addition, we over-plot the VISIR long-slit width and mean orientation during observation. In Appendix A, we individually introduce each target in our sample.

3 Observations and processing

We used VISIR-NEAR in long-slit spectroscopy configuration, whose spectral resolution is ≈300 and spatial resolution is 07076/pixel. The slit width is 0."4. Our observations were conducted from 15–17 September and 12–19 December 2019. An overview of the observations is listed in Table 2. Observing setup and data processing stages are discussed in Sects. 3.1 and 3.2, respectively.

3.1 Observing setup

The long-slit spectrograph of VISIR covers the full N band in one setting at a spectral resolution of approximately 300 (see the VISIR user manual¹).

Observations with NEAR could be performed only in pupil-tracking mode, with a predetermined instrument rotation angle that cannot be changed (the NEAR Lyot stop and VLT pupil are then aligned, and the AO control matrix has been calculated only for this angle). As a result, the field rotated during observations. This was uncritical for our observations, as they were relatively short, and the field typically rotated by only a few degrees during an observation. Moreover, knowing the field orientation during each integration allowed us to extract subsets of exposures if the amount of field rotation was large enough to cause significant spatial smearing for some targets. For a given target, we achieved this by stacking only the observations for which field rotation caused a spatial offset at the spatial edges (which we limited to 1".5 from the central source on either side of the slit, as for all sources except for HD 97048 - we detected no significant spatially resolved emission with the analysis presented in Sect. 4.1, when applied to the average of all available exposures of each object,) of the spectra by less than the median full-width half maximum (FWHM) of the 7.9 μm PSF (the narrowest PSF of a significant spectral feature) from the mean slit orientation angle.

We used chopping with the deformable secondary mirror at the maximum throw of 5".5. The stability of the PSF provided by NEAR ensured operational consistency near the diffraction limit. Through Gaussian fitting, we measured the average FWHM of the PSF at 7.9 μm (the shortest wavelength in the analyzed spectra) for all observed calibrators (163 exposures) to be ≈1.15 ± 0.04 λ/D. At longer wavelengths, the PSF quality of science target observations cannot be directly measured in general since the disks that dominate the N-band emission of the science targets may be somewhat spatially resolved in the continuum and more so in the PAH bands. Since the NEAR AO wavefront sensor was specifically designed for a planet detection experiment in the α Cen system, it was not optimized for performance on much fainter targets. While our calibrator stars are much brighter than the adopted brightness limit for the AO system (I ≤ 8.5 mag) and photon shot noise is thus not a limiting factor in the AO performance, the science targets' brightnesses are much closer to the magnitude limit and may be in the regime where the limited number of photons starts to degrade the S/N of the wavefront signal and hence the quality of the PSF. Thus, a marginally resolved disk source observed with a perfect PSF may look very similar to an unresolved disk source observed with a slightly poorer PSF – and the two scenarios cannot be easily distinguished².

We assessed the PSF quality at fainter magnitudes with HD 95881 – which is known to have a very compact N-band continuum emission unresolved with VISIR. At a brightness of I ≈ 8.0 mag, the source is indeed near the limit for the AO wave-front sensor. The FWHM of the PSF ranges between ≈ 1.2 λ/D at 7.9 μm and ≈1.05 λ/D at 14 μm. Therefore, we find that the performance of VISIR-NEAR remains close to diffraction-limited over the brightness range of our science targets (the faintest of which is I ≈ 8.1 mag).

Fig. 2 Mean slit orientation (yellow rectangle) overlaid on scattered-light images observed with SPHERE for all targets in our sample except HD 95881.

3.2 Data processing

We implemented spectral calibration for our targets by observing spectrophotometric standard stars (Table 2), expressing an observed spectrum (of either a calibrator or the science target) as (1)

where S (λ) denotes the raw extracted signal from the source, I_* (λ) is the intrinsic spectrum of the observed object, Α(λ) is the atmospheric extinction profile (which is a strong function of, e.g., integrated water-vapor column and airmass), and R(λ) is the spectral instrument response curve. Since Α(λ) and R(λ) are similar for both the calibrator and the science object (up to a difference in observing conditions), and a simulated stellar spectrum exists for the calibrator star, a calibrated science object spectrum is achieved straightforwardly: (2)

where I_science(λ) is the calibrated spectrum of the science object and I_*,int(λ) is the intrinsic spectrum of the calibrator. However, the science object and its calibrator are not observed under identical atmospheric conditions. The resulting difference between the calibrator's atmospheric extinction profile Α(λ) should also be considered. We compute this difference in the following manner: we simulate two telluric absorption profiles using SkyCalc³, Sky Model Calculator by ESO, each with the respective atmospheric conditions measured during the observation of the science object and calibrator (which vary as a function of airmass and integrated water-vapor column). The ratio thereof makes the correction we introduce into Eq. (2) in order to account for the difference between the A (λ) of the science object and the calibrator. We used the spectral templates presented in Cohen et al. (1999) for the intrinsic spectra of the calibrators.

Finally, we superimposed all snapshot exposures of each object (Sect. 3.1) to receive a final calibrated 2D spectrum, where one axis represents the radial dimension (along the long axis of the slit) and the other the spectral dimension. In Fig. 3 we present a calibrated 2D spectrum of HD97048, where spatially resolved emission in the PAH bands can be noticed.

Fig. 3 Calibrated 2D spectrum of HD 97048. The x-axis is the spectral dimension perpendicular to the long axis of the slit, and the y-axis is the spatial axis parallel thereto. The source is positioned in the center of the y-axis. Spatially resolved emission in the PAH bands is highlighted with white arrows. Note the image was square-root-stretched for demonstration purposes.

4 Modeling

4.1 Preliminary excess flux detection in the 8.6 and 11.3 μm PAH bands

We performed a preliminary analysis of the spatially resolved spectra for all sources in our sample in an attempt to detect spatially resolved excess flux in two prominent PAH emission-bearing spectral bands – 8.6 and 11.3 μm (e.g., van Boekel et al. 2005; Lagage et al. 2006; Maaskant et al. 2014) – by defining two spectral windows, each centered on the emission peak of the 8.6, and 11.3 μm PAH templates, respectively, and whose spectral width is the FWHM of a Gaussian fitted to each template. We then define reference (continuum) spectral windows. In the case of the 8.6 μm aperture, the continuum spectral aperture is of similar width. It is located 2 FWHM (usually approximately 0.5 μm) down the red part of the spectrum (since the 7.9 μm PAH feature often dominates the blue part). The procedure is similar in the case of the 11.3 μm aperture, except that we average two reference spectra on either spectral side of the PAH aperture. Each spectral aperture is then summed over the spectral dimension and normalized to obtain a 1D normalized flux profile, representing the spatially resolved emission at these spectral windows. Finally, we subtract the reference flux profile from the PAH aperture. This procedure yields positive residuals if the PAH emission is more spatially extended than the emission at the reference wavelengths. We note that this qualitative analysis remains insensitive to flux differences in the vicinity around the center of the object because the difference between the two bands there is obscured by the fact that the shape of the PSF dominates the central emission. In Fig. 4, we plot an example of this procedure applied to the 8.6 μm PAH band of HD 97048, and in Fig. 5 we plot the result of this analysis on all object and both bands.

Fig. 4 Excess flux analysis of the 8.6 μm PAH band of HD 97048. Top panel: Gaussian fitting of the 8.6 μm PAH template, with a resulting standard deviation of σtemplate ≈ 0.14 μm. The aperture of the PAH spectral window is defined as ±1.17 σtemplate, and the continuum aperture is set to be two FWHMs down the red part of the spectrum (the chosen widths for each are plotted atop the panel). Bottom panel: Fluxes integrated over the continuum and PAH apertures, and the computed excess flux (blue, orange, and green curves, respectively) and their respective ±1σ uncertainties (see Appendix C), indicated by curve widths.

4.2 Fitting spatially resolved spectra

4.2.1 Rationale

Following our preliminary analysis of the spatial extent of PAH emission at the 8.6 μm and 11.3 μm bands (Sect. 4.1), we analyze the spatial dependence of the spectral shape in more detail for selected objects. We extract spectra at different distances from the star in the along-slit direction, which are the rows of our 2D spectra (see Sect. 3). For each object, we consider 2n+1 spectra, where n is the number of spectra sampled on each star side in the along-slit direction, and one additional spectrum samples the on-star location. We then perform a spectral decomposition and fitting analysis of each spectrum. Because the telescope PSF is sampled by multiple pixels, the spectra of adjacent rows are correlated. Since the ≈diffraction-limited PSF width scales ≈linearly with wavelength, the spatial correlation length is longer on the red side of the N band than on the blue side.

In the fitting procedure, we followed the reasoning and methodology of van Boekel et al. (2005). We performed a decomposition of the observed spectra into the opacity curves of five common mineral dust species, with two grain sizes for each component, plus an empirical template for each of the significant PAH bands and a gray-body continuum component. The chosen mineral species represent the most abundant dust species in circumstellar material that exhibit spectral structure in the N band (e.g., Bouwman et al. 2001). The adopted grain sizes of the mineral components have volume equivalent radii of 0.1 μm ("small") and 1.5 μm ("big"). This choice of dust species and grain sizes was found to provide a suitable set of base functions for reproducing the observed range of shapes of the ~ 10 μm silicate features (van Boekel et al. 2005).

van Boekel et al. (2005) demonstrated that such a comparatively simple model is adequate for a limited spectral range (e.g., only the N band, as in the current study); in this case, more sophisticated models such as that by Juhász et al. (2010), containing considerably more degrees of freedom describing the disk structure, are insufficiently constrained due to the small spectral coverage and would constitute an over-fitting of the data.

Unlike the single N-band template for all the PAH emissions, which was used for the modeling of van Boekel et al. (2005), we modeled four PAH emission bands independently rather than as a single template. We generated the PAH templates for each band as follows: 7.9 μm was fitted as a power-law for each source independently since the peak of its emission is not within the modeled spectral range. For 8.6 μm, we generated an individual high-S/N template of this band for the following sources: HD 97048, HD 100453, and HD 169142. For the other sources, we assumed a template that is the average. For 11.3 μm, we repeated a process similar to that of the 8.6 μm band. The 12.6 μm band exhibits the weakest signal, for which it was impossible to generate a satisfactory high-S/N template from our data. Instead, we compared the highest-S/N instance of this emission band (in HD 100453) to two functions with optimized parameters: a Gaussian and a Lorentzian, and used the Kolmogorov-Smirnov (K-S) test to determine which of the two functions resembles it more. We found that the Lorentzian (with optimized parameters df and γ) scored higher on the K-S test and was subsequently adopted as our template for this emission feature.

Fig. 5 Normalized flux excess in the 8.6 μm and 11.3 μm PAH bands (top and bottom panels, respectively). The width of the lines represents the ±1σ residuals of the PAH and continuum windows (see Sect. 4.1). The x-axis is parallel to the slit (i.e., its spatial dimension), and the zero point is the object's center (the star). Note that negative values indicate that the flux in the continuum windows is greater than that in the PAH window. This is caused by silicate emission in the continuum window.

4.2.2 Parametric model

For each model N-band spectrum, we fit 18 free parameters, accounting for the following: (1) a decoupled continuum component; and (2) an additional continuum component that is coupled to ten silicate opacity curves of five silicate species, each of two different grain sizes, and four PAH emission components. The model is described analytically as (3)

where λ is the range of modeled wavelengths (≈8-13 μm), y is the simulated model. The degrees of freedom are color-coded in orange and are as follows: D and C are the amplitudes of the decoupled and coupled continuum components, respectively, T₁ and T₂ are the effective temperatures thereof, w_i is the relative abundance of the ith silicate component opacity curve (S_i), with (Σ_i w_i = 1), and A_j is the amplitude of the jth PAH component (Pj).

We plot the constituents of our parametric model in Appendix B.

4.2.3 Optimization With MultiNest

For each N-band spectrum, we perform nonlinear optimization of a parametric spectral model and derive its constituents' uncertainties using MultiNest (Feroz et al. 2009), a multimodal MCMC algorithm with implementation in Python (Buchner et al. 2014). MultiNest implements nested sampling (Skilling 2004), a Monte Carlo method targeted at the efficient calculation of the Bayesian evidence. Using MultiNest allows for exploring regions of interest in the parameter space with a reduced risk of being trapped in local likelihood maxima.

We assume uniform, hardbound priors for all parameters. Priors of the fractional amplitudes of the silicate opacity curves (w_i) are always normalized by the sum thereof, the priors of each of the PAH template amplitudes (Α_j) are bound to a range of 0–100, and the priors of the effective temperatures of the Planckians (T₁, T₂) span across a range of 50–600 K.

That said, from examining the posterior distributions of our spectral fitting procedure (see Appendix D), we find that degeneracies between any combination of parameters are insignificant, in agreement with van Boekel et al. (2005), who found the silicate components to be largely uncorrelated.

We describe the statistical model used for the MCMC optimization routine in Appendix C.

5 Results

5.1 Fitting unresolved (integrated) spectra

In Fig. 6 we show the calibrated long-slit spectra of our targets, and in Table 3 we list the median values and ±50 ± 34.1 percentiles of each constituent of our parametric model (see Appendix B) for each target. The large central panels show the total spectra (black curves), summed over the spatial dimension. The smaller panels on either side show the spectra at different spatial positions along the spatial dimension of the slit. From bottom to top we show spectra extracted at the central source position, 0."9 away from the center, 2."1 away from the center (all extracted over aperture approximately one FWHM of the PSF in size), and integrated over 2."7 to 3."3 (an aperture approximately three PSFs in size), as labeled in the left (or right) column of the object at the top of the figure. We plot the posterior distributions of the model parameters for each object in Appendix D.

The total emission is strongly dominated by the flux from the central positions. We do see emission at the off-axis positions, in spectral shape often reminiscent of the spectrum at the central location, but with notably bright PAH emission in some of the targets (HD 97048, HD 100546, HD 169142). For the off-axis apertures, the expected temperatures for dust in thermal equilibrium are typically too low to lead to significant thermal emission; the off-axis continuum emission may therefore be dominated by "contamination" from the central source, that is, light emerging from within the central disk regions that is redistributed in the focal plane by the PSF. We investigate this hypothesis in the following.

We first created a PSF model as described in Appendix F. We then estimated the flux from the inner disk regions by performing 1D Gaussian fits to the spatial profile at each wavelength, where the flux in each spectral bin was taken to be the volume of the Gaussian fit. We then assumed this light to come from a spatially unresolved source and convolved it with our PSF model to estimate the contamination from the inner disk regions at all spatial locations, and extracted the expected contamination from the respective apertures; the resulting spectra are shown in blue in Fig. 6. We observe the following.

For most off-axis apertures, the contamination from the central source (blue lines) accounts for much or all of the observed flux, except for the PAH bands. We note that the calibrators are optically much brighter than the science targets. Hence, the AO system always operates in a regime where the PSF quality is not limited by photon noise from the star. The science targets are closer to the limiting magnitude of the AO system (I ≈ 8.5 mag), and therefore they will be in the range where the PSF quality depends on the target brightness in the optical. The PSF quality may be somewhat degraded for the science targets compared to the calibration observations from which the PSF model was created.

For AB Aurigae and HD 100546, we see that in the apertures above and below the central source, there is somewhat more flux than the redistribution model predicts; these systems are dominated by (silicate) dust in thermal equilibrium, with a transition-disk geometry. At larger radii, the observed fluxes approach those predicted by the contamination model again. Thus, the silicate emission appears slightly resolved in the VISIR-NEAR observations.

For HD 97048, the continuum emission in the off-axis apertures is substantially higher than predicted for contamination by the central source, particularly on the southwest side. The large spatial extent in the continuum was already noted by van Boekel et al. (2004); this emission stems from nonthermal equilibrium dust.

The PAH bands in the off-axis apertures for HD 97048 are much stronger than can be explained by contamination from the central source; the off-axis equivalent widths of the PAH bands are much larger than at the central position. For HD 100546, the PAH emission is also clearly stronger off-axis than in the central aperture. For HD 169142, this may be marginally the case.

Finally, for HD 169142, the central peak of the spatial profiles is shaped somewhat differently from the other sources. Whereas in the other sources, the continuum profile is consistent with strongly centrally peaked emission broadened by the PSF, in HD 169142, the observed profile is broader and more triangular. We qualitatively discuss this in Sect. 5.2.5.

Fig. 6 Calibrated long-slit spectra of our targets. Left panel: Spatially resolved spectra. The black curves are down-slit spatially resolved spectra. Each spectrum represents an average of a ≈ λ/D-wide spectral bin (every three rows of the reduced 2D spectrum, whose centroids are separated by 0″.228 and are arranged bottom up, where the bottom spectrum is that of the central object). The uppermost spectrum, farthest from the star, is averaged over three such apertures. The blue curves are the expected flux contamination coming from the unresolved central source, extracted as described in Sect. 5.1. The off-axis direction is indicated by the two letters in the top or bottom panels (e.g., NE for northeast). Middle panel: Fitted integrated spectra. The solid black line is the observed spectrum, the thicker orange line is the best-fit spectrum, and the filled gray area is the ±1σ error-range, where σ is the flux-wise standard deviation of spectra generated with all best-fit parameters within the 50 ± 34 percentiles range. Right panel: Similar to the left panel, but on the opposite side of the star with respect to the slit axis. Note that spectra marked with an asterisk are identical in both the left and right panels. Note also that the spectral artifact visible at ≈9.6 μu is caused by the atmospheric ozone absorption band. Continued on the next page.

5.2 Object-wise spatial intensity distribution of PAHs

We performed a nonlinear optimization analysis (Sect. 4.2) of the spectra of five targets that exhibit spatially resolved PAH emission: AB Aurigae, HD 97048, HD 100546, HD 163296, and HD 169142. The remaining objects were not analyzed for the following reasons: HD 100453 did not display any spatially resolved PAH structure across either band, and HD 95881 and HD 179218 displayed complex spectral structures of silicate emission that could not be accurately modeled with our templates.

In Fig. 7 we present our results of spatially resolved PAH emission of four spectral bands of five select targets. In the following, we discuss the results of each object individually.

5.2.1 AB Aurigae

From Fig. 6, it is evident that PAHs across all bands contribute very little to the overall emission and that silicates (primarily olivine and pyroxene) dominate the emission, in agreement with Bouwman et al. (2000), who used an Infrared Space Observatory (ISO) spectrum to determine the composition of the circum-stellar material. Indeed, the silicates are so dominant in the integrated spectra that the PAHs are hardly seen, except for a slight peak at 11.3 μm. However, the spatially resolved spectra show weak PAH emission (primarily also of the 11.3 μm band) that peaks at some distance from the star on either side.

The PAH intensity (Fig. 7) shows that its distribution across all bands exhibits a ring-like structure peaking relatively uniformly at ±0."3 and ±0."4 at the NE/SW sides, respectively, with varying widths. However, when observing the integrated spectrum (Fig. 6), the marginal flux of the PAH bands cannot be significantly determined due to the much-brighter silicate emission. This is in agreement with van Boekel et al. (2005), who reported that no PAH emission was detected in the integrated N-band spectrum of this object, as was observed with TIMMI2 on the 3.6 m telescope. As we are aware, we present the first-ever detected PAH emission from within this object. Amorphous silicate emission dominates in the central region close to the star, and no PAH emission is detected. Unlike in the case of HD 100546 (see Sect. 5.2.3), we detect no degeneracy between the amorphous silicate templates and the PAH templates (where here the dominant species are olivines).

Fig. 7 Radial profile of the normalized PAH emission in four spectral bands for select sources. Sources are ordered column-wise, and the PAH bands are ordered row-wise and are color-coded (from 7.9 μm in the uppermost row to 12.6 μm in the bottom row). The filled area represents the ±1σ error range around the median value of the emission amplitude posterior distribution of the relevant PAH template (see Appendix C). The orientation of the spatial dimension (parallel to the long axis of the slit) is denoted in the upper-right corner of the uppermost panel. Additionally, the dashed gray line is the high-S/N synthesized PSF (Appendix F), plotted for reference.

5.2.2 HD 97048

As the spectrum of this object is dominated by a featureless continuum with strong PAH emission on top, discerning this emission is not as challenging as it is for objects with substantial silicate emission (e.g., AB Aurigae, HD100546, and HD163296). Here, a few interesting features of the PAH emission are evident: (1) for all three bands of significant S/N (that excludes the 12.6 μm band), the PAH emission exhibits radial asymmetry, consistent with the SPHERE intensity profile and the derived projected geometry of the disk (Lagage et al. 2006; Doucet et al. 2007; van der Plas et al. 2017). (2) the intensity profiles of all bands detected at high S/N (mainly the 8.6 μm and 11.3 μm) exhibit substructure that is reminiscent of that seen in the SPHERE intensity profile, indicating that the PAH emission follows the gas density in the high disk regions as traced by the small dust grains, well coupled to the gas, that is seen in scattered light. (3) when considering Figs. 6 and 7, it is evident that for this object, the 11.3 μm feature becomes relatively strong compared to the 7.9 μm and 8.3 μm bands at larger distances from the star. This, in turn, implies that the fraction of PAHs with neutral charge increases with increasing radial distance; that is, the ionization fraction is the highest close to the star; this is consistent with the observations and analysis by Maaskant et al. (2014).

5.2.3 HD 100546

In this object, we observe apparent differences in the spatial intensity profiles between the different PAH bands. The most prominent example is the 11.3 μm feature, which is centrally peaked whereas the other three bands show a ring-like profile, albeit with different radial decay rates (see Fig. 7). At first sight, this is puzzling, and cannot be explained with, for example, a specific spatial distribution of neutral and ionized PAHs: the emission in both the 11.3 and 12.6 μm band is dominated by neutral PAHs, so they should show a similar behavior if this is the underlying effect governing the observed profiles. However, the 12.6 μm band, as do the 7.9 and 8.6 μm bands, show a ring-like structure.

We suspect that the apparent deviant behavior in the 11.3 μm band occurs for methodological rather than astrophysical reasons. The 11.3 μm band spectrally overlaps with a strong emission feature of small forsterite grains of similar wavelength and spectral width. If the template of forsterite we used in our spectral decomposition is imperfect (i.e., its spectral shape does not exactly match that of the astrophysical grains), the minimization procedure may compensate for this by adding 11.3 μm PAH emission. We have two sources in our sample in which forsterite is abundant, HD 100546 and HD 163296. In both these (and only in these), our spectral decomposition yields a strongly centrally peaked 11.3 μm PAH profile, whereas all other PAH bands show a ring-like structure and are depleted near the central star. While from such a small number of stars, no hard conclusion can be drawn, this observation is consistent with our hypothesis that in disks rich in forsterite, our spectral decomposition mistakes the crystalline silicate emission close to the star for 11.3 μm PAH emission.

In their analysis of properties of forsterite in HD 100546, Mulders et al. (2011) find forsterite to be the most abundant of the crystalline silicates in HD 100546 and find it to be primarily concentrated in the central ~20 au vicinity around the star. In addition, they find a grain-size distribution of 0.1-1.5 μm to fit the observed emission features optimally. More importantly for our discussion, Juhász et al. (2010) and Mulders et al. (2011) show that the central wavelength The emission feature at 11.3 μm exhibits a shift of several percent as a function of grain-size distribution and can thus degenerate with the 11.3 μm PAH feature. They furthermore demonstrate this dependence on the central wavelength 69 μm emission feature (e.g., iron content and temperature); however, in the case of extensive laboratory experiments, Zeidler et al. (2015), who explored the temperature dependence of such emission spectra, demonstrate that they exhibit very little spectral shift for the temperatures that we expect (on the order of 200 K). In the integrated spectrum, we find that the contribution of the 11.3 μm PAH emission band is approximately 2% in the relevant spectral window. Thus, a systematic discrepancy in the fitting procedure of the 11.3 μm PAH feature might have occurred because our forsterite templates fail to reproduce the correct shape of its emission since our model considers only two distinct grain sizes. In contrast, more should be considered for a robust fit (Mulders et al. 2011). This would then lead our fitter to systematically compensate the difference in the spectral shapes of the forsterite emission and our templates by increasing the 11.3 μm PAH emission - which is unlikely to be detected in the posterior distribution analysis. Thus, we decide not to consider the 11.3 μm PAH contribution in the unresolved center of the object as a solid detection.

The three remaining bands do not exhibit such proximity to crystalline emission features and exhibit a relatively consistent asymmetric ring-like structure of the spatial PAH intensities. In Sect. 6.2.3, we compare our derived PAH intensity profiles with that of the SPHERE image to characterize it further. Here, the brighter side in these bands coincides with the far side of the disk, as expected from the larger projected disk surface area visible in the SPHERE images (Fig. 2).

Finally, in their work, Maaskant et al. (2014) find that gas flowing through disk gaps can contribute significantly to the observed ionized PAH emission. This could manifest as an increase in emission in the 7.9 μm and 8.6 μm (approximately ionized) bands, relative to the 11.3 μm and 12.6 μm (~neutral) bands, corresponding to the angular size of a gap. Here, we find that while the two bluer PAH bands exhibit considerably different spatial distribution, they are wider than the 12.6 μm bands -the only one probing the predominantly neutral PAH emission. That said, their respective widths also vary significantly and do not seem to represent the distribution of ionized PAHs identically. A more informative parameterization of the radial behavior of the charge state of PAHs may be through ratios of PAH band intensities, as discussed in Maaskant et al. (2014), rather than the individual intensities themselves. We hope to address the probing of the spatial distribution of the charge state of PAHs in future work.

5.2.4 HD 163296

In this object, we notice a discrepancy similar to that in HD 100546 (Sect. 5.2.3), where the 11.3 μm PAH band profile is centrally peaked, whereas the two asymmetric bluer bands exhibit a ring-like shape. However, the bluer bands in HD 163296 exhibit a roughly similar radial decay. The total flux of the PAH bands is considerably lower (by a factor of ≈4) than that of HD 100546, and the profiles are noisier, especially the 12.6 μm band that is of too low S/N to constrain significantly. We suspect that here, too, the centrally peaked shape of the 11.3 μm PAH profile is caused by a degeneracy with the forsterite templates and does not represent the real intensity profile. In the case of the two bluer bands, the asymmetry in the ring-like lobes of the intensity distribution is in agreement with the projected geometry seen in the SPHERE image (Fig. 2).

This is interesting in light of recent work by Petit dit de la Roche et al. (2021), who had also observed this object with VISIR-NEAR (albeit in the imaging mode) through two different N-band filters, where one corresponds to the 8.6 μm PAH band (PAH1 filter, centered at λ₀ = 8.58 μm with Δλ = 0.41 μm) and the other contains the 11.3 μm PAH band (NEAR filter, centered at λ₀ = 11.25 μm with Δλ = 2.25 μm). They found, however, that in both cases-the bulk of the intensity is contained within an unresolved inner region of the disk. This is probably because of the much brighter silicate emission originating from this unresolved inner region, which also dominates much of the N-band spectral range (see Fig. 6) and is difficult to decouple from the PAH emission without spectral decomposition. In our work, the indication that the silicate emission originates in an unresolved inner disk region is also manifested in the aforementioned degeneracy with the 11.3 μm PAH band. Thus, by probing the individual PAH bands, we can resolve the disk and trace prominent substructure features better spatially.

5.2.5 HD 169142

In most sources, the emission ("along-slit" intensity) profile in the continuum is consistent with strongly centrally peaked emission broadened by the PSF. The observed profile of HD 169142 is broader and more triangular. It suggests overall more spatially extended but less centrally peaked emission, marginally resolved and broadened by our PSF. We qualitatively test this scenario, inspired by the SPHERE image that shows a bright ring of ≈0".717 in radius (Bertrang & Avenhaus 2018), and with a radial width of roughly one-quarter of the radius.

We model the intrinsic emission profile as a simple ring of radius R and a radial width of R/4 of uniform intensity. This mimics the profile observed with SPHERE and is a reasonable choice for the functional form of the intensity distribution. Still, we do not spatially resolve the width of the ring in the VISIR data, and other choices of parametrization may be equally valid. The situation is simplified by the system's nearly perfectly face-on orientation; therefore, we do not incline the ring. We convolve the model image with a model PSF and project it onto the 0."4 VISIR slit. Then, we integrate over the across-slit dimension to obtain a spatial, along-slit profile of our model that can be combined directly with the observed spatial profile. We then optimize the ring radius R to match the observed profile best. This is illustrated in Fig. 8. As a model PSF, we use the following approximation: we take a perfect diffraction-limited PSF of a circular aperture with a diameter of 8.2 m and a central obscuration of 1.2 m, following the geometry of the VLT unit telescopes. We slightly broaden and "soften" the PSF by convolving it with a Gaussian kernel of width σ = 0."075 to account for the PSF being close to but not quite diffraction-limited. The kernel width was chosen so that the synthetic PSF profile (black curve in Fig. 8) best matches the core of the observed PSF profile (red curve in Fig. 8). We see that the central core of the model PSF matches the observed one nearly perfectly (note that the model curve has a much denser spatial sampling), but the in the wings, the observed curve has slightly more power. We use this approximation because we need to convolve the 2D model image with the PSF before projecting it onto the slit. Still, we do not directly measure the PSF in our long-slit spectroscopy data. Since the current discussion is qualitative only, we deem this sufficient.

The ring radius for which the modeled profile best matches the observations is 0."14 (gray curve in Fig. 8); here shown for the 11.3 μm PAH band. The model curve reproduces the observed profile within the observational uncertainties out to a radius of ≈0."35 on either side of the center. In contrast, at larger radii, the observed emission levels are somewhat higher than those given by the model. This is partly due to the true PSF having more power than the PSF model and, in part, to the source showing low-intensity PAH emission beyond the ring in our simplified model.

This experiment shows that the observed PAH emission profile in HD 169142 can indeed have a marginally resolved ring-like shape, plus some extra lower intensity emission at larger radii. There is no dominant contribution from a centrally peaked source. Interestingly, the best-fit ring radius of R = 0."14 is smaller than the one observed in near-infrared scattered light with SPHERE (≈0."17) and also smaller than the continuum size of ≈0."16 estimated from Q-band observations (18.8 and 24.5 μm) by Honda et al. (2012). In Fig. 8, we also show the best-fit profile for a ring with R = 0."17 with the gray, dashed curve. This model yields a worse match and is not consistent with the observations.

Thus, the emerging picture is that the central disk gap in HD 169142 is mostly devoid of PAH emission. The transition from the gap to the optically thick disk region is smooth, with gradually increasing gas density. In such a configuration, the PAH emission traces the lowest density surface layers of the inner wall of the optically thick disk, analogous to the profiles shown in Fig. 1. We expect the PAH emission to trace layers similar to the highly optically thick tracer ¹²CO. Indeed, recent ALMA observations show a ¹²CO profile peaking around 0."12 (Yu et al. 2021), which is at a slightly smaller radius than our best-fit value for the PAH ring.

This picture agrees qualitatively with the study's results by Devinat et al. (2022), who analyzed archival observations of HD 169142 taken with VLT-Nasmyth Adaptive Optics System -Near-Infrared Imager and Spectrograph (NACO)(AO-supported) and VLT-VISIR (seeing-limited but with image quality close to diffraction-limited). Using radiative transfer modeling, they find that the PAH intensity must be strongly suppressed close to the central star. Their observations can both be reproduced with an inner cavity (radius 0."17 = 20 au) that is entirely devoid of PAH emission (i.e., all PAH emission arises from the optically thick disk part beyond 0."17) and with a model where PAHs are still present in the 5–20 au region but at a much-reduced abundance compared to the region beyond 20 au. Our simple model, which has all PAH emission confined to a narrow radial region by design, would favor the latter option since a model in which the PAH intensity inside 20 au is zero would yield a broader profile than observed (see Fig. 8). We also observe PAH emission at much larger radii, which is not included in our simple model, as is evident in Fig. 8 where the measured intensities are higher than the modeled ones.

Finally, Petit dit de la Roche et al. (2021) compared the radial intensity profiles of this object, as observed through the PAH1 and NEAR spectral filters (see Sect. 5.2.4) with ProDiMo-generated models at these bands (Woitke et al. 2009). While the ProDiMo model predicted a roughly similar apparent size for both bands (43 and 45 AU, respectively), they find that the object's apparent size, as observed through the PAH1 filter, is considerably smaller (24 ± 1 AU). This agrees with observations by Okamoto et al. (2017), who attribute >9 μm emissions to the inner wall of the disk, rendering the apparent size larger. They suggest a potential explanation for this discrepancy based on the results of Maaskant et al. (2014) (see Sect. 5.2.3): it is possible that gas flowing through the disk gap could manifest in increased emission at the 8.6 μm PAH band relative to the 12.6 μm PAH band, corresponding to the angular size of the gap. Assuming that the neutral PAH emission is centered in the gap wall, one would expect a more compact apparent size seen through the PAH1 filter instead of the NEAR filter. Since the ProDiMo model does not account for this difference, this leads to the aforementioned discrepancy.

In this work, however, we find that the object's apparent size is identical across all four PAH bands to within 1σ, such that we do not seem to identify the different spatial distributions of PAHs as a function of their charge state. Therefore, we suspect that the effect detected by Petit dit de la Roche et al. (2021) may reflect a difference in continuum emission (as suggested by Okamoto et al. (2017)), implicitly integrated within the filter apertures (especially the wider NEAR filter), rather than a difference in the distribution of ionized versus neutral PAHs.

Fig. 8 Spatial emission profile of HD 169142 in the 11.3 μm PAH band. The red curve shows the observed PSF profile (see Appendix F), and the black curve shows the employed approximation thereof used in the 2D convolution of the model (see Sect. 5.2.5). The solid gray curve shows the ring model that best fits the observed central peak (ring radius 0."14). The dashed gray curve shows the best-fit model with the ring radius fixed to the value observed in the SPHERE data (0."17).

6 Discussion

6.1 Spatial PAH intensity profiles: Ring-like versus centrally peaked

Of the five sources in Fig. 7, three (AB Aurigae, HD 100546, HD 163296) have a well-resolved ring-like geometry in the PAH emission. Interestingly, these all have strong emission features from (sub)micron-sized silicate grains. Two sources (HD 97048, HD 169142) show more centrally peaked PAH emission, and these both do not show silicate emission; in HD 169142, the PAH emission may also have a ring-like geometry but with a much smaller radius than in the "silicate" sources (Sect. 5.2.5). While the number of sources is too small (and their diversity in properties too large) to draw any firm conclusions, we may establish that: (1) the spatial structure of the disks, in particular radial gaps, causes differences in the emergent spectra by excluding specific temperature regimes (see Fig. 7). Our targets have diverse geometries, and all but HD 163296 have large gaps in their disks' inner ≈20 au (see Fig. 2). However, substructure alone cannot account for the observed dichotomy of spatial PAH emission profile regimes, as, for example, HD 97048 and HD 100546 are fairly similar in disk substructure (i.e., both have an inner disk then a large radial gap and an outer disk; e.g., Ginski et al. 2016; Mendigutía et al. 2017, respectively) but exhibit the opposite behavior of the PAH profiles. (2) The correlation between the geometry of the PAH emission and the presence of small silicate grains is suggestive. With that in mind, potential explanations are as follows.

First, small silicate grains are absent in the surface layers of HD 97048 and HD 169142 due to grain growth and settling. This has two observational consequences: (1) spectral features of small silicates naturally disappear, and (2) PAHs can absorb more UV photons for lack of competition from small silicate grains, thus increasing the brightness of PAH emission features (Dullemond & Dominik 2004). In the other three targets, small silicate grains remain in the inner disk regions, absorbing stellar radiation and providing much surface area for PAHs to coagulate onto, thereby dimming PAH emission from the inner disk regions. Why the small silicate grains have disappeared from the inner disk regions in two sources but not in the other three remains an open question.

Second, the composition of the inner disks of HD 97048 and HD 169142 may be intrinsically different (more carbon-rich; Van Kerckhoven & Waelkens 2002; Devinat et al. 2022) from the other objects and as such may contain predominantly carbonaceous dust and little or no silicates. However, the required depletion of oxygen as a leading-order explanation for the lack of silicate emission within the inner regions of disks where we detect a ring-like PAH emission profile appears chemically unlikely. The bulk composition of the disk material is unknown, but the starting material (interstellar medium) is expected to be oxygen-rich (C/O ≈ 0.52 Anders & Grevesse 1989). The photospheres of both stars do not suggest a strong depletion in oxygen in the accreted material (for HD 169142, C/O ≈ 0.68 (Folsom et al. 2012), for HD 97048 we could not find a carbon abundance in the literature, but the oxygen abundance is roughly solar (Acke & Waelkens 2004). If anything, carbon rather than oxygen is expected to be depleted in the solid phase in the hot inner disks, as was strongly the case in the solar nebula.

We note that in the sample of Folsom et al. (2012) HD 169142 has a somewhat higher than average photospheric C/O ratio, but that of AB Aurigae is even higher (approaching unity). That disk has one of the strongest 10 micron silicate features known. In summary, attributing the absence of a silicate feature to an intrinsically oxygen-depleted dust population containing little or no silicates appears currently not supported by observations.

Finally, PAHs may be depleted in the inner disks due to coagulation onto other dust. In their work, Geers et al. (2009) find that the absence of PAH emission can be explained by the depletion of PAHs in the early phases of star formation, by being trapped in ices and or coagulating with dust or other PAHs. This process likely continues well into the class II phase, where an abundance of small dust grains in the disk to coagulate exists. Therefore, one would expect such depletion to manifest in inner regions of disks (where dust densities are highest), where small silicate grain emission is observed.

6.2 Comparison with SPHERE intensities

We then performed a qualitative comparison of the derived PAH intensities for select sources (Sect. 5.2) with a 1D intensity profile of the SPHERE images (Fig. 2). This comparison is worthwhile because both these observations probe the uppermost vertical layer of the disk essentially, albeit different components thereof, which might adhere to different physics (e.g., Lagage et al. 2006; Ginski et al. 2016). Another critical difference between the two is the inherent modulation of the intensity distribution of polarized radiation due to the projected geometry of the disk. Similarly, the PAH intensity profile is expected to scale with the amount of UV radiation reaching the disk surface; this would also yield a 1/r² intensity profile, modulated by disk structure effects and spatial variations in the abundance and molecular properties of the PAHs.

To compare the two, we generated 1D intensity profiles from the SPHERE images as follows: first, we overplotted the proper dimensions and orientation of the VISIR long-slit as illustrated in Fig. 2. We proceeded to average the intensity of all pixels within the slit, parallel to its long-axis, to receive the average intensity of the SPHERE image along the long-axis of the slit, corresponding to the spatial dimension of our 2D spectra (Sect. 3.2). We then normalized and convolved the resulting intensity profile with a high-S/N VISIR PSF (see Appendix F) at the wavelengths of the respective PAH bands for comparison. Finally, we scaled the intensity profile with r², assuming that scattered-light intensity decays as an inverse square law (except for the polarized light intensity modulation caused by projected geometry) with the stellar radiation field. This scaling is only approximate as it ignores inclination effects and is strictly accurate only for disks seen close to face-on. We note that the SPHERE observations have been performed with an apodized Lyot coronagraph with a radius of 93 mas. Therefore, any scattered light arising at smaller radii will be strongly suppressed in the SPHERE data. Hence, we may underestimate the scattered light intensities within approximately the central resolution element of our VISIR data for the 7.9 and 8.6 μm PAH bands and about the central 2/3 of a resolution element for the 11.3 and 12.6 μm bands.

In Figs. 9a-13 it is evident that overall-PAH intensities decrease less steeply with radius than the scattered light profiles. In their work on HD 163296, Muro-Arena et al. (2018) compare the strength of scattered light emission of the rings to ALMA images, and find that its outer ring is fainter than expected, and potentially attribute it to shadowing or settling. This also suggests that the scattered light intensity might show a different radial dependence than the diluted stellar radiation field. This observation, coupled with the results of Dullemond & Dominik (2004), may indicate that in all objects in our sample, dust-settling or shadowing effects may cause the discrepancy in the spatial scattered light and PAH intensities.

Another potential explanation for this discrepancy is that in regions where the different PAH/scattered-light profiles decay rates are observed, PAH molecules are saturated by incoming exciting stellar radiation (i.e., their cooling timescales are comparable to- or longer than- the photon absorption rate), in which case the PAH emission is not coupled to the diluted stellar radiation field, but instead maintains its decay power-law (Draine et al. 2021). However, if the PAH excitation is in the single photon regime (i.e., the average time between the absorption of individual exciting stellar photons is long compared to the time over which the PAH molecule cools through the emission of infrared radiation), then the intensity is expected to be directly proportional to the local stellar radiation field, and the spectral shape is the same everywhere (potentially modulated by disk structure and spatial variations in PAH abundance and molecular properties) (Draine & Li 2001; Li & Draine 2001).

To probe the relevance of this effect to the observed discrepancy, we model the heating and cooling processes of C₆₀H₂₀ and its cation () in the disk around AB Aurigae at different radial distances from the central star. The specific PAH molecule C₆₀H₂₀ is selected since it represents the typical PAH in protoplanetary disks around Herbig Ae/Be stars (Seok & Li 2017). In Fig. 9b we plot the temperature probability distribution functions and emission spectra of C₆₀H₂₀ and for two nominal projected radial distances from the central star of AB Aurigae, which are approximately the innermost and outermost disk regions of each object, probed by VISIR-NEAR. In Sect. 6.2.1–6.2.5, we discuss the comparison between the PAH and scattered-light profiles, as well as the result of this simulation for each object individually. We plot the temperature probability and emission spectra for the remaining four objects in Appendix E.

We note that while additional considerations are likely necessary to compare the two thoroughly, our purpose here is to qualitatively examine whether or not prominent morphological features (e.g., projected asymmetry, substructure) can be identified in both and address large-scale discrepancies, as both diagnostics trace similarly high vertical layers in the disk. We find the morphological comparison of the scattered light and PAH intensity profiles to be visually most evident if we scale the former by r², but not the PAH profiles, though the visual match between the profiles is not equally good in all sources.

6.2.1 AB Aurigae

The SPHERE image of this object (Fig. 2) makes it readily clear that its circumstellar disk is extremely rich in substructure, spanning across a range of different scales, and is also azimuthally asymmetric. Therefore, a 1D representation of its intensity distribution does not represent its complexity, and indeed blurs some of the subtler spatial features that fell within the slit. The rather small inclination of the disk (30° according to Boccaletti et al. 2020) negates the differences in the polarized intensity (PI). During observation, the slit was positioned roughly along the minor axis of the disk (i.e., the axis of minimal PI but maximal asymmetry).

In Fig. 9a we compare all PAH bands with a r²-scaled and convolved 1-D profile of the SPHERE intensity distribution contained within the slit (see Sect. 6.2.1). Unlike in cases where a comparison in this configuration yields the greatest similarity between the two intensity distributions (Sect. 6.2.3), here the two are not alike. The ring-like PAH distribution evident in the three bluer bands is considerably more compact than the substructure manifested in the SPHERE profile. The slight asymmetry, regarding which both the PAH and the SPHERE profiles were consistent for other objects (Sects. 6.2.2, 6.2.3, and 6.2.4), is equivocal in this case (despite the asymmetry being the smallest with respect to the aforementioned cases, with the exception of the nearly pole-on HD 169142).

As shown in Fig. 9b, PAHs will be stochastically heated if present in the disk around AB Aurigae. Therefore, the observed discrepancy between the two profiles cannot be explained by a potential saturation state of the emitting PAHs and remains an open question. Comparison between scattered light and millimeter continuum imaging of AB Aurigae is made in Tang et al. (2017) and Currie et al. (2022).

6.2.2 HD 97048

The PAH emission of this object serves as the optimal case of our analysis. It is very well spatially resolved, its spectral quality is high and there is no silicate emission present to potentially hamper the robustness of the constraints placed on the PAH emission.

Here, the slit was once again positioned almost in parallel to the short axis of the rather inclined object (≈41^° according to van der Plas et al. 2017). Such a positioning of the slit, while incidental, is optimal for our purposes, as it captures a considerable measure of the projected asymmetry of the object, thus allowing a comparison thereof between the SPHERE intensity profile and that of the PAHs.

Indeed, all bands (with the exception of the insignificantly constrained 12.6 μm band) exhibit a consistent behavior of centrally peaked emission with considerable spatial extent. Additionally, the spatial asymmetry and the consistent bump(s) located at ≈0."5 (0."75) in all three bands, are features to look for in the SPHERE intensity profile as well. We find PAHs to undergo stochastic heating in the disk around HD 97048, and therefore the radial profiles of the PAH emission bands (after being scaled by r⁻²) directly reflect the spatial distribution of PAHs in the disk (see Appendix E).

When considering Fig. 10, two important observations are made: (1) Both intensity profiles display identical asymmetry, where the right (far) disk side exhibits more prominent emission that decays more slowly. (2) The bump visible in the PAH profiles at ≈0."5 is also present in the SPHERE profile, roughly at the same location, and represents an increase in the surface density of the dust. One can argue, therefore, that here, PAH molecules' column density is coupled with that of the dust. The keen-eyed reader may notice that a second, albeit weaker, bump at ≈0."75 may also be correlated, particularly visible at the stronger, 11.3 μm band.

This object presents an outstanding potential for a detailed analysis of the properties of PAHs (e.g., size distribution and charge state) as a function of distance from the star. We hope to explore this in future work. A comparison between scattered light and millimeter continuum imaging of HD 97048 is made in van der Plas et al. (2017).

Fig. 9 PAH intensity profiles and model calculations for AB Aurigae. Panel a: Comparison of the normalized radial profiles of the three prominent PAH bands (7.9 μm, 8.6 μm and 11.3 μm, and 12.6 μm, from bottom to top, respectively) with a radial intensity profile of the reduced SPHERE image (Fig. 2). Panel b: Temperature probability distributions (left panel) and emission spectra (right panel) of neutral C60H20 (solid lines; labeled "PAH") and its cation (dashed lines; labeled "PAH+") for two nominal distances: 0."228, which corresponds to ≈40AU from the central star (red lines), and 1."14, which corresponds to «200AU from the central star (blue lines), respectively, in the innermost and outermost disk regions probed by VISIR-NEAR. The temperature probability distribution functions for C60H20 and are broad even for the innermost disk region, indicating that PAHs undergo stochastic heating by individual stellar photons even in the innermost disk region, where the photon absorption rate is ~25 times more frequent compared to that in the outermost disk region. The stochastic-heating nature of PAHs in the disk around AB Aurigae is also indicated by the close resemblance of the model emission spectra calculated for C60H20 and at 40AU (solid and dashed red lines) to that at 200AU, when scaled to the same radial distance, (solid and dashed blue lines).

Fig. 10 Comparison of the normalized radial profiles of the three prominent PAH bands of HD 97048 (similar to Fig. 9a).

6.2.3 HD 100546

This object displays the best agreement between the PAH intensity profiles and that of the r² -scaled SPHERE image (with the exception of the 11.3 μm band for reasons explained in Sect. 5.2.3).

Here too, similarly to AB Aurigae and HD 97048 (Sect. 6.2.1, 6.2.2), the slit is positioned roughly in parallel to the minor axis of the inclined circumstellar disk (≈41° according to Ginski et al. 2016), and the minimal PI is compensated by the fact that the minor axis is in fact that of most interest to this comparison, as it potentially exhibits the greatest asymmetry of the projected disk.

When considering Fig. 11, it is evident that the similarity between the PAHs and the SPHERE intensity is considerably greater than for AB Aurigae and HD 97048. Indeed, the 7.9 μm nearly completely overlaps the SPHERE profile and reproduces the two smeared bumps extending along the left lobe of the profile until ≈−1"(albeit not perfectly). Surprisingly, it is also considerably more spatially extended than the 8.6 μm and 12.6 μm profiles, which one might expect to follow a similar spatial pattern (as seen in the previous two cases). The 8.6 μm and 12.6 μm bands, while more centrally centered and compact, still roughly match the inner region of the SPHERE intensity profile, although they decay spatially much faster.

This behavior may be reminiscent, for example, of the measure of coupling of a dominant charge state to the dust. Exploring this possibility further is, unfortunately, out of the scope of this study.

This object emphasizes best the peculiar observation that maximum similarity between the two intensity profiles is reached when the SPHERE profile is scaled by r² while the PAHs are not (esp. for the 7.9 μm band). A comparison between scattered light and millimeter continuum imaging of HD 100546 is made in Pineda et al. (2019).

6.2.4 HD 163296

The PAH profiles of this object are of relatively low S/N and allow only a limited qualitative comparison with the SPHERE profile. In this comparison, we only consider the 7.9 μm and 8.6 μm PAH bands, as the 12.6 μm is of too low S/N and is not constrained, and the 11.3 μm band may experience the same bias as in HD 100546 (see Sect. 5.2.4).

To our advantage serves, again, the fact that the object is inclined by ≈46^° with a position angle of ≈−88.8^° (Diep et al. 2019). As such, the slit was positioned nearly parallel to the minor axis of the circumstellar disk during observation, thus allowing us to probe the axis of maximal asymmetry.

When considering Fig. 12 it is evident that while the PAH profiles are devoid of much detail, the asymmetry is detectable in both profiles, and is roughly in agreement concerning the position of the dominant ring-like structure. However, it is also evident that the SPHERE image is much richer in substructure extending farther from the star, whereas the PAH emission decays abruptly, and we cannot detect spatial features beyond ≈0."5 from the star. Unfortunately, the S/N of the PAH profiles does not allow for comparison with the scattered light profiles or investigation of spatial variations in the PAH properties; the two bluest bands, which we detect significantly, represent the same dominant charge state. As for the other stars, multiple-photon excitation is not important on the scales resolved by the VISIR-NEAR observations.

A comparison between scattered light and millimeter continuum imaging of HD 163296 is made in Muro-Arena et al. (2018).

Fig. 11 Comparison of the normalized radial profiles of the three prominent PAH bands of HD 100546 (similar to Fig. 9a).

Fig. 12 Comparison of the normalized radial profiles of the three prominent PAH bands of HD 163296 (similar to Fig. 9a).

6.2.5 HD 169142

The PAH emission in this object appears centrally peaked close to the star, though the spatial profile is notably wider than the VISIR-NEAR PSF and has a distinct shape. In Sect. 5.2.5 we showed that the PAH profile is consistent with a ring-like emission, analogous to that seen in scattered light, but with a smaller ring radius (see Figs. 2 and 8).

Here, the circumstellar disk is nearly pole-on (with an inclination of ≈ 13° (Macías et al. 2019), so we expect essentially no asymmetry to be detected. Indeed, in Fig. 13, we see that nearly all PAH emission (> 1σ) is contained within the innermost cavity of the circumstellar disk. We find PAHs to be stochastically heated in the inner disk and emit prominently in the 6.2, 7.7, 8.6, and 11.3 μm bands (see Appendix E).

A comparison between scattered light and millimeter continuum imaging of HD 169142 is made in Macías et al. (2019); Pérez et al. (2020).

7 Summary

We conducted an AO-supported N-band long-slit spectroscopy survey of eight protoplanetary disks and spatially resolved the PAH emission in five objects: AB Aurigae, HD 97048, HD 100546, HD 163296, and HD 169142. Using a nonlinear optimization routine, we fit the spatially resolved spectra, which we parameterized with a continuum component, five silicate species of two nominal grain sizes, and four PAH bands (7.9, 8.6, 11.3, and 12.6 μm). For AB Aurigae, the three bluer bands exhibit a ring-like emission with a radius of roughly 0."3. For HD 97048, we find an asymmetric extended, centrally peaked emission with varying decay rates for each of the three bluer PAH bands. For HD 100546, we find a central cavity in the emission profiles of all PAH bands except for the 11.3 μm band, where a degeneracy exists between the emission in the PAH feature and that of the abundant crystalline forsterite, leading to a false-positive centrally peaked intensity profile. For HD 163296, we similarly find that the extracted 11.3 μm PAH profile is centrally peaked and possibly confused with forsterite emission, whereas the two bluest bands exhibit a roughly identical ringlike structure. For HD 169142, we observe a widened centrally peaked uniform PAH emission across all four bands, which is consistent with a ring-like geometry with a radius that is somewhat smaller (≈0."41) than that of the ring observed in scattered light (0."17).

Within our small sample, we find that in sources with strong silicate features, the PAH emission appears to be ringlike (AB Aurigae, HD 100546, and HD 163296), whereas, in sources without silicate emission, the emission tends to be centrally peaked (most notably in HD 97048), which can be due to the competition of the silicate grains and PAHs over UV photons (and, subsequently, a considerably stronger emission feature of the first), dust-settling and shadowing effects, and/or the coagulation of PAHs on the small dust grains that are abundant in the inner regions of disks with strong silicate emission.

Finally, we compared the extracted spatial profiles of the PAH emission of each object with a 1D intensity profile of a scattered light image obtained with SPHERE, as the two probe the highest vertical layers of the disk and are therefore worthwhile to compare. The main conclusions of this analysis are as follows:

• –

  PAH emission arises from the highest disk surface layer (similar to ¹²CO). Scattered light observations trace a somewhat deeper layer, but one that is still spatially close to the PAH emitting layer.

• –

  Naively, we expect both the scattered light and PAH intensity profiles to dilute approximately as r⁻² with the stellar radiation field. This is modulated by inclination, disk substructure, and spatially variable dust-settling and shadowing effects.

• –

  We find that the PAH intensities decrease radially less steeply than the scattered-light profiles. We considered two potential reasons for this: (1) We investigated through modeling whether it could be due to the PAH excitation being saturated in the inner disk regions. However, we find that, on the scales probed by our observations, PAHs undergo stochastic heating by individual stellar photons even in the innermost disk regions in all objects in our sample. (2) Dustsettling and shadowing effects impact dust grains more so than PAHs, thus causing a steeper radial intensity decay of the former. Objects with a ring-like PAH intensity profile exhibit richness in silicate emission, whereas the centrally peaked HD 97048 lacks significant silicate emission.

• –

  Despite the difference in intensity decay, steepness, scattered-light, and PAH intensity profiles agree on prominent projected and intrinsic geometric features, such as inclination effects, significant gaps, and rings.

Fig. 13 Comparison of the normalized radial profiles of the three prominent PAH bands of HD 169142 (similar to Fig. 9a).

Acknowledgements.

We thank Dr. Christian Ginski for the SPHERE images and the anonymous referee for offering valuable input concerning the editing process of this work.

Appendix A Target sample: Detailed introduction

Appendix A.1 AB Aurigae

The disk around this object, extending up to ≈550 au (Hashimoto et al. 2011) is extremely rich in substructures: gaps at distances of ≈30, 70 au were detected in, for example, millimeter CO isotopologue observations with the Institut de radioastronomie millimetrique (IRAM) interferometer (Piétu et al. 2005), near-infrared polarimetry with the Subaru telescope (Hashimoto et al. 2011) and dust continuum and ¹²CO observations with ALMA (Tang et al. 2017). These observations also detected spiral structures on different scales. In addition, exquisite images of the circumstellar disk as observed with SPHERE showcase complex spiral structures, as recently reported by Boccaletti et al. (2020), and further constrain a "bump" in an inner spiral, possibly induced by a ≈4-13 M_Jup companion at a distance of ≈40 au (Poblete et al. 2020).

The object's spatially unresolved N-band spectrum was modeled by, for example, van Boekel et al. (2005). While abundant silicate emission was observed and modeled, PAH emission was not detected. Later, Juhász et al. (2010) detected several PAH bands in the Spitzer spectrum of the object.

Appendix A.2 HD95881

Verhoeff et al. (2010) performed robust fitting of a radiative-transfer model to reproduce multiple observables of the pro-toplanetary disk around this object across different bands-spanning from Q-band imaging, infrared spectroscopy, and K-and N-band interferometry. The derived density and temperature structure of the disk indicates that it consists of a thick puffed-up inner rim and an outer region that has a flaring gas surface and is relatively void of "visible" dust grains, and conclude that the disk is in a transition phase from a gas-rich flaring disk to a gas-poor self-shadowed disk.

More recently, Järvinen et al. (2018) detected weak mean longitudinal stellar magnetic fields using the High Accuracy Radial velocity Planet Searcher polarimeter (HARPSpol) attached to the ESO 3.6 m telescope at the La Silla observatory.

Appendix A.3 HD 97048

The well-known circumstellar disk around HD 97048 has been thoroughly investigated. This young object (~5 Myr according to Fairlamb et al. (2015)) exhibits a particularly extended disk (up to ~1000 AU), which was observed at an early stage with the Hubble Space Telescope (Doering et al. 2007) and since then on multiple additional occasions. Its extended emission (e.g., Lagage et al. 2006; Ginski et al. 2016), substructure (e.g., Ginski et al. 2016; van der Plas et al. 2017), and chemistry (e.g., Siebenmorgen et al. 2000; Meeus et al. 2012; Booth et al. 2019) have been probed.

Observations of, for example, millimeter-sized continuum emission observed with ALMA (van der Plas et al. 2017; Walsh et al. 2016), tracing the continuum emission up to ≈355 au, and polarized scattered light observations with SPHERE (Ginski et al. 2016) reveal gaps on distance scales of ≈40-350 au, as well as tracing the continuum emission to ≈355 au.

A study by Quanz et al. (2012) focusing on the inner region of the disk (~15-160 au) with the VLT/NACO polarimetric differential imager was able to detect and resolve dust continuum emission in the H and K bands. Studies aimed at specifically probing the properties of PAHs (van Boekel et al. 2004; Lagage et al. 2006) were able to resolve the distribution thereof spatially. Finally, Pinte et al. (2019) report a planet discovered by its kinematic signature in the ALMA line observations data.

Appendix A.4 HD 100453

The disk around this object possesses ample and diverse substructures, such as spiral arms and shadows (e.g., Benisty et al. 2016), a gap approximately of the size of ~1-20 au (Wagner et al. 2015; Khalafinejad et al. 2016), and finally-an M dwarf companion at a distance of ≈120 au (e.g., Wagner et al. 2015; Dong et al. 2015), with which binary-disk interaction (namely the spiral arms of the outer disk) is evident (Rosotti et al. 2020). Later, Rosotti et al. (2020) presented clear evidence for that hypothesis.

While no direct detection of planets embedded in the disk has been made, Nealon et al. (2020) posit a potential existence thereof as an explanation for the mutual inclination between the inner and outer disk components (Bohn et al. 2022). Modeling the near-infrared intensity distribution of the disk based on inter-ferometric observations, Klarmann et al. (2017) found evidence for stochastically heated particles.

Appendix A.5 HD 100546

Previous observations and modeling of the disk around this object traditionally divide it into an inner and outer component (Bouwman et al. 2003; Panić et al. 2014), separated by a gap between approximately 1 and 10 au (e.g., Menu et al. 2015; Jamialahmadi et al. 2018)), with a possible misalignment of the two (Kluska et al. 2020). Additional more delicate substructure features, such as spirals (Boccaletti et al. 2013; Follette et al. 2017) and inflow or jets (Mendigutía et al. 2017) had been postulated and later found by Schneider et al. (2020) to be most likely a jet. The same is true for posited existence of companions at different distances throughout the disk (e.g., Quanz et al. 2013; Brittain et al. 2014; Cugno et al. 2019; Casassus & Pérez 2019; Pérez et al. 2020). van Boekel et al. (2004) were able to spatially resolve the N-band PAH emission of this object with the TIMMI2 mid-infrared instrument on the ESO 3.6 m telescope at La Silla Observatory, albeit with a seeing-limited spatial resolution of only ≈0."85.

Appendix A.6 HD 163296

Substructure in the disk of this object had been identified across a wide spectral range of observations, predominantly in the submillimeter with ALMA (e.g., Isella et al. 2018) and in near-infrared polarimetry (e.g., Garufi et al. 2014). These constrained the presence of gaps on distances around ≈ 140, 80, 50, and 10 au. Subsequently, companions were suggested as a potential mechanism to induce these gaps (e.g., Liu et al. 2018; Teague et al. 2018; Pinte et al. 2018). More recently, evidence for a temporally variable sub-structure on sub-au scales was found in VLT Interferometer (VLTI) observations using Multi AperTure mid-Infrared Spectroscopic Experiment (MATISSE) (Varga et al. 2021) and GRAVITY (Sanchez-Bermudez et al. 2021). Finally, Petit dit de la Roche et al. (2021) recently observed this object with VISIR-NEAR in imaging mode and spatially resolved images in two different PAH narrowband filters.

Appendix A.7 HD 169142

This object has been studied and observed with various methods across a wide swath of the electromagnetic spectrum-from the submillimeter and millimeter observations with ALMA (Macías et al. 2019), near-infrared (Fukagawa et al. 2010) and mid-infrared (Okamoto et al. 2017). Additionally, scattered light polarimetry surveys with SPHERE (Bertrang & Avenhaus 2018) produced images of the multiple rings (up to four Macías et al. (2017); Pérez et al. (2019)), and interferometric surveys in the near-infrared were also conducted with the VLTI-The Precision Integrated-Optics Near-infrared Imaging Experiment (PIONIER) instrument in order to constrain the inner-region substructure of the disk (Chen et al. 2018, 2019). The innermost ring (inner disk rim, Bertrang & Avenhaus (2018)) and the outer three rings are located at distances of ≈20 and ≈29, 55, 74, and 87 au, respectively. While planets were not directly detected in the disk, their potential presence at different distance scales or that of ongoing planet formation was postulated by, for example, Osorio et al. (2014), Reggiani et al. (2014), and Biller et al. (2014), but are inconclusive. More recently, a scattered light angular differential imaging study by Gratton et al. (2019) with SPHERE suggested the presence of a ≈2 M_Jup planet at a distance of ~ 100-150 au. Strong spatially unresolved PAH emission features at 3.3, 6.2, 7.7, 8.6, 11.3, and 12.6 μm were detected by Seok & Li (2016) and were recently spatially resolved by Petit dit de la Roche et al. (2021) with VISIR-NEAR in the imaging mode.

Appendix A.8 HD 179218

While a less studied object, its circumstellar disk manifests mainly through its infrared excess (Meeus et al. 2001) and is subsequently classified as a group-Ia source (i.e., of possible flared geometry), of which PAHs are an optimal tracer. Indeed, Meeus et al. (2001) and Bouwman et al. (2001) both detect PAH emission at the ≈8.6 and 11.3 μm regions, which was later confirmed by Juhász et al. (2010). Additionally, Meeus et al. (2001) and Bouwman et al. (2001) detected silicate emission, which was later confirmed and modeled by van Boekel et al. (2005). More recently, Taha et al. (2018) observed the object with the Canari-Cam (Telesco et al. 2003) high-resolution mid-infrared imager and spectrograph and were able to resolve the N-band PAH emission spatially, concluding that similarly to HD 97048–PAHs extend far across the flared disk's surface layer in a predominantly ionized state. Recent MATISSE-VLTI preliminary low spectral resolution observations of this object could not infer the presence or not of a PAH emission in the inner disk region (< 10 au). Stochastically heated carbon dust particles are, however, probably populating this inner disk region (Kokoulina et al. 2021).

Appendix B Parametric model: Template plots

Fig. B.1 Templates used in the parametric model. Upper Panels: Silicate species templates (Si) of five silicate species (Table B.1), where for each specie there exit two templates–corresponding to two prescribed grain-sizes of 0.1 μm ("small") and 1.5 μm ("big"). Bottom Panels: PAH component templates (Pj).

Appendix C Statistical model

We evaluated the goodness of fit to the reduced N-band spectra for each simulated model with parameter set Θ assuming a normal distribution, whose log is (C.1)

where there are m observed spectral measurements (1024 pixels) x_i with respective uncertainties σ_i, and y_i is the simulated N-band emission spectrum associated with the model parameters Θ. As the measurement uncertainty of the observed flux values, σ_i, we adopt the standard deviation of the flux over all objects' exposures.

Our derived best-fit values of the model parameters are the median values of the posterior distribution, and the uncertainties are the 50 ± 34.1 percentiles ranges of the posterior distribution, computed after removal of the MCMC "burn-in" phase (points with relative likelihood < 10⁻³ times that of the best-fit). We refer to the range of uncertainties as σ, defined irrespective of the distribution. For a normal distribution, σ corresponds to the standard deviation. We note that the uncertainty derived for each parameter in the model from its respective posterior distribution, which we refer to as σ, is not the same as the σ_i in Eq. (C.1), which represents the uncertainty in the observed fluxes.

Appendix D Unresolved spectra posterior distributions

In Figs. D.1–D.7, we plot the posterior distribution of all model constituents (see Sect. 4.2.2) of the unresolved spectra fitting procedure (see Sect. 5.1) for all objects in our sample.

Fig. D.1 AB Aurigae: unresolved fitting posterior distributions.

Fig. D.2 HD 95881: unresolved fitting posterior distributions.

Fig. D.3 HD 97048: unresolved fitting posterior distributions.

Fig. D.4 HD 100453: unresolved fitting posterior distributions.

Fig. D.5 HD 100546: unresolved fitting posterior distributions.

Fig. D.6 HD 163296: unresolved fitting posterior distributions.

Fig. D.7 HD 169142: unresolved fitting posterior distributions.

Fig. D.8 HD 179218: unresolved fitting posterior distributions.

Appendix E Temperature probability distribution and PAH model emission spectra

In Figs. E.1–E.4, we plot the temperature probability distribution and PAH model emission spectra (see Sect. 6.2) for HD 97048, HD 100546, HD 163296, and HD 169142, respectively.

Fig. E.1 HD 97048: temperature probability distributions (left panel) and emission spectra (right panel) of neutral C60H20 (similar to Fig. 9b).

Fig. E.2 HD 100546: temperature probability distributions (left panel) and emission spectra (right panel) of neutral C60H20 (similar to Fig. 9b).

Fig. E.3 HD 163296: temperature probability distributions (left panel) and emission spectra (right panel) of neutral C60H20 (similar to Fig. 9b).

Fig. E.4 HD 169142: temperature probability distributions (left panel) and emission spectra (right panel) of neutral C60H20 (similar to Fig. 9b).

Appendix F PSF model

The inner disk regions dominate the total flux of our targets but are either spatially unresolved or marginally resolved with VISIR-NEAR. The PSF redistributes the flux arising from these regions in the VISIR focal plane, leading to a small fraction of the light landing in the off-axis apertures. In order to estimate this fraction we create a high-S/N PSF model as follows.

For each of our calibrator observations, we created a 2D spectrum (dimensions: across-slit (dispersion) and along-slit (spatial) direction), which we normalized such that the flux in each spectral bin is unity. The spectral resolution of VISIR in the mode employed here is approximately R = 300; hence, we measured a quasi-monochromatic PSF at each wavelength.

We over-sampled the 2D spectra of the calibrators by a factor of 2 in the spatial dimension, aligned them, and combined them into one high-S/N 2D spectrum. We then resample the spatial profiles in each spectral bin to a common grid in terms of λ/D (where D is the VLT primary mirror diameter of 8.2 m), over-sampling by a factor of 2.

We spectrally averaged these quasi-monochromatic PSFs to produce a single, 1D "super PSF" along the spatial dimension, which is displayed in Fig. F.1. More than ten diffraction rings are visible.

There is an optical ghost approximately 4 arcsec below the target on the VISIR detector (toward the left of the profile in Fig. F.1, as the gray part of the curve left of ≈ −10λ/D). We replaced this part of our PSF profile with the continuation of a linear fit between 9λ/D and 4λ/D. We note that the ghost is far enough away from the target to not interfere with any potential spatially extended emission that is the topic of this study. Finally, we resampled the super PSF at each wavelength back onto the instrumental spatial grid of the VISIR detector.

Fig. F.1 VISIR-NEAR quasi-monochromatic PSF as recorded on our calibrators.

Appendix G Potential PSF-widening due to atmospheric conditions

To affirm the assumption that our detected spatial distributions of the PAHs represent the intrinsic structure of the observed objects, rather than a misinterpretation of a PSF widened due to sub-optimal atmospheric conditions, we conducted a series of tests to detect potential nonlinear widening of the science targets' PSF (i.e., not scaled as ≈ λ/D) at particular spectral bins where such widening might be expected, for example, the 9-10 μm ozone band and the 7.7-7.9 μm water vapor band. In the former, no PAH features are present, so a widening of the PSF in that band is not necessarily indicative of the intrinsic distribution of PAHs within the object, whereas in the latter the opposite is true. Thus, if widening as a function of the PSF is detected in both bands, one would have reason to suspect the derived PAH spatial distributions. However, the absence of such widening in the ozone band would serve to show reinforce our conviction in the data. Additionally, we perform these tests on the calibrators as well, which should ideally exhibit no widening of the PSF at all, if it indeed reflects the intrinsic structure of the observed objects. We plot our results for the calibrators in Fig. G.1 and for the science targets in Fig. G.2.

Fig. G.1 indicates that with the exception of a single faulty observation of HD 92305 (which was not used as a calibrator in the processing stage), the PSF residual widths of calibrators do not seem to exhibit a strong correlation with the airmass nor the PWV in either band. The same can be said about the residual widths of the science targets' PSF in the ozone band (HD 169142 is spatially resolved in the continuum and therefore deviates). By observing the water-vapor bands of the ozone band, one might assume, however, a strong correlation and be reluctant to attribute those widths to intrinsic energy distributions within the objects themselves. However, the first three cases serve to show that there seems to be no correlation between the given atmospheric conditions under scrutiny and the PSF residual widths, which serves to reinforce the assumption that the observed intensity distributions in the water-vapor (or 7.9 μm PAH) band of the science targets is not an atmospheric effect, but rather an intrinsic property of their structure.

Fig. G.1 Calibrator PSF residual width (in pixels) vs. atmospheric conditions. The residual widths are calculated as follows: (1) each column (i.e., spectral coordinate; see Fig. 3) of the detector is fitted with a Gaussian to estimate the width of the PSF. (2) The fitted widths are then fitted with a first-order polynomial (ideally, λ/D). (3) The fitted polynomial is then subtracted from the fitted widths to generate a residual widths vector -which represents the nonlinear deviations of the fitted widths. Left Panel: PSF residual width vs. airmass during observation for the water vapor and ozone bands. Right Panel: PSF residual width vs. precipitable water vapor (PWV) during observation for the water vapor and ozone bands.

Fig. G.2 Science target PSF residual width (in pixels) vs. atmospheric conditions (similar to Fig. G.1).

References

All Tables

All Figures

Fig. 1 Variations in the effective emission height of different tracers for a 3D disk demonstrated for HD 97048: (1) scattered-light as observed with SPHERE (Ginski et al. 2016), (2) PAH emission as observed with VISIR (Lagage et al. 2006), and (3) 13/12CO as observed with ALMA (Walsh et al. 2016). The inset on the bottom right is Fig. 4 from Lagage et al. (2006), demonstrating the 2D spatial distribution of the 8.6 μm PAH emission intensity of HD 97048.

In the text

	Fig. 2 Mean slit orientation (yellow rectangle) overlaid on scattered-light images observed with SPHERE for all targets in our sample except HD 95881.
In the text

	Fig. 3 Calibrated 2D spectrum of HD 97048. The x-axis is the spectral dimension perpendicular to the long axis of the slit, and the y-axis is the spatial axis parallel thereto. The source is positioned in the center of the y-axis. Spatially resolved emission in the PAH bands is highlighted with white arrows. Note the image was square-root-stretched for demonstration purposes.
In the text

Fig. 4 Excess flux analysis of the 8.6 μm PAH band of HD 97048. Top panel: Gaussian fitting of the 8.6 μm PAH template, with a resulting standard deviation of σtemplate ≈ 0.14 μm. The aperture of the PAH spectral window is defined as ±1.17 σtemplate, and the continuum aperture is set to be two FWHMs down the red part of the spectrum (the chosen widths for each are plotted atop the panel). Bottom panel: Fluxes integrated over the continuum and PAH apertures, and the computed excess flux (blue, orange, and green curves, respectively) and their respective ±1σ uncertainties (see Appendix C), indicated by curve widths.

In the text

Fig. 5 Normalized flux excess in the 8.6 μm and 11.3 μm PAH bands (top and bottom panels, respectively). The width of the lines represents the ±1σ residuals of the PAH and continuum windows (see Sect. 4.1). The x-axis is parallel to the slit (i.e., its spatial dimension), and the zero point is the object's center (the star). Note that negative values indicate that the flux in the continuum windows is greater than that in the PAH window. This is caused by silicate emission in the continuum window.

In the text

Fig. 6 Calibrated long-slit spectra of our targets. Left panel: Spatially resolved spectra. The black curves are down-slit spatially resolved spectra. Each spectrum represents an average of a ≈ λ/D-wide spectral bin (every three rows of the reduced 2D spectrum, whose centroids are separated by 0″.228 and are arranged bottom up, where the bottom spectrum is that of the central object). The uppermost spectrum, farthest from the star, is averaged over three such apertures. The blue curves are the expected flux contamination coming from the unresolved central source, extracted as described in Sect. 5.1. The off-axis direction is indicated by the two letters in the top or bottom panels (e.g., NE for northeast). Middle panel: Fitted integrated spectra. The solid black line is the observed spectrum, the thicker orange line is the best-fit spectrum, and the filled gray area is the ±1σ error-range, where σ is the flux-wise standard deviation of spectra generated with all best-fit parameters within the 50 ± 34 percentiles range. Right panel: Similar to the left panel, but on the opposite side of the star with respect to the slit axis. Note that spectra marked with an asterisk are identical in both the left and right panels. Note also that the spectral artifact visible at ≈9.6 μu is caused by the atmospheric ozone absorption band. Continued on the next page.

In the text

Fig. 7 Radial profile of the normalized PAH emission in four spectral bands for select sources. Sources are ordered column-wise, and the PAH bands are ordered row-wise and are color-coded (from 7.9 μm in the uppermost row to 12.6 μm in the bottom row). The filled area represents the ±1σ error range around the median value of the emission amplitude posterior distribution of the relevant PAH template (see Appendix C). The orientation of the spatial dimension (parallel to the long axis of the slit) is denoted in the upper-right corner of the uppermost panel. Additionally, the dashed gray line is the high-S/N synthesized PSF (Appendix F), plotted for reference.

In the text

Fig. 8 Spatial emission profile of HD 169142 in the 11.3 μm PAH band. The red curve shows the observed PSF profile (see Appendix F), and the black curve shows the employed approximation thereof used in the 2D convolution of the model (see Sect. 5.2.5). The solid gray curve shows the ring model that best fits the observed central peak (ring radius 0."14). The dashed gray curve shows the best-fit model with the ring radius fixed to the value observed in the SPHERE data (0."17).

In the text

Fig. 9 PAH intensity profiles and model calculations for AB Aurigae. Panel a: Comparison of the normalized radial profiles of the three prominent PAH bands (7.9 μm, 8.6 μm and 11.3 μm, and 12.6 μm, from bottom to top, respectively) with a radial intensity profile of the reduced SPHERE image (Fig. 2). Panel b: Temperature probability distributions (left panel) and emission spectra (right panel) of neutral C60H20 (solid lines; labeled "PAH") and its cation (dashed lines; labeled "PAH+") for two nominal distances: 0."228, which corresponds to ≈40AU from the central star (red lines), and 1."14, which corresponds to «200AU from the central star (blue lines), respectively, in the innermost and outermost disk regions probed by VISIR-NEAR. The temperature probability distribution functions for C60H20 and are broad even for the innermost disk region, indicating that PAHs undergo stochastic heating by individual stellar photons even in the innermost disk region, where the photon absorption rate is ~25 times more frequent compared to that in the outermost disk region. The stochastic-heating nature of PAHs in the disk around AB Aurigae is also indicated by the close resemblance of the model emission spectra calculated for C60H20 and at 40AU (solid and dashed red lines) to that at 200AU, when scaled to the same radial distance, (solid and dashed blue lines).

In the text

	Fig. 10 Comparison of the normalized radial profiles of the three prominent PAH bands of HD 97048 (similar to Fig. 9a).
In the text

	Fig. 11 Comparison of the normalized radial profiles of the three prominent PAH bands of HD 100546 (similar to Fig. 9a).
In the text

	Fig. 12 Comparison of the normalized radial profiles of the three prominent PAH bands of HD 163296 (similar to Fig. 9a).
In the text

	Fig. 13 Comparison of the normalized radial profiles of the three prominent PAH bands of HD 169142 (similar to Fig. 9a).
In the text

	Fig. B.1 Templates used in the parametric model. Upper Panels: Silicate species templates (Si) of five silicate species (Table B.1), where for each specie there exit two templates–corresponding to two prescribed grain-sizes of 0.1 μm ("small") and 1.5 μm ("big"). Bottom Panels: PAH component templates (Pj).
In the text

	Fig. D.1 AB Aurigae: unresolved fitting posterior distributions.
In the text

	Fig. D.2 HD 95881: unresolved fitting posterior distributions.
In the text

	Fig. D.3 HD 97048: unresolved fitting posterior distributions.
In the text

	Fig. D.4 HD 100453: unresolved fitting posterior distributions.
In the text

	Fig. D.5 HD 100546: unresolved fitting posterior distributions.
In the text

	Fig. D.6 HD 163296: unresolved fitting posterior distributions.
In the text

	Fig. D.7 HD 169142: unresolved fitting posterior distributions.
In the text

	Fig. D.8 HD 179218: unresolved fitting posterior distributions.
In the text

	Fig. E.1 HD 97048: temperature probability distributions (left panel) and emission spectra (right panel) of neutral C60H20 (similar to Fig. 9b).
In the text

	Fig. E.2 HD 100546: temperature probability distributions (left panel) and emission spectra (right panel) of neutral C60H20 (similar to Fig. 9b).
In the text

	Fig. E.3 HD 163296: temperature probability distributions (left panel) and emission spectra (right panel) of neutral C60H20 (similar to Fig. 9b).
In the text

	Fig. E.4 HD 169142: temperature probability distributions (left panel) and emission spectra (right panel) of neutral C60H20 (similar to Fig. 9b).
In the text

	Fig. F.1 VISIR-NEAR quasi-monochromatic PSF as recorded on our calibrators.
In the text

Fig. G.1 Calibrator PSF residual width (in pixels) vs. atmospheric conditions. The residual widths are calculated as follows: (1) each column (i.e., spectral coordinate; see Fig. 3) of the detector is fitted with a Gaussian to estimate the width of the PSF. (2) The fitted widths are then fitted with a first-order polynomial (ideally, λ/D). (3) The fitted polynomial is then subtracted from the fitted widths to generate a residual widths vector -which represents the nonlinear deviations of the fitted widths. Left Panel: PSF residual width vs. airmass during observation for the water vapor and ozone bands. Right Panel: PSF residual width vs. precipitable water vapor (PWV) during observation for the water vapor and ozone bands.

In the text

	Fig. G.2 Science target PSF residual width (in pixels) vs. atmospheric conditions (similar to Fig. G.1).
In the text