© The Authors 2026

1 Introduction

Photodissociation regions (PDRs) are the transition between the ionized and molecular regions in the interstellar medium (ISM), appearing at the edges of interstellar clouds when they are irradiated by a far-ultraviolet (FUV) source. Typically, there is a sharp HII-to-HI transition or ionization front (IF), and an HI-to-H₂ transition or dissociation front (DF), where the location and length scale are driven by the attenuation of the FUV radiation field and the heating, cooling, and photochemical processes it influences (Hollenbach & Tielens 1997, 1999; Wolfire et al. 2022). In nearby PDRs (distance ≲1000 pc) with an edge-on orientation, these transitions can be spatially resolved, which is why they are considered ideal laboratories for studying interstellar gas and dust and its response to changes in the environment (i.e., radiation field, density, temperature, and molecular fraction).

Typically, PDRs are very bright in the near-infrared (NIR) and mid-infrared (MIR), as the UV excitation of small carbonaceous species in the ISM drives radiative cascades through their vibrational states, resulting in a group of emission bands spanning from around 3 to 20 µm (Leger & Puget 1984; Allamandola et al. 1985). Most of these bands are associated with vibrational modes of bonds found in polycyclic aromatic hydrocarbons (PAHs) and they are also referred to as the aromatic infrared bands (AIBs). The brightest emission features of this nature are found at 3.3, 5.2, 6.2, 7.7, 8.6, 11.0, 11.3, 12.7, and 16.4 µm, and each of these bright bands appears to exhibit shape variations, or consists of multiple components (Peeters et al. 2002, 2004; van Diedenhoven et al. 2004; Chown et al. 2024; Van De Putte et al. 2024). JWST revealed one of the most detailed views of the AIB spectrum with its observations of the Orion Bar PDR, revealing numerous weaker subcomponents that blend with the main bands and affect the shape of their emission profiles (Chown et al. 2024). Diagnostics based on the intensities of these bands are used in observational studies to detect the evolution of their physical properties, such as the size or charge of the carriers, and these changes are commonly interpreted as tracers for the environment, such as the radiation field or density (e.g., Galliano et al. 2008; Pilleri et al. 2012; Stock & Peeters 2017). Highly spatially resolved studies of edge-on PDRs minimize the blending between such different conditions within the beam size or along the line of sight.

NGC 7023 is an isolated reflection nebula for which HD 200755, the B2Ve binary star at its center, is the single driver of the illumination and the PDR physics. It is of particular interest for the study of fullerenes (C₆₀) in space, observed in the central cavity (Sellgren et al. 2007, 2010; Berné & Tielens 2012). Several PDRs and their DFs are observed as bright illuminated ridges bordering this cavity. Integral field spectroscopy mosaics were observed with NIRSpec and MIRI as part of the JWST Guaranteed Time Observations (GTO) program #1192 (Misselt et al. 2025). In the context of the available JWST observations of PDRs, NGC 7023 is driven by a soft radiation field of moderate intensity (T_eff ∼ 18 000 K, G₀ ∼ 2.6 × 10³; Chokshi et al. 1988). Comparatively, the radiation field driving the Horsehead is hard with a low intensity (T_eff ∼ 33 000 K, G₀ ∼ 100; Habart et al. 2005), while that driving the Orion Bar is hard with a high intensity (T_eff ∼ 39 000 K, G₀ ∼ 2 × 10⁴; Marconi et al. 1998; Peeters et al. 2024), Here, G₀ is the radiation field in units, as defined by Habing (1968), and G₀ = 1.7 represents the average radiation field in the local ISM (Draine 1978).

The above radiation field properties result in an extended atomic hydrogen region of the PDR, while an IF and the typical H recombination lines are not observed in our field of view. The emission bands are much brighter compared to cases such as the Horsehead (Misselt et al. 2025), while the thermal dust continuum longward of 15 µm is much weaker than that observed in the Orion Bar (Van De Putte et al. 2024). The above properties allow for a highly detailed decomposition of the AIBs over the entire 3–20 µm wavelength range for individual spatial elements of the integral field units (IFUs), and the limited dust continuum emission makes NGC 7023 especially suitable to study the 16–18 µm complex. The general importance of the 16–18 µm complex lies in their connection to C-C-C bending modes, as discussed in previous Spitzer-based works (Peeters et al. 2004; Boersma et al. 2010). Some studies in NGC 7023 revealed that the 16.4 and 17.4 µm spatial distributions may be different, possibly related to the presence of different ionization stages of PAHs (Shannon et al. 2015, 2016). The spectro-spatial data and emission feature maps from the Spitzer era (Werner et al. 2004; Sellgren et al. 2010) also showed that C₆₀ is present closer to the star, based on the spatial distribution of its main emission feature at 18.9 µm, which peaks in the central cavity. Another characteristic C₆₀ feature at 17.4 µm is blended with the AIB at the same wavelength and, as such, it exhibits two local maxima: one near the peak of the 18.9 µm emission and one near the bright northwest (NW) filaments (the DFs of the PDR). The region where C₆₀ is observed does not fall within the field of view of the JWST spectroscopy presented here, although intermediate species may be observed between the C₆₀ region and the DF (Berné & Tielens 2012). However, at the spatial resolution of Spitzer, the separation between the emission (sub)components is not clearly established.

We aim to leverage the JWST spectroscopy to investigate the spatial distributions of the carbonaceous emission features in detail. We placed a particular focus on the 16.4 and 17.4 µm bands and made use of the entire 3–20 µm wavelength range to decompose the other prominent feature complexes. We show which bands or components have spatial distributions that match these 16–18 µm emission subcomponents. In Sect. 2, we summarize our processing of the JWST NIRSpec and MIRI IFU data. The methods of Sect. 3 decompose the spectra of each spaxel using the PAHFIT tool and the results are presented as maps of each subcomponent in Sect. 4. In Sect. 5 we discuss which spatial distributions are the most similar to one another and which might therefore originate from the same populations. In this paper, we compare these findings to previous works to offer a number of possible interpretations.

2 Data

We made use of the NGC 7023 data from guaranteed time observation program #1192 (GTO-1192) for JWST (Misselt et al. 2025). The observing strategy is similar to that of the PDRs4All program for the Orion Bar (Berné et al. 2022), where NIRSpec and MIRI IFU mosaics spatially resolve the variations among the spectral features. The footprint of the NIRSpec mosaic is illustrated in Fig. 1, with a MIRI F1130W image of the same program displayed in the background. The additional rectangles shown in Fig. 1 represent two regions that had previously been studied using Spitzer data (Berné & Tielens 2012; Mackie et al. 2015); in those studies, the 18.9 µm emission of C₆₀ was detected on the star-facing side of the respective rectangles. The field of view of the JWST spectroscopy does not extend as far toward the central star and, therefore, no detection of the 18.9 µm feature is expected.

A detailed overview of the spectroscopic observing strategy, the data reduction, and a spectral inventory are given in Misselt et al. (2025). In summary, the authors employed a few extra steps including an improved flagging of outliers for a cleaner end result and a world coordinate system (WCS) adjustment based on the imaging data. To apply our all-wavelength, per-spaxel spectral decomposition (Sect. 3.3), we need data cubes whose spatial grids are identical. We built a total of 13 separate data cubes (NIRSpec IFU G395M and MIRI MRS Channels 1 through 4; SHORT, MEDIUM, and LONG bands) by starting from the CRF files produced by Misselt et al. (2025) and running the cube building step of the JWST pipeline with custom arguments. The NIRSpec IFU cubes have the smallest spatial extent, while the MIRI MRS Channel 1 covers a similar area. To ensure full wavelength coverage and to avoid oversampling the longer wavelength channels, each cube was built to the resolution of the MIRI MRS channel 4 (the largest spaxel size, 0.″35), while the analysis was limited to the NIRSpec IFU footprint (3.″3 × 21″ or 9 × 52 spaxels). This combination of spaxel size and footprint dimensions results in around 500 spectral extractions for our analysis. The footprint and the size of the spaxel apertures are illustrated in the right panel of Fig. 1.

To correct the astrometry based on NIRCam data, Misselt et al. (2025) applied a shift of around 6 spaxels to the NIR-Spec data, parallel to the long axis of the mosaic and toward the star. At the NIRSpec data cube resolution (0.″1_), this shift corresponds to around 0.″6_, or about 2 spaxels in our lower res olution (0.″35) cube. To apply the same correction, we modified the WCS of each CRF file prior to applying the cube build step. The shift for the MIRI IFU cubes was minimal. In practice, this WCS correction reduces the gap between the emission peaks at the DF when the NIRSpec and MIRI IFU maps are compared.

We also extracted spectra for two larger apertures: the “ATM” (atomic PDR) and “DF1” (dissociation front 1) regions defined by Misselt et al. (2025) and shown in the right panel of Fig. 1. We used these spectra, characterized by a high signal-to-noise ratio, to illustrate the spectral variations in the top panels of Figs. 2 and 3. We also use them to demonstrate the decomposition method explained in the next section. To apply the decomposition over the complete wavelength range of interest (3.2–26 µm), a single data cube covering the combined NIR-Spec and MIRI wavelength range is needed. The jumps in flux between the channels are small (a few %) and to combine the cubes that are already on a common spatial grid by construction (as described above), we used an algorithm that smoothly blends the spectral segments by using a sliding weighted mean in the overlap regions (Van De Putte et al. 2024).

Fig. 1 Left: position of HD 200775 relative to the NGC 7023 NW filament and footprints of the NIRSpec observations. The background is the F1130W map from the MIRI imager observations of the same program (Misselt et al. 2025). The dashed rectangle indicates the coverage of Spitzer-based spectroscopy and feature maps (Werner et al. 2004; Shannon et al. 2015), which includes coverage closer to the star, where C60 was observed. The labeled footprints are the regions of interest for the studies by Berné & Tielens (2012) and Mackie et al. (2015). The diamond (a 3.6″ square) indicates the typical spectral extraction size for Spitzer-based studies. Right: zoom-in of the same image, showing the two extraction apertures of Misselt et al. (2025) used for the demonstration in Figs. 2 and 3, and the size of a single 0.35″ spaxel of our spectral extractions (Sect. 2). The edge of the mosaic is around 37″ away from the star.

3 Method

3.1 Spectral decomposition with PAHFIT

To decompose the spectra, we made use of the latest version of the PAHFIT tool (Smith et al. 2007). While it was originally written in IDL (interactive data language), PAHFIT was recently redeveloped in Python to enhance its usability for the analysis of JWST data¹. For this work specifically, we used the same version of PAHFIT as presented by Van De Putte et al. (2025). A crucial part of applying PAHFIT to a given spectrum is the use of a suitable “science pack,” which is the list of components to be fit and their parameters. The PAHFIT model represents the AIBs with a combination of Drude profiles for which the amplitudes are fit independently, with the central wavelengths and the full width at half maximum (FWHM) listed in the science pack. To set up a suitable AIB configuration, we started with the parameters provided by the “PDR pack” presented by Van De Putte et al. (2025), which was set up based on JWST spectra of the Orion Bar (Peeters et al. 2024; Chown et al. 2024; Van De Putte et al. 2024). For a detailed account of how the PDR science pack parameters were chosen, we refer to the appendix in Van De Putte et al. (2025). Since the signal-to-noise and AIB-to-continuum ratios of the 16–18 µm complex are much higher in NGC 7023 compared to the Orion Bar, a more detailed view of these emission bands can be revealed. We revised this part of the PDR pack to obtain reliable fits of the 17.4 µm feature (see Sect. 4.2). These changes will be included in an updated version of the PDR pack and each mention of the PDR pack hereafter refers to the one updated in this work.

The PDR pack is supported by additional PAHFIT components that represent the thermal continuum of the dust. The standard PAHFIT approach uses modified blackbody (MBB) curves for which the intensity scales as λ⁻²B_ν(λ, T), where B_ν is the Planck function evaluated in F_ν units (e.g., MJy/sr) at a wavelength λ for a certain temperature, T (Smith et al. 2007). Van De Putte et al. (2025) introduced an alternate continuum shape which includes the silicate emission around 10 and 20 µm, to model the steep 15–26 µm continuum of the Orion Bar. For NGC 7023, we found that using a set of standard PAHFIT MBB of temperatures T_i (intensity ∼ λ⁻²B_ν(λ, T_i)) produces an acceptable continuum model that does not suffer of the mismatch near 15 µm, as the continuum of NGC 7023 is not as steeply increasing toward longer wavelengths. We include continuum components with the following values for T_i: 35, 40, 45, 50, 58, 65, 78, 90, 113, 135, 168, 200, 250, 300, and 400 K. A stellar continuum component was not included. We note that the PAHFIT attenuation model was not included in the fit, even though the molecular region is likely affected by dust extinction, as would be expected given the presence of H₂O and CO₂ ice absorption (Misselt et al. 2025). As noted in previous works (Peeters et al. 2017; Van De Putte et al. 2025), the inclusion of attenuation can lead to degeneracies that cause large jumps in the attenuation between different spatial elements.

Analogous to Van De Putte et al. (2025), we removed the emission lines from the data, rather than including them in the model and fitting all of them with PAHFIT, as this can improve the convergence and consistency. This was done by masking out small parts of the wavelength range, typically within two times the FWHM of each line.

Fig. 2 Zoomed-in view of the emission complexes in the ATM and DF1 regions of NGC 7023 and the corresponding PAHFIT decompositions. The aperture-extracted spectra (defined in Sect. 2 and Fig. 1) are shown in the top panels with a normalization and offset that match the base and peak of a certain feature (first and second columns scaled to 3.3 µm, third column to 5.2 µm). These still have the emission lines included. The gray lines in the bottom two rows are the modified ATM and DF1 spectra (emission lines removed), used as input for PAHFIT. The components of the decomposition are shown as labeled in the legend.

3.2 Integrated feature power definitions

Most bands are decomposed into multiple subcomponents by the PDR pack. The PAHFIT model yields the power of each fitted Drude profile, integrated over the width of the profile in frequency space (surface brightness in W m⁻² sr⁻¹). In the figures and sections that follow, we often refer to the “total integrated power” for the major emission bands. With this, we mean a sum of the aforementioned integrated Drude power, over all sub-components considered part of a band, as defined by Van De Putte et al. (2025). Considering the updates to the PDR pack presented in Sect. 4.2, we provide an updated version of these sum definitions in Table 1.

3.3 Application of PAHFIT to the spectral data cube

The scripts to apply PAHFIT to each spaxel were bundled and made available as the PAHFITcube utility². This freely available tool also keeps track of the input and output data to construct maps, and offers functions for parallel processing, resuming a series of fits after an interruption, and making overviews of the maps in a similar way to what is shown in Fig. A.2.

3.4 Limitations

This paper is an initial overview of the carbonaceous emission bands in our NGC 7023 data. Here, we address a number of issues that may need to be addressed in followup work. Since this work applies the PDR pack with PAHFIT, all the same caveats listed by Van De Putte et al. (2025) are relevant here as well. In addition, we highlight a few subjects that might require a more careful or dedicated analysis.

First, we note a few wavelength ranges where one should be careful when interpreting the power of the features as determined by these PAHFIT results. For the weaker features in the 3.2–3.6 µm range, the fit is affected by the lack of description of the nebular continuum, as was noted in the section on Limitations in Van De Putte et al. (2025). For the 3.3 µm band, it is notable that it does not seem to widen toward the blue side as was observed for the Orion Bar spectra, but this effect is subtle and a good continuum model is needed to quantify it. For similar reasons, there is a degeneracy between the assumed FWHM of the broader features in the 3.8–4.7 µm range (Table 2) and the assumed continuum level. Typically, a continuum estimation that is too low will be compensated by a tuning with wider features. The 6.0 µm feature at the base of the blue wing of the 6.2 µm band will be underestimated if there is a mismatch between the model and this blue wing. For this work, the FWHM of the 8.6 µm band was adjusted to 0.3 µm, though this band also appears to have an additional narrow component, as can be seen in the ATM spectrum (Fig. 3), while the PDR pack approaches this with just one broad feature. The 11.0 µm feature is positioned on the steep blue wing of the 11.3 µm band, while the 11.0 µm feature itself also likely has a steep blue wing and a wider red wing. For the Orion Bar, the 11.0 µm profile and amplitude are sufficiently approximated by a single Drude profile (Van De Putte et al. 2025), but for NGC 7023, this feature is much stronger relative to 11.3 µm. It seems that multiple components are needed to reproduce the shape more closely (Shannon et al. 2016). The smaller features redward of 13.5 (e.g., 14.2 µm) tend to be underestimated due to the influence of the wide component added to fit the red side of the 11–14 µm complex (centered around 13.2 µm), which is prominent in the DF1 spectrum. Finally, the feature at 17.7 µm appears rather complex, with a flat peak and some substructure, and it is approximated by a single Drude profile.

To interpret the region beyond the DF1, the presence of ice absorption needs to be taken into account. For the 3.3 µm band, the H₂O absorption spanning 3.2–3.3 µm (Misselt et al. 2025) will have an effect. Additional ice absorption bands are found at 4.2–4.3 µm (CO₂ stretching), at 14.8–15.4 µm (CO₂ bending), and at 4.64–4.69 µm (CO). For this reason, this region (MO1, DF2 and MO2 in Misselt et al. 2025) is not the main focus of this work. In addition, the presence of this ice indicates that a dust attenuation correction will be necessary to properly compare the inter-band ratios across the mosaic.

Finally, we note that this work does not make full use of the available resolution, since the data cubes were built on a common resolution matching MRS channel 4 (0.″35, see Sect. 2). The phenomena we demonstrate are still clear with the larger spatial scales probed by our choice of spatial grid. For the pruposes of a future work that focuses on a specific wavelength range, data cubes that exploit the finer spatial resolution of the shorter channels could be built. This will require point spread function (PSF) matching or interpolating appropriately. Alternatively, the tools for fitting entire spectral cubes may need some form of regularization to address the difference in spatial resolution across the channels.

Fig. 3 Analogous to Fig. 2, for 6–13 µm and 13–23 µm. The normalization in the top panels is based on the 7.7 and 16.4 µm features, respectively. The updated tuning for the 16–18 µm complex better captures the shape changes, resulting in a more robust extraction of the power of the 16.4 and 17.4 µm features.

4 Results

We focus our discussion on the brightest emission complexes with the strongest profile variations and show three views of the results throughout this section. We start with a full-spectrum overview of the emission profile differences and how these are captured by the PAHFIT decomposition in Sects. 4.1 and 4.2 (Figs. 2 and 3). We then discuss the spatial distributions of the main complexes (Fig. 5) and their subcomponents (Figs. 7–9) in Sect. 4.3.

Fig. 4 Detailed view of the updated decomposition in the 16–18 µm range. The light gray line shows the ATM spectrum, with bumps at 16.7 and 17.2 µm that raise the need for three components to model the central “plateau.” The individual Drude profiles (colorful curves and labels) are plotted after adding the continuum components, and the total model is the solid black line.

4.1 Main bands

Before discussing the spatial variations seen among the emission profiles and the corresponding spatial distributions, we provide a brief overview of common interpretations for each complex. The main features in the NIRSpec range are the 3.3 µm aromatic C-H stretching and 3.4 µm aliphatic C-H stretching features (Fig. 2, left panels). As already shown by Misselt et al. (2025), the 3.4/3.3 µm ratio that traces the aliphatic material abundance is higher for the DF, analogous to what a similar decomposition reveals for the Orion Bar (Peeters et al. 2024). Two weaker features between 3.8 and 4.8 µm (Fig. 2, center panels) were attributed to analogous stretching modes of deuterium-substituted material: the deuterated 4.4 µm aromatic and 4.7 µm aliphatic features (Misselt et al. 2025). The identification of the weaker features at 5.2, 5.7, 5.9, and 6.0 µm (Fig. 2, right panels) has not been established as precisely. Previously, Chown et al. (2024) and Van De Putte et al. (2025) identified and quantified the shape changes of the 5.7 µm feature between the Atomic PDR and DF regions of the Orion Bar. Through overtones or combination modes, a connection was proposed with the aromatic C-H out-of-plane bending modes in the 10–15 µm range (Allamandola et al. 1989; Boersma et al. 2009; Mackie et al. 2015; Chown et al. 2024).

The wavelength range shown in Fig. 3 contains a number of major complexes. The bands near 6.2 µm and 8.6 µm are typically enhanced for PAH cations and, hence, they are considered to be tracers of cationic aromatic species (Allamandola et al. 1999; Peeters et al. 2002). The 6.2 µm band is associated with aromatic C-C stretching, and the 8.6 µm band with aromatic C– H in-plane bending. The 7.7 µm complex, exhibits narrow bands attributed to mixed-character aromatic C-C stretching or C–H in-plane bending modes, also associated with cations. These bands appear to be perched on top of a broader component near 7.9 µm, which has previously been associated with the presence of PAH clusters or very small grains (VSGs; Berné et al. 2007; Pilleri et al. 2012; Peeters et al. 2017; Stock & Peeters 2017; Khan et al. 2025). The decomposition with the PDR pack attempts to separate these narrow and broad contributions, while the broad-to-total ratio determined based on the PAHFIT results was shown to increase with the distance from the star in the Orion Bar (Van De Putte et al. 2025).

The 10–15 µm range, contains various aromatic C–H out-of-plane bending modes, with the strongest peaks near 11.3 and 12.7 µm. The total 11.3 µm flux is often used in ratios (such as 6.2/11.3 µm) to trace neutral PAHs, while there are components associated with cationic species as well, such as the 11.0 µm feature (e.g., Hudgins & Allamandola 1999) and a blue-side component of the 11.3 µm band identified as 11.207 µm by Khan et al. (2025). We note that Fig. 3 shows that the 11.0 µm feature is significantly stronger in the ATM spectrum. These modes can also be interpreted in terms of hydrogen adjacency at the edges of the aromatic material (solo, duo, trio, quartet; e.g., Hony et al. 2001; Bauschlicher et al. 2008; Khan et al. 2025). For example, the 11.3 µm band is dominated by solo hydrogens, and the 12.7 µm band by duo and trio hydrogens.

The group of bands in the 16–18 µm range consists of aromatic C–C–C bending modes, also known as skeletal modes (Allamandola et al. 1989; Van Kerckhoven et al. 2000). These likely contain information about the structure on a larger scale, since they involve the bending of larger groups of atoms (Peeters et al. 2004; Boersma et al. 2010), and the ratios of the components have connections to the charge state as well (Peeters et al. 2012; Shannon et al. 2015). The spectra in the top right panel of Fig. 3 have been normalized to the 16.4 µm band, which emphasizes the difference in the relative 17.4 µm feature strength.

Fig. 5 Power maps of the main bands based on sums of the subcomponents contributing to each band as defined in Table 1. The color scale represents the spectrally integrated power (W m−2 sr−1) normalized to the brightest spaxel. The dashed black vertical line corresponds to DF1 (more precisely: the maximum of 16.4 µm when averaged vertically) and the ATM region is in the leftmost 3″ of the image.

Fig. 6 Maps of the individual components of the 5.7 µm band decomposition. The lower panel is analogous to Fig. 5. The top panel is a legend illustrating the subcomponents labeled A, B, and C. The light gray line shows the ATM spectrum data, the solid black line is the best fitting total model, and the dashed line shows the sum of all other components (continuum and wings of nearby Drude profiles), to which the three 5.7 µm components were added.

4.2 Decomposition and PDR pack updates

To emphasize the profile variations, the extracted spectra for the ATM and DF1 regions are shown in the top panels of Figs. 2 and 3. They were normalized by applying a vertical shift and a normalization factor using reference points at the peak and the base of a particular feature (as stated in each caption). Below, we describe a number of updates that were made to the PDR pack, based on the details and differences between ATM and DF1 that are revealed by this comparison. A detailed description of the decomposition as specified by the PDR pack is provided in the appendix of Van De Putte et al. (2025) and not included here for brevity.

The 16.4 µm profile is similarly shaped for the ATM and DF1 spectra of NGC 7023, while the rest of the 16–18 µm complex exhibits significant differences. Compared to the Orion Bar spectra of Van De Putte et al. (2025), these differences are more strongly pronounced in NGC 7023 owing to its fainter continuum. For this reason, the decomposition by the original PDR pack could not sufficiently capture the flux of the 17.4 µm band and the details of the emission profile variations. Therefore, we updated the PDR pack for this wavelength range using a procedure analogous to Van De Putte et al. (2025). A set of central wavelengths and FWHM values were determined for both the ATM and DF1 spectrum by applying a least-square fit to the 15.7–18.0 µm range, and a continuum subtraction based on a preliminary PAHFIT model. The results for this update are summarized in Table 2. The main difference is that the total number of Drude profiles has increased by 2, to better describe the shape of the plateau between the 16.4 and 17.4 µm bands. The new features were motivated based on two bumps seen in the ATM spectrum (Fig. 4): one on the red side of the 16.4 µm feature, and one near the blue side of the 17.4 µm feature. In the DF1 spectrum, the plateau peaks near 17.0 µm instead, indicating a third feature. The version of the PDR pack based on the Orion Bar data could not sufficiently capture these details, as it defined only two wide components centered at 16.8 and 17.1 µm to model this plateau.

These data also provide an opportunity to determine a set of Drude profiles near 3.8–5.0 µm for the PDR pack. This wavelength range is problematic in the Orion Bar data, because of the wavelength gap in the G395H observation and the larger number of emission lines. Our tuning for the features near 3.8, 4.0, 4.35, and 4.75 µm is given in Table 2. We note that the features in Table 2 are likely too wide, as their tuning partially compensates for the lack of a continuum description (see ‘Limitations’ section of Van De Putte et al. 2025). A realistic continuum description is challenging as there are possible contributions from nebular emission or PAH-related quasi-continua (Allamandola et al. 1989; Boersma et al. 2023).

The updated PDR pack was applied to the full 3.2–26 µm range and two resulting fits are shown in the lower two panels of Figs. 2 and 3. As stated in Sect. 3.1, we did not use the PAH-FIT dust attenuation model for the final results. As a test, we performed the decomposition both with and without the classic PAHFIT dust attenuation model (Smith et al. 2007). The attenuation resulting from the fit was zero for both the ATM and DF1 spectra. We note that Fleming et al. (2010) applied IDL PAHFIT to Spitzer spectra of NGC 7023, and found a negligible optical depth at the peak of the silicate extinction feature at 9.7 µm (τ_9.7 < 0.2). Past the DF, some level of extinction is expected, considering that ice absorption is seen in the NIRSpec data (Misselt et al. 2025). Therefore, as stated in Sect. 3.4, this region is not the main focus in this work.

Fig. 7 Maps of the 7.7 µm decomposition and legend analogous to Fig. 6. The sum of the narrow components (as in Table 1) is shown in addition to the individual maps.

Fig. 8 Analogous to Fig. 6 for the 11.0 and 11.3 µm bands.

Fig. 9 Analogous to Fig. 6 for the 12.7 µm band. Note the similarity of 12.7A to 11.3A in Fig. 8.

4.3 Spatial distributions

A selection of spatial distributions resulting from our method (Sect. 3.3) is displayed in Figs. 5 through 9. A complete overview of maps for all individual PDR pack components is provided in Appendix A. We limited the scope of the discussion to the bands at 3.3, 3.4, 5.2, 5.7, 6.2, 7.7, 8.6, 11.0, 11.3, 12.7, 16.4 and 17.4 µm. While the PDR pack also fits numerous weaker features, those fits are more sensitive to being biased by a mismatch in the continuum model and such features require an individualized extraction approach.

One of the main points of this discussion, is how individual subcomponents have spatial distributions that differ from the integrated power in each band (Sect. 3.2, Table 1). We start with an overview of the total power maps in this section. In what follows for the 5.7, 7.7, 11.3, and 12.7 µm bands, Figs. 6, 7, 8, and 9 zoom in on the subcomponents and show maps of the individual components alongside the corresponding total (see Sect. 4.3.2). In the following sections, we describe the similarities and differences observed between these spatial distributions, further clarified by the one-dimensional cuts shown in Fig. 10. For each cut, the data point are medians over the vertical range shown in the bottom panel. We set the lower limit of the vertical axis to zero to emphasize that all emission maps are nonzero on both sides of DF1 and we normalized each one-dimensional spatial profile to its peak brightness.

4.3.1 High and low emission in the ATM region

A common element among all maps in Fig.5 is a maximum at DF1, with the 11.0 µm band being the only exception. In the figures that follow as well as the appendix (Fig. A.2), several other maps lacking the DF1 peak are seen: the smooth decrease of 11.0 µm with no trace of the DF1 peak also appears for 10.8 µm (Fig. A.2). Other maps lacking the DF1 peak do not show the same smooth behavior; they are noisier (e.g., 5.44), more biased by the continuum fit (e.g., 14.0), or suffer from degeneracy with another feature (e.g., 7.7C, 7.7D). In addition, they show a region with zero power near DF1 as a result.

The most conspicuous difference between the maps is found in the star-facing half of the field of view (which contains ATM): a wide variety of emission strengths relative to the DF1 peak. Hence, we categorized the features into three types based on their relative brightness in the ATM region compared to their brightness at DF1. The ratio P_ATM/P_DF1 can be low (<0.5, type I), intermediate (>0.5, <0.8 type II), or high (>0.8, type III). While there are more details that characterize the spatial profiles, this grouping into types based on the ATM/DF1 intensity ratio will aid the discussion. Considering the total band emission shown in Fig. 5, and the cuts in Fig. 10, the bands at 3.3, 3.4, 5.2, 5.7, 7.7, and 11.3 µm are of type I, those at 6.2, 7.7, 8.6, 12.7, and 16.4 µm are of type II, while the 11.0 and 17.4 µm features are of type III. Spatial cuts for several type I features are grouped in the top panel of Fig. 10. The 3.3 and 11.3 µm (total) cuts are nearly identical and have matching peak positions despite being observed with different instruments. The center panel shows several type II spatial profiles. For the spatial cuts of 6.2 and 8.6 µm, we note the parallel nature of the curves on the ATM side, despite the intensity difference in that region.

We emphasize the difference between the 16.4 and 17.4 µm bands. The 17.4 µm cut is shown in both the type II and type III panels of Fig. 10 to facilitate a comparison to the 16.4 µm cut. The 16.4 and 17.4 µm cuts have matching peak positions and behave similarly beyond DF1. The spatial profiles are nearly parallel in the ATM region as well, though the 17.4 µm feature exhibits a secondary maximum nearly as strong as the DF1 peak, while the strength of 16.4 µm is suppressed relative to DF1. The nature of these two features is discussed further in Sect. 5.1.

Fig. 10 Spatial profiles along cuts perpendicular to the DF for selected features, normalized to the peak of each cut. The cuts are organized by spatial distribution types I, II, and III, as introduced in Sect. 4.3. The 17.4 µm cut is shown twice to directly compare it to the type II and type III cuts. The dashed lines in the bottom panel indicate the region over which the maps were median-stacked along the vertical axis to obtain the cuts, with the 17.4 µm map from Fig. 5 in the background. The shaded areas indicate the widths of the ATM and DF1 regions.

4.3.2 Individual subcomponents

The top panels of Figs. 2 and 3 show the band profile differences between ATM and DF1 and these variations are similar to those seen in the Orion Bar (Chown et al. 2024). The fits with the PDR pack capture some of these profile differences if the relative contributions of the subcomponents to each band are considered. The maps of the individual subcomponents of the 5.7, 7.7, 11.3, and 12.7 µm complexes (Figs. 6–9) reveal additional spatial distributions of type III (bright in ATM). These figures include a panel with a detailed view of the decomposition, showing the ATM spectrum and its PAHFIT model to clarify the interpretation of the subcomponent maps.

For the 7.7 µm maps (Fig. 7), we have added maps of two partial power sums: “7.7 narrow” and “7.9 broad” (Table 1, Sect. 3.2). A spatial cut of the “7.7 narrow” combination is shown in the lower (type III) panel of Fig. 10 and it is remarkably similar to the 11.0 µm feature. As noted in the previous section, the 11.0 µm cut has less substructure; there is no peak at DF1. The top panel of Fig. 8 shows that the 11.0 µm feature is not very well approximated by a single Drude profile, making it likely that the 11.0 µm profile has a pronounced red tail similar to 11.3 µm. We verified the 11.0 µm spatial distribution result by applying an ad-hoc 11.0 µm measurement that consists of fitting a local linear continuum and numerically integrating the data. This approach results in a nearly identical spatial distribution.

The 11.3 and 12.7 µm bands both exhibit profile changes, and the PDR pack explains this as their bluest subcomponents being of type III. The 11.3 µm band is described by four subcomponents (A, B, C, and D) all of which are shown in Fig. 8. The peak of the 12.7 µm band is decomposed into two features, although the blue wing contains a significant contribution by the 12.6 µm component (Fig. 9). We show cuts for the bluest subcomponents called 11.3A and 12.7A (see also Table 1) in the type III panel of Fig. 10, showing the similarity with the 17.4 µm band. Analogous behavior seems to be present for the blue component of the 5.7 µm band, though the spectra for this band features and the resulting maps of the components are noisier (Fig. 6). The cuts are also compared to 11.0 µm, and we note how 11.3A and 12.7A still have more substructure than 11.0 µm. Therefore, 11.0 µm could be considered a special or extreme case of the type III spatial distribution.

Conversely the reddest subcomponents tend to exhibit type I spatial distributions, being suppressed in the ATM region and having a strong peak at DF1. The clearest examples are 5.7C (Fig. 6), 7.7 broad (Fig. 7), 11.3D (Fig. 8) and 12.7B (Fig. 9). The 7.7 broad spatial cut was added to the type I (top) panel of Fig. 10, showing its similarity to the 3.4 µm aliphatic C-H feature, with nearly identical peak positions and widths. To further emphasize the similarity of the spatial behavior of the profile variations at 5.7, 7.7, 11.3, and 12.7 µm, we show subcomponent-to-total ratio maps (e.g., 5.7A/5.7) in Fig. 11. We note how these values show a positive gradient (i.e., the profiles become more blue-side dominated) in the direction toward the star, that persists up to the left edge of the field of view. We propose by extrapolation that the profiles keep evolving beyond this boundary, for regions progressively closer to the star (see also Sect. 5.5).

Summarizing, the subcomponent maps reveal that the bluest subcomponents at 5.7, 7.7, 11.3, and 12.7 µm have type III spatial distributions, even though the total band emission is of type I or type II. We provide an overview of the spatial distribution types of the features in Table 3, together with the P_ATM/P_DF1 ratio, in which the feature strengths P_ATM and P_DF1 are measured based on the fits shown in Figs. 2 and 3, and the strengths of summed features are computed as defined in Table 1.

Fig. 11 Maps of ratios describing variations in the band profiles. The color bars represent a linear scale between the minimum and maximum. The ratios represent the relative contributions of subcomponents to the total as defined in Table 1. The ratios of the bluest components of each complex are shown alongside the 17.4/16.4 ratio to illustrate their similarity, and how the steepest gradient is in the same location.

5 Discussion

There are numerous previous studies that mapped the carbonaceous emission bands in NGC 7023 using Spitzer data (Werner et al. 2004), cited throughout the discussion that follows. With the typical 1.″8 resolution of Spitzer-based spectral maps, details such as the structure and the precise location of the DF could not be resolved. The spatial resolution of our JWST-based maps reveals that most of the changes in the carbonaceous band characteristics occur near DF1 or at around 3″ ahead of this DF (see ratio maps of Fig. 11). In the following sections, we start by discussing the differences between the spatial distributions of the 16.4 and 17.4 µm features. This is followed by a discussion of possible relations with tracers based on the other bands and interpretations from Spitzer-based work.

5.1 Spatial distribution of 16.4 and 17.4 µm

Originally, Sellgren et al. (2010) reported that the 16.4 and 17.4 µm spatial distributions were similar in NGC 7023. In later works, it is typically indicated that two spatial distributions exist for this emission complex, one of which extends further toward the star (Rosenberg et al. 2011; Boersma et al. 2013, 2014, 2015; Shannon et al. 2015; Sidhu et al. 2023). For example, Shannon et al. (2015) discussed the morphology of the 16.4 and 17.4 µm features, and noted how maps of the 17.4 and 11.0 µm emission are both brighter close to the star, and that maps of the 12.7 and 16.4 µm emission are similar as well. Our maps and cuts presented in Sect. 4.3 clarify that the 16.4 and 17.4 µm bands have highly similar spatial distributions on the far side of the DF, while in the direction of the ATM region, the 17.4 µm emission is enhanced (Figs. 5 and 10). The high resolution of JWST reveals that a change in the 17.4/16.4 ratio takes place near the boundary between DF1 and the more irradiated atomic PDR region (Fig. 11). It is therefore likely that UV processing plays a role, where the resulting changes to the emitting population (molecular size, charge, or structure) lead to different contributions by the relevant vibrational modes.

As we discuss in Sect. 4.1, the vibrational modes linked to these bands are skeletal C–C–C modes, as revealed by computed spectra of large PAHs (Allamandola et al. 1989; Van Kerckhoven et al. 2000). More specifically, the computational results of Ricca et al. (2010) show contributions at 16.4 µm by in-plane elongation and compression modes, along with contributions to 17.4 µm by modes involving C–H out-of-plane bending. Based on spectra in the NASA Ames PAH database (PAHdb), the presence of pendent rings was suggested as an explanation for the 16.4 µm band (Boersma et al. 2010), while larger (but compact) PAHs (Ricca et al. 2012) in the 50–130 carbon atoms range consistently show an emission band near 17.4 µm.

Peeters et al. (2012) organized the components contributing to the 16–18 µm emission into three groups, based on the correlation between the bands as observed in NGC 2023. The 16–18 µm plateau and 15.8 µm band correspond best to the 11.3 µm emission, the 16.4 µm band closely matches the 12.7 µm emission, and the 17.4 µm emission appears to represent a more spatially diffuse component. Peeters et al. (2012) associated the 17.4 µm band to doubly ionized PAHs or a certain subset of cationic PAHs. Previously, it was concluded that the (normalized) 16– 18 µm bands are mostly uncorrelated among themselves and that (besides the 12.7–16.4 µm connection) there is no obvious correlation with the MIR features at 7–9 µm and 11–15 µm (e.g., Boersma et al. 2010). The JWST spectra, our decomposition, and the spatial distributions of the individual subcomponents reveal the opposite result (Sect. 4.3.2). In particular, it appears that the 17.4 µm spatial distribution is particularly similar to the blue-side and narrower components of the 7.7, 11.3, and 12.7 µm bands, as we demonstrate via correlation diagrams in the next section.

Fig. 12 Correlation plots exploring the relation between the 16.4 and 17.4 µm bands and a selection of other bands. The data points represent the spaxels of the maps, and their colors represent the offset from the edge of the mosaic, with the locations of ATM and DF1 indicated on the color bar. The black lines show linear fits through the origin (y = ax, solid) and with a variable intercept (y = bx + c, dashed). These lines are typically closer for 17.4 µm, and more separated for 16.4 µm, with the exception of the 12.7 µm panels. This indicates that the 7.7 narrow, 8.6, 11.3A, and 12.7A features relate better with 17.4 µm, while the total 12.7 µm flux relates beter with 16.4 µm (Sect. 5.2).

5.2 Correlations of 16.4 and 17.4 µm with feature ratios

In the previous section, we note that it was unclear which features were best correlated with 17.4 µm (Boersma et al. 2010), while the spatial position of the peak emission in NGC 2023 indicates that the 17.4 µm band is grouped with the 8.6 and 11.0 µm bands, and the 16.4 µm band with the 6.2, 7.7, and 12.7 µm bands (Peeters et al. 2017). Our results, presented in Sect. 4.3, reveal a visual similarity between the spatial distributions of the 17.4 µm band and the narrower bluer subcomponents of the 7.7, 11.3, and 12.7 µm bands, as well as the 11.0 µm band.

In this section, we discuss the correlation diagrams (Fig. 12) comparing the 16.4 and 17.4 µm bands to the aforementioned bands and their blue subcomponents. All quantities shown have been normalized to the total 11.3 µm band. Although this point is not discussed in this work (see Sect. 3.4), we included the molecular region (“Mol”) in the diagrams for completeness (shown on the right side of DF1 in the maps). As a visual aid, we added two linear fits to each panel of Fig. 12: one which is a simple scaling relation through the origin (y = ax), and one which allows a non-zero intercept (y = bx + c). When these linear fits are parallel, then extrapolating the relationship predicts a behavior where both quantities approach zero simultaneously, indicating that the emitting species are cospatial and behave as a single population with no other contributions. Conversely, a larger separation between the linear fits visualizes a non-zero intercept. Even if the correlation is strong, this means that at least one of the two bands in the diagram has significant contributions from another population with a different spatial distribution.

For the total 8.6, 11.0, and 12.7 µm bands, we find a result that is consistent with the grouping of Peeters et al. (2017): the 17.4 µm feature relates best to the 8.6 and 11.0 µm features, while the 16.4 µm band appears to relate best to the total of 12.7 µm. We conclude this from the linear fits of Fig. 12: in each quadrant, compare the lines for 16.4 µm (top row) and 17.4 µm (bottom row). In most cases, the lines are closer together for the 17.4 µm diagram, while the total 12.7 µm band is the exception, establishing its relation with the 16.4 µm band instead. The total 7.7/11.3 ratio shows a strong bifurcation in the diagram, because of the varying contribution by the 7.7 broad component, which gets stronger on the shielded side of the DF. This becomes even more apparent if the molecular region is included (blue points). This indicates that both the 16.4 and 17.4 µm features are likely not related to the broad 7.7 µm component.

Inspecting the subcomponents at 11.3 and 12.7 µm reveals additional relationships with the 17.4 µm band. While the total 11.3 µm emission cannot be shown as it is the normalization, we can still show the blue subcomponent ratio 11.3A/11.3 here, which describes the band profile. Large 11.3A values correspond to a sharper profile, similar to that of the “atomic PDR” template spectrum of the Orion Bar (Chown et al. 2024). There is some difference between the two linear fits for 16.4 µm, while the lines almost overlap for 17.4 µm. This shows that while both 16.4 and 17.4 µm correlate strongly with 11.3A and 12.7A, the 17.4 µm band is co-spatial with 11.3A while the 16.4 µm band is not. For the 12.7 µm band, considering the blue 12.7A subcomponent individually has a major effect on the relationship with the 17.4 µm band. The 12.7A feature relates to 17.4 µm analogously to 11.3A, despite the fact that the total 12.7 µm emission relates best to the 16.4 µm band.

These somewhat subtle differences in the correlation diagram originate from differences that are more obvious when the spatial maps (Fig. 11) or the cuts (Fig. 10) are inspected visually, emphasizing the importance of these spatially resolved studies. Summarizing, the diagrams of Fig. 12 revealed a key difference between the 16.4 and 17.4 µm bands: 17.4 µm emission is co-spatial with the 7.7 narrow, 8.6, and 11.0 features and the blue 11.3 and 12.7 µm subcomponents, while only part of the contributions to the 16.4 µm band come from the same population. In the following, we discuss possible interpretations of these connections.

5.3 Charge tracers

Several of the features highlighted in the previous section have the following in common: they are considered tracers of charged PAHs. The 6.2/11.3, 7.7/11.3, 8.6/11.3, and 11.0/11.3 ratios are considered key quantities to trace the ionization state of PAHs, where it is assumed that the 11.3 µm emission is dominated by emission from neutrals, while the 6.2, 7.7, 8.6, and 11.0 bands are dominated by cationic PAHs (Hudgins & Allamandola 1999; Hony et al. 2001). These tracers were used in the past to interpret the changing band ratios across NGC 7023 in terms of the charge state (Shannon et al. 2015; Stock et al. 2016; Boersma et al. 2014; Peeters et al. 2017; Sidhu et al. 2022). In Fig. 12, the region including ATM up to DF1 appears in the top right of the correlation diagrams, indicating that charged emission carriers are more abundant, as expected from previous work. Our maps (Figs. 5–9) and spatial profiles (Fig. 10) paint a spatially resolved picture about the contribution of charged and neutral carriers to the different bands. We suggest that the type III features are dominated by a charged population that itself has a type III spatial distribution, while the neutral population is distributed according to type I. The mixed (type II) features then contain contributions by both charged and neutral carriers in varying ratios.

Previously, Boersma et al. (2014) revealed that the spectra can be grouped into morphological zones that run parallel to the NGC 7023 DF, and concluded that both cationic and neutral PAHs contribute significantly to the 7.7 and 12.7 µm bands. Other works include blind signal separation techniques, which revealed three representative spectra and spatial distributions in NGC 7023 (Berné et al. 2007; Rosenberg et al. 2011), which the authors associated with charged PAHs, neutral PAHs, and VSGs. Boersma et al. (2015) discuss the use of the “traditional” charge tracers and concluded that the 7.7/11.3 ratio is not a reliable charge tracer near the DF. Our work shows that using the 7.7 µm band as a charge tracer requires carefully separating the narrower (type III) subcomponents from the broad (type I) 7.7 µm emission. Studies of the effectiveness of certain band ratios as charge tracers, such as the one by Boersma et al. (2015), could be revisited while taking the individual subcomponents under consideration. An in-depth comparison of charge tracers may also reveal the selective destruction of certain subpopulations within the cation population (Singh et al. 2025).

Via the correlation with 11.0 µm, Shannon et al. (2015) attributed the 17.4 µm feature to cationic PAHs and suggested that the 16.4 band arises from both cationic and neutral PAHs. Our results support that the 16.4 µm band is of a more mixed nature (type II) than the 17.4 µm feature (type III) and we showed that the latter was more strongly correlated to the ionized PAH tracers (Fig. 12). Shannon et al. (2016) studied the 11.3 and 12.7 µm profile variations in NGC 7023 and NGC 2023, reporting a cationic contribution to 11.3 µm. Similarly, they found a potentially neutral contribution to the 11.0 µm band by decomposing it into two components. Their two components for the 12.7 µm band have different spatial distributions as well, and one of them is a sharper peak on the blue side. Our JWST-based decomposition and maps now resolve these spatial distribution differences: the 12.7A and 12.7B components have type III and type I spatial distributions respectively, with most of the profile variation taking place over a distance of less than 1″ as revealed by the 12.7A/12.7 ratio map (Fig. 11). A similar effect is found for the 11.3 µm band: while the total 11.3 µm emission is often used to trace neutral PAHs, our work reveals that the 11.3A component contributes up to 20% to this band.

5.4 Profile diagnostics and comparison to the Orion Bar

The PDRs4All program (Berné et al. 2022, ERS-1288) observed the Orion Bar PDR with a similar observing strategy as our NGC 7023 observations. Several previous works analyzing the PDRs4All spectra made use of the “template spectra” that were extracted from five apertures representing the key zones in this PDR (Peeters et al. 2024; Chown et al. 2024; Van De Putte et al. 2025, available via the CDS). They include the ionized region ahead of the IF (named “HII”, but it also has a PDR in the background), the neutral region between the IF and the DF (“atomic PDR”), and three DF (“DF1”, “DF2”, “DF3”). Van De Putte et al. (2025) applied the PDR pack to these five spectra and discussed how PAHFIT captures a number of band profile variations which can serve as diagnostics for the photochemical evolution of the emission carriers. One of the results includes a diagram comparing the 6.2 µm FWHM to the 7.7 broad/7.7 total ratios, showing that both of these increase for the more shielded DF1, DF2, and DF3 regions. This was interpreted as follows: the 6.2 µm profile, associated with pure aromatic C–C stretching modes and enhanced for cationic species (Allamandola et al. 1999; Peeters et al. 2002), broadens due to the effects of anharmonicity (Mackie et al. 2022) when small PAHs with higher excitation temperatures contribute to the emission. The 7.7 broad contribution to the total 7.7 µm emission is thought to originate from a population of VSGs (Berné et al. 2007; Pilleri et al. 2012; Peeters et al. 2017; Stock & Peeters 2017; Khan et al. 2025). Both of these components are expected to be vulnerable to the UV radiation field and, therefore, they are enhanced in the shielded regions.

We reproduced this diagram in Fig. 13 and applied the same diagnostics to the NGC 7023 spectra, to clarify how the observed profiles and conditions compare to each of the Orion Bar template regions. The 6.2 µm FWHMs were calculated analogously to Van De Putte et al. (2025), by evaluating the 6.2 µm profile of the PAHFIT model (sum of Drude profiles, Table 1) and determining the maximum and half-maximum points numerically. The Orion Bar templates were ordered in the following way for all diagnostics examined by Van De Putte et al. (2025): atomic PDR, DF1, HII (background PDR), DF2, and DF3. The point cloud for NGC 7023 in Fig. 13 is color-coded by the distance from the edge of the mosaic, and a similar sequence is found, where DF1 and the region beyond show enhanced 6.2 FWHM and 7.7 broad/7.7 total values. The DF1 spectra of NGC 7023 (gray dots as marked with the black line on the color bar) appear to be most similar to the conditions of DF2 in the Orion Bar. The narrowing of the 6.2 µm band toward the star tends to the same FWHM value as in the Orion Bar.

In our decomposition, the 7.7 µm broad feature has a strikingly similar spatial distribution to that of the red wing of the 11.3 µm band (11.3D, Fig. 8). Three spectral signals were extracted over the 10.0–19.5 µm range from NGC 7023 Spitzer spectra by Rosenberg et al. (2011), and they were characterized using PAHdb. They showed a first component with a pronounced 11.0 µm band, a second component with a pronounced but narrower 11.3 µm band, and a third component with a broader red wing at 11.3 µm, interpreted as PAH cations, neutral PAHs, and VSGs respectively. For the Orion Bar, Khan et al. (2025) decomposed the 11.3 µm band into two profiles extracted from the data, both of which have an asymmetric nature, and find that the red wing of the 11.3 µm band behaves more independently from the main PAH emission bands, suggesting VSGs or PAH clusters as a possible origin. Boersma et al. (2014) include the following candidates for the red wing of 11.3 µm: PAH anions, anharmonicity of hot PAHs, VSGs and PAH clusters, or potentially superhydrogenation near DF1, which can play a role in the H₂ formation taking place (Rauls & Hornekær 2008; Thrower et al. 2012; Jensen et al. 2019). Another suggestion for the shift of the 7.7 µm band (here facilitated by the broad contribution) is the presence of aliphatic side groups, which are expected to be removed when the UV radiation field increases (Bregman & Temi 2005; Sloan et al. 2007; Sadjadi et al. 2015a; Shannon & Boersma 2019). That could explain the spatial similarity of 7.7 broad with the 3.4 µm band (aliphatic C-H stretch, Fig. 10), which was previously demonstrated for NGC 7023 by combining the Spitzer data with AKARI spectroscopy (Pilleri et al. 2015).

The regions close to the star exhibit much lower values of 7.7 broad/7.7 total; even when compared to the Orion Bar atomic PDR (Fig. 13). The narrow features continue getting more pronounced as the star is approached, while the broad component is further suppressed. The ratios of 11.3A/11.3, 12.7A/12.7, and 17.4/16.4 in Fig. 11 exhibit a similar gradient, and this cohesive behavior between the band profiles suggests that there are at least two populations, one which becomes gradually more dominant closer to the star. We also note that the 11.0 µm band in the NGC 7023 ATM region (Fig. 8) is stronger relative to 11.3 µm, compared to what is seen in the Orion Bar (Khan et al. 2025). All of the above shows that the differences in the emitting populations are more pronounced than what is observed in the Orion Bar, and we speculate about the cause as follows.

The harder radiation field driving the Orion Bar PDR results in the formation of an IF, which is not the case for the softer radiation field in NGC 7023. This may result in different behavior of the PAH size and charge as a function of the radiation field, for which the intensity is determined by the distance to the star or the shielding by the gas and dust when entering the PDR. The harder radiation field will dehydrogenate and destroy PAHs and ionize the surviving population, which may lead to a situation where only larger species survive in the region with a significant ionized PAH fraction. The PAH database fitting approach by Maragkoudakis et al. (2026) finds that only large PAHs (>70 C atoms) survive near the IF of the Orion Bar. While the estimated PAH size in NGC 7023 also increases closer to the star (e.g., Croiset et al. 2016), the softer radiation field may allow a substantial population of smaller PAH to survive in the region where PAHs are ionized, and such smaller PAH are dehydrogenated more quickly (Andrews et al. 2015). The presence of these smaller PAHs may lead to a different response of the charged/neutral PAH tracers (such as the narrow 7.7 µm features or larger 11.0/11.3 ratio of NGC 7023).

Besides the radiation field, the hydrogen density will also affect the dehydrogenation balance, as higher densities restore the hydrogenation at a faster rate. The different regions in NGC 7023 exhibit a broad range of densities, from the order of 10² cm⁻³ in the central cavity, to 10⁴ cm⁻³ near the PDR surface (Pilleri et al. 2012), up to 10⁵ or 10⁶ cm⁻³ for local clumps and filaments (Köhler et al. 2014). Chemical modeling of NGC 7023 by Murga et al. (2022) showed that including the dehydrogenation balance changes the PAH evolution depending on the assumed age of the nebula: a major fraction of the PAHs survive after 10⁴ years, while most are destroyed by 10⁵ years. Detailed photochemical modeling (Montillaud et al. 2013; Andrews et al. 2015; Murga et al. 2022), consistent with the constraints provided by our data, will be necessary to understand the interplay between the processes and how the environmental conditions lead to the observed outcome.

Fig. 13 Profile diagnostics diagram based on the 6.2 and 7.7 µm profiles, comparing our results for NGC 7023 (colorful circles) with those for the five Orion Bar template spectra (labeled markers, Van De Putte et al. 2025). The color bar indicates the position in the NGC 7023 mosaic and the black marker bar is the location of DF1. In the atomic PDR regime, the NGC 7023 spectra exhibit much more strongly pronounced 7.7 µm narrow features compared to the Orion Bar “Atomic” template. The 6.2 µm FWHM tends to a constant value, and the upper panel shows a map of the 6.2 µm FWHM with DF1 marked by the dotted line. For the more shielded regime (NGC 7023 beyond DF1; Orion Bar HII, DF2, and DF3), the data points follow a similar path in the diagram for both objects.

5.5 Outlook: The path to fullerenes

The central region of NGC 7023 is known to exhibit fullerene (C₆₀) emission within a distance of ∼25″ from HD 200775 (Sellgren et al. 2010; Berné & Tielens 2012). This detection with Spitzer was based on the main features of C₆₀ at 7.05, 7.5, 17.4, and 18.9 µm (Cami et al. 2010). None of these characteristic features are observed in the spectra of this GTO program, confirming that no C₆₀ is formed this close to the DF. While we do observe a 17.4 µm feature that gets stronger toward the star, this feature can also be explained by PAHs that emit at this wavelength (Sellgren et al. 2010; Shannon et al. 2016). In this work we emphasize the difference between the spatial distributions of the 17.4 µm and 16.4 µm features, and how the 17.4 µm band becomes gradually more prominent while the profiles of several other bands evolve alongside it (5.7, 6.2, 7.7, 11.3, 12.7 µm). The band ratio gradients shown in Fig. 11 suggest that the changes in the band profiles will continue beyond the edge of our spectroscopic field of view. This suggests that observing regions progressively closer to the star would trace further changes in the carbonaceous emission bands alongside the appearance of the C₆₀ emission, which would allow a characterization of the photoprocessing stages that eventually lead to the formation of fullerenes. Below, we briefly summarize a previously suggested fullerene formation pathway, where PAHs are processed into fullerenes (Berné & Tielens 2012; Andrews et al. 2015). Then we provide an outlook with some suggestions as to which spectral details could contain information about the processing stages taking place.

On the irradiated side of DF1, the “grandPAH” scenario may be taking place (Andrews et al. 2015), in which PAHs are created when the UV radiation field evaporates them from other existing materials such as VSGs or PAH clusters (Rapacioli et al. 2005; Pilleri et al. 2012). The initial PAH species are processed into a population of “grandPAHs”, a smaller set of species consisting of more photostable PAHs which are larger, more compact, and more symmetric. This scenario was proposed to explain how the carbonaceous band spectrum appears mostly insensitive to the radiation field parameter (G₀) driving different PDRs (Andrews et al. 2015).

The fullerene formation scenario of Berné & Tielens (2012) involves further photoprocessing of the initial population of larger PAHs. For the conditions in the area where C₆₀ is observed, the PAH photoprocessing models by Montillaud et al. (2013) suggest that PAHs with a typical size of 50–70 carbon atoms can be fully dehydrogenated through photodissociation of the H atoms, leaving a population of similarly sized carbon clusters. The next photoprocessing steps for the formation of fullerenes then involve the removal of carbon atoms from the hexagonal grid, which introduces five-membered rings that allow the fragments to curve (Berné & Tielens 2012; Berné et al. 2015). When such fragments form closed cages, lab experiments confirm that they can shrink via the removal of C-C pairs, until a stable fullerene is formed (Zhen et al. 2014).

An alternative fullerene formation scenario involves the processing of hydrogenated amorphous carbon grains (HAC), observed as certain spectral plateaus at 6–9 and 10–13 µm in planetary nebulae that contain fullerenes (Bernard-Salas et al. 2012). Photoprocessing the larger HAC grains creates structures of a mixed aromatic and aliphatic nature, and five-membered rings are formed by dehydrogenating some of those structures (García-Hernández et al. 2010; Micelotta et al. 2012). Such species have been named “proto-fullerenes” and may emit spectral lines that are normally associated with C₆₀. The 16.4 µm band is also considered an indicator of five-membered rings (Duley & Hu 2012).

While the 16.4 and 17.4 µm bands are both associated with skeletal C-C-C modes of PAHs, their spatial distributions are different (Sect. 5.1). One possible interpretation was via their correlation with typical PAH charge tracers (Fig. 12). One could expect that there must be certain structural differences between the neutral and charged population, so that the charged one dominates the 17.4 µm emission. Nonetheless, Boersma et al. (2010) reported that the PAH charge is not expected to strongly affect the 16–18 µm spectra. Their search within their spectral database for species with strong bands near 16.4 or 17.4 µm did not reveal a structural class that preferentially emits at 17.4 µm. Combining recent expansions of the database (Bauschlicher et al. 2018; Mattioda et al. 2020) with the depth and spatial resolution of JWST at 16–18 µm may lead to different conclusions.

Additional information about the processing may be available from the AIBs in the 10–15 µm range associated with C-H out-of-plane bending modes of PAHs. These bands are assigned to vibrational modes involving the different hydrogen adjacency classes (solo, duo, trio, quartet, quintet), which trace the edge structure of PAHs (Hudgins & Allamandola 1999; Hony et al. 2001). The ratios of these modes can be modified by changes in the edge structures of the carbon skeleton or by dehydrogenation, potentially a first step in the fullerene formation process (Berné & Tielens 2012). Fleming et al. (2010) produced maps of the same region using the PAHFIT IDL routine, and interpreted the 12.7/11.3 ratio and the associated hydrogen adjacency as an indicator of the dehydrogenation state. They found that this ratio is higher both at the bright filament and close to the star, and reaches a minimum about 10″ away from the DF in the star direction. However, other work shows that even though 11.3 µm band is dominated by solo modes, and 12.7 µm by duo or trio modes (e.g., Khan et al. 2025), the 12.7/11.3 variation in NGC 7023 is dominated by charge effects (Rosenberg et al. 2011; Peeters et al. 2012; Boersma et al. 2015; Shannon et al. 2016). The blue side of the 11.3 µm band seems to be related to charge tracers (Sect. 5.3), but other interpretations are possible, such as the temperature or mass distribution (Candian & Sarre 2015), or the presence of molecules containing oxygen or magnesium (Sadjadi et al. 2015b).

Changes in the PAH size are also a piece of the puzzle, as the carbon fragments in the scenario by Berné & Tielens (2012) need to be larger than 60 atoms. A variety of size tracers will be explored in upcoming work (E. Roscoe et al., in preparation). Previous PAH size investigations in NGC 7023 suggest that small PAHs are broken down in the cavity and only large species survive (Boersma et al. 2013, 2014; Croiset et al. 2016; Boersma et al. 2026). Our map of the 17.4 µm band would support this, as previous work based on PAHdb found that compact PAHs of 50–130 carbon atoms consistently show a feature near this wavelength (Boersma et al. 2010). In Sect. 5.4 an initial look into the size was provided through the 6.2 µm FWHM map of Fig. 13, where the narrower 6.2 µm band was interpreted as larger PAH sizes. It appears that the 6.2 µm FWHM tends to a constant minimum in the direction of the star, though this may also indicate that this tracer becomes insensitive to the PAH size at that point.

A JWST program that completed during Cycle 4 (ID: GO-8007; PI: Tielens) has observed three regions in the cavity of NGC 7023 with MIRI MRS, with a primary goal of characterizing the C₆₀ cations. These spectra will also reveal how much further the AIB profiles keep evolving, further probing the processing sequence of the AIB carriers between the NGC 7023 NW filament and the star.

6 Summary and conclusions

We applied the PAHFIT spectral decomposition tool to extract the brightness of individual subcomponents of the carbonaceous emission bands observed with the NIRSpec and MIRI IFUs in NGC 7023. The spatially resolved maps of the individual features presented in this work reveal that the emission band subcomponents display distinct spatial distributions. We organized the features into types I, II, and III, according to their spatial distributions, where the roman numeral is an indicator of the brightness in the atomic PDR region (ATM) relative to the main dissociation front (DF1), in ascending order (Sect. 4.3.1, Table 3).

The 17.4 µm band in particular is of type III, which distinguishes it from the type I or II nature of most other emission bands when only the sum of all their subcomponents is considered. However, mapping the individual subcomponents of the 5.7, 7.7, 11.3, and 12.7 µm bands reveals that their sub-components are characterized by different spatial distributions. Consistently, blue-side subcomponents such as 11.3A are of type III, even when the redder subcomponents or the total emission of each band is of type I or II. The spatial distributions indicate the existence of at least two particular populations of carbonaceous emission band carriers, one of which is strongly emitting in the ATM region, while the other is suppressed relative to DF1. Normalizing the bands to the total 11.3 µm emission, reveals strong correlations of the 17.4 µm band with these subcomponents, as well as with the 8.6 and 11.0 µm features, which are considered to be tracers of ionized PAHs. This indicates a connection between the 17.4 µm band and a population of ionized carbonaceous emitters.

The changes in the populations may relate to the photochemical evolution of PAHs, which is driven to a more advanced stage by the increasing radiation field as the star is approached. Maps of a set of emission band profile indicators (5.7A/5.7, 7.7 narrow/7.7 broad, 11.3A/11.3, 12.7A/12.7; Fig. 11) all exhibit a positive gradient in the direction of the star that persists up to the edge of the field of view, indicating that the profiles will keep evolving for smaller distances. Since C₆₀ emission is known to be present closer to the star, these profile changes may indicate processing steps that are part of a fullerene formation pathway. Recently completed JWST spectroscopy along this spatial direction will reveal the continuation of the evolution sequence revealed in this work.

As JWST observes more spectra of similar quality, the variations in the emission profiles and their relationship to the environmental conditions will be strongly constrained. It is therefore increasingly important that theoretical models are equipped to explain profile variations or the spatial distribution of individual subcomponents (rather than the ratios of the total band intensities). Furthermore, the ensuing interpretation, in terms of the underlying populations, should be consistent with the environment (radiation field, density, presence of atomic, or molecular hydrogen). In particular, determining which specific photochemical processes near DF1 lead to the changes in the 17.4/16.4 ratio will require models that can accurately predict the emission features in this wavelength range.

Data availability

Machine-readable versions of the maps presented in Figs. 5, 6, 7, 8, 9, and A.2 are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/710/A60.

Acknowledgements

This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with GTO program #01192. D.V.D.P., K.D.G., and A.N.-C. are partially supported by NASA grant 80NSSC21K1294. K.M. is supported by JWST–NIRCam contract no. NAS5-02015 to the University of Arizona. MB acknowledges funding from the Belgian Science Policy Office (BELSPO) through the PRODEX project “JWST/MIRI Science exploitation” (C4000142239). Part of this work was performed at the French MIRI center with the support of CNES and the ANR-labcom INCLASS between IAS and ACRI-ST, and also supported by the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES. This work made use of Astropy³: a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration 2013, 2018, 2022).

References

Appendix A Overview of all maps

For completeness, in Fig. A.1, we illustrate the spatial distribution of the dust continuum emission as derived from the PAHFIT results. The results themselves consist of a set of best-fit scaling factors for the MBB component at each temperature. The continuum maps shown in Fig. A.1 were created by taking the sum of the continuum model functions, and evaluating the resulting continuum model at the given wavelengths. Additional feature intensity maps are shown in Fig. A.2. The wavelength and FWHM of each feature are given in the science pack published by Van De Putte et al. (2025), and modified by the updates in Table 2. A machine-readable version of these maps is available via the CDS. (Sect. Data availability).

Fig. A.1 Maps of the continuum model evaluated at different wavelengths (as indicated in by the labels in µm). Each map is normalized to its maximum value.

Appendix B Fit residuals

In Fig. B.1, we show the residuals of the DF1 fit of Figs. 2 and 3. The model approximates the spectrum within 10%, with larger deviations up to 15% taking place in the wavelength ranges of 4.5-5.2 µm and 9-11 µm. Besides fitting performance, these figures also indicate the potential presence of additional profile variations or subcomponents, such as the peaks near 6.9 and 9.8 µm.

Fig. A.2 Full overview of spatial distributions extracted from PAHFIT results. The maps show normalized quantities analogous to Fig. 5: the frequency-integrated surface brightness (W m−2 sr−1), divided by the maximum value of the map.

Fig. B.1 Residuals of the DF1 fit (same as in Figs. 2 and 3), defined as (data - model)/model. The largest deviations are in the continuum near the wings of the main emission complexes.

All Tables

All Figures

Fig. 1 Left: position of HD 200775 relative to the NGC 7023 NW filament and footprints of the NIRSpec observations. The background is the F1130W map from the MIRI imager observations of the same program (Misselt et al. 2025). The dashed rectangle indicates the coverage of Spitzer-based spectroscopy and feature maps (Werner et al. 2004; Shannon et al. 2015), which includes coverage closer to the star, where C60 was observed. The labeled footprints are the regions of interest for the studies by Berné & Tielens (2012) and Mackie et al. (2015). The diamond (a 3.6″ square) indicates the typical spectral extraction size for Spitzer-based studies. Right: zoom-in of the same image, showing the two extraction apertures of Misselt et al. (2025) used for the demonstration in Figs. 2 and 3, and the size of a single 0.35″ spaxel of our spectral extractions (Sect. 2). The edge of the mosaic is around 37″ away from the star.

In the text

Fig. 2 Zoomed-in view of the emission complexes in the ATM and DF1 regions of NGC 7023 and the corresponding PAHFIT decompositions. The aperture-extracted spectra (defined in Sect. 2 and Fig. 1) are shown in the top panels with a normalization and offset that match the base and peak of a certain feature (first and second columns scaled to 3.3 µm, third column to 5.2 µm). These still have the emission lines included. The gray lines in the bottom two rows are the modified ATM and DF1 spectra (emission lines removed), used as input for PAHFIT. The components of the decomposition are shown as labeled in the legend.

In the text

	Fig. 3 Analogous to Fig. 2, for 6–13 µm and 13–23 µm. The normalization in the top panels is based on the 7.7 and 16.4 µm features, respectively. The updated tuning for the 16–18 µm complex better captures the shape changes, resulting in a more robust extraction of the power of the 16.4 and 17.4 µm features.
In the text

	Fig. 4 Detailed view of the updated decomposition in the 16–18 µm range. The light gray line shows the ATM spectrum, with bumps at 16.7 and 17.2 µm that raise the need for three components to model the central “plateau.” The individual Drude profiles (colorful curves and labels) are plotted after adding the continuum components, and the total model is the solid black line.
In the text

	Fig. 5 Power maps of the main bands based on sums of the subcomponents contributing to each band as defined in Table 1. The color scale represents the spectrally integrated power (W m−2 sr−1) normalized to the brightest spaxel. The dashed black vertical line corresponds to DF1 (more precisely: the maximum of 16.4 µm when averaged vertically) and the ATM region is in the leftmost 3″ of the image.
In the text

Fig. 6 Maps of the individual components of the 5.7 µm band decomposition. The lower panel is analogous to Fig. 5. The top panel is a legend illustrating the subcomponents labeled A, B, and C. The light gray line shows the ATM spectrum data, the solid black line is the best fitting total model, and the dashed line shows the sum of all other components (continuum and wings of nearby Drude profiles), to which the three 5.7 µm components were added.

In the text

	Fig. 7 Maps of the 7.7 µm decomposition and legend analogous to Fig. 6. The sum of the narrow components (as in Table 1) is shown in addition to the individual maps.
In the text

	Fig. 8 Analogous to Fig. 6 for the 11.0 and 11.3 µm bands.
In the text

	Fig. 9 Analogous to Fig. 6 for the 12.7 µm band. Note the similarity of 12.7A to 11.3A in Fig. 8.
In the text

Fig. 10 Spatial profiles along cuts perpendicular to the DF for selected features, normalized to the peak of each cut. The cuts are organized by spatial distribution types I, II, and III, as introduced in Sect. 4.3. The 17.4 µm cut is shown twice to directly compare it to the type II and type III cuts. The dashed lines in the bottom panel indicate the region over which the maps were median-stacked along the vertical axis to obtain the cuts, with the 17.4 µm map from Fig. 5 in the background. The shaded areas indicate the widths of the ATM and DF1 regions.

In the text

Fig. 11 Maps of ratios describing variations in the band profiles. The color bars represent a linear scale between the minimum and maximum. The ratios represent the relative contributions of subcomponents to the total as defined in Table 1. The ratios of the bluest components of each complex are shown alongside the 17.4/16.4 ratio to illustrate their similarity, and how the steepest gradient is in the same location.

In the text

Fig. 12 Correlation plots exploring the relation between the 16.4 and 17.4 µm bands and a selection of other bands. The data points represent the spaxels of the maps, and their colors represent the offset from the edge of the mosaic, with the locations of ATM and DF1 indicated on the color bar. The black lines show linear fits through the origin (y = ax, solid) and with a variable intercept (y = bx + c, dashed). These lines are typically closer for 17.4 µm, and more separated for 16.4 µm, with the exception of the 12.7 µm panels. This indicates that the 7.7 narrow, 8.6, 11.3A, and 12.7A features relate better with 17.4 µm, while the total 12.7 µm flux relates beter with 16.4 µm (Sect. 5.2).

In the text

Fig. 13 Profile diagnostics diagram based on the 6.2 and 7.7 µm profiles, comparing our results for NGC 7023 (colorful circles) with those for the five Orion Bar template spectra (labeled markers, Van De Putte et al. 2025). The color bar indicates the position in the NGC 7023 mosaic and the black marker bar is the location of DF1. In the atomic PDR regime, the NGC 7023 spectra exhibit much more strongly pronounced 7.7 µm narrow features compared to the Orion Bar “Atomic” template. The 6.2 µm FWHM tends to a constant value, and the upper panel shows a map of the 6.2 µm FWHM with DF1 marked by the dotted line. For the more shielded regime (NGC 7023 beyond DF1; Orion Bar HII, DF2, and DF3), the data points follow a similar path in the diagram for both objects.

In the text

	Fig. A.1 Maps of the continuum model evaluated at different wavelengths (as indicated in by the labels in µm). Each map is normalized to its maximum value.
In the text

	Fig. A.2 Full overview of spatial distributions extracted from PAHFIT results. The maps show normalized quantities analogous to Fig. 5: the frequency-integrated surface brightness (W m−2 sr−1), divided by the maximum value of the map.
In the text

	Fig. B.1 Residuals of the DF1 fit (same as in Figs. 2 and 3), defined as (data - model)/model. The largest deviations are in the continuum near the wings of the main emission complexes.
In the text