The PAH emission characteristics of the reflection nebula NGC2023

TL;DR

This study maps and analyzes PAH, C60, and H2 emissions in NGC2023, revealing spatial distributions, spectral features, and potential PAH evolution influenced by radiation exposure.

Contribution

It provides detailed spectral maps of PAH emissions, explores their spatial and spectral variations, and investigates PAH size, charge, and structural effects using computational models.

Findings

PAH emission bands show correlated spatial distributions.

Distinct components contribute to 7-9 micron PAH emission.

PAH features are independent of underlying plateaus.

Abstract

We present 5-20 micron spectral maps of the reflection nebula NGC2023 obtained with the Infrared Spectrograph SL and SH modes on board the Spitzer Space Telescope which reveal emission from polycyclic aromatic hydrocarbons (PAHs), C60, and H2 superposed on a dust continuum. We show that several PAH emission bands correlate with each other and exhibit distinct spatial distributions revealing a spatial sequence with distance from the illuminating star. We explore the distinct morphology of the 6.2, 7.7 and 8.6 micron PAH bands and find that at least two spatially distinct components contribute to the 7--9 micron PAH emission in NGC2023. We report that the PAH features behave independently of the underlying plateaus. We present spectra of compact oval PAHs ranging in size from C_66 to C_210, determined computationally using density functional theory, and investigate trends in the band positions and relative intensities as a function of PAH size, charge and geometry. Based on the NASA Ames PAH database, we discuss the 7--9 micron components in terms of band assignments and relative intensities. We assign the plateau emission to very small grains with possible contributions from PAH clusters and identify components in the 7--9 micron emission that likely originates in these structures. Based on the assignments and the observed spatial sequence, we discuss the photochemical evolution of the interstellar PAH family as they are more and more exposed to the radiation field of the central star in the evaporative flows associated with the PDRs in NGC2023.

Tables

Table 1. Table 1 Observation log

	north position		south position
map coordinatesa	5:41:40.65, -2:13:47.5		5:41:37.63, -2:16:42.5
mode	SL	SH	SL	SH
PID	50511		30295	20097
AORs	26337024		17977856	14033920
cycle x ramp time	2x14s	2x30s	2x14s	1x30s
pointings $\parallel$	1	7	3	12
step size $\parallel$	26′′	5′′	26′′	5.65′′
pointings $\perp$	20	15	18	12
step size $\perp$	1.85′′	2.3′′	3.6′′	4.7′′
Backgrounda,b	5:42:1.00, -2:6:54.5		5:40:26.21, -2:54:40.4

Table 2. Table 2 The number of solo and duo CH groups for the molecules studied.

molecule	solo	duo	molecule	solo	duo
C1 C24H12	0	12	C3 C96H24	12	12
O1 C32H14	2	12	O3 C112H26	14	12
A1 C40H16	4	12	A3 C128H28	16	12
T1 C48H18	6	12	T3 C144H30	18	12
C2 C54H18	4	12	C4 C150H30	18	12
O2 C66H20	8	12	O4 C170H32	20	12
A2 C78H22	10	12	A4 C190H34	22	12
T2 C90H24	12	12	T4 C210H36	24	12

Table 3. Table 3 List of emission features by groupings presented in Section  3.3 and schematically illustrated in Figure  21 , their aromatic hydrocarbon vibrational assignments, likely carriers, and dominant charge. Note that the carriers of the features are not exclusively due to species in one charge state, but are often dominated by species in a particular charge state; size ranges are quite broad but do show population maxima (Schutte et al., 1993 ) ; and that the features that comprise the groups deduced here are not confined only to the zones depicted in Figure  21 , but do reach their maxima in these zones.

Group	Feature	Vibrational Assignments	Carriers1	Predominant
				Charge
Molecular Cloud Edge
1	continuum		VSG
	10–15 $\mu m$  Plateau	Blend of aromatic CHoop bends from	Small PAHs, large PAHs with irregular edges,	0
		solo, duo, trio and quartet H atoms	PAH clusters, PAH-rich VSGs and small amorphous C particles
	8.1 $\mu m$  extreme
2	5–10 $\mu m$  Plateau	Blend of weak aromatic CC stretches and CHip bends	Small PAHs, large PAHs, PAH clusters, PAH-rich VSGs	0
	G8.2 $\mu m$  component	Aromatic CC stretches and CHip bends	Large PAHs, PAH clusters, PAH-rich VSGs	0
	G7.8 $\mu m$  component	Aromatic CC stretches and CHip bends	Large PAHs, PAH clusters, PAH-rich VSGs	0
PDR Front2
3	15–18 $\mu m$  Plateau		Large PAHs	0
	15.8 $\mu m$  band	?	Large PAHs	0
	11.2 $\mu m$  band	CHoop bend from solo H atoms	Large PAHs	0
PDR Front2 & Cavity
4	17.8 $\mu m$  band	?	Large PAHs	+
	12.7 $\mu m$  band3	Blend of CHoop bends from duo and trio H atoms	Large PAHs	+4
5	12.7 $\mu m$  band3	Blend of CHoop bends from duo and trio H atoms	Large PAHs	+4
	16.4 $\mu m$  band	?	Large PAHs	+
	7.7 $\mu m$  band	CC stretches and CHip bends	Large PAHs	Roughly 50/50 +,0
	6.2 $\mu m$  band	CC stretches	Large PAHs	+
Cavity
6	8.6 $\mu m$  band (GS)	CHip bends	Large PAHs	+
	G7.6 $\mu m$  component	Aromatic CC stretches and CHip bends	Large PAHs	+
	G8.6 $\mu m$  component	CHip bends	Large PAHs	+
7	8.6 $\mu m$  band (LS)	CHip bends	Large PAHs	+
	11.0 $\mu m$  band	CHoop bend from solo H atoms	Large PAHs	+
	17.4 $\mu m$  band	?	Large PAHs	+
8	7.35 $\mu m$  extreme
Vicinity of Exciting Star5
	17.4 $\mu m$  band	CC stretches and bends	C60	0
	19.2 $\mu m$  band	CC stretches and bends	C60	0

Table 4. Table 4 The intrinsic intensity ratios a for the ions for an excitation of 8 eV.

Name Formula	cation					anion
	$I_{8.6}$/$I_{6.2}$	$I_{7.6}$/$I_{6.2}$	$I_{8.2}$/$I_{6.2}$	$I_{7.8}$/$I_{6.2}$		$I_{8.6}$/$I_{6.2}$	$I_{7.6}$/$I_{6.2}$	$I_{8.2}$/$I_{6.2}$	$I_{7.8}$/$I_{6.2}$
C1 C24H12	0.06	0.81	0.12	0.01		0.14	0.97	0.28	0.95
C2 C54H18	0.16	0.53	0.19	0.28		0.40	0.97	0.29	0.90
C3 C96H24	0.48	1.01	0.10	0.28		1.15	1.05	0.14	1.10
C4 C150H30	1.41	0.94	0.08	1.30		2.28	0.84	0.16	2.12
O1 C36H16	0.39	0.65	0.18	0.31		0.08	2.44	0.16	0.95
O2 C66H20	0.34	1.13	0.21	0.56		0.68	2.05	0.26	2.38
O3 C112H26	0.68	1.35	0.10	1.21		1.51	1.85	0.87	4.00
O4 C170H32	1.86	0.97	1.07	3.81		2.83	1.62	2.43	8.03
A1 C40H16	0.19	1.79	0.81	0.38		0.10	2.61	0.32	1.30
A2 C78H22	0.72	2.52	0.46	1.50		0.42	1.87	0.54	4.37
A3 C128H28	0.78	1.47	0.38	4.06		0.82	1.11	0.78	5.62
A4 C190H34	1.20	0.79	0.83	4.15		2.56	1.33	1.89	5.79
T1 C48H18	0.35	2.96	1.00	0.35		0.25	2.31	0.48	1.40
T2 C90H24	1.44	3.21	0.49	2.98		0.61	1.18	0.40	2.36
T3 C144H30	1.76	1.45	0.61	3.71		1.87	1.67	0.62	3.16
T4 C210H36	1.80	0.96	0.64	3.28		3.61	2.09	1.23	4.53
Avg 1 C1+O1+A1+T1	0.21	1.28	0.40	0.21		0.14	2.02	0.30	1.13
Avg 2 C2+O2+A2+T2	0.43	1.27	0.27	0.82		0.51	1.38	0.34	2.07
Avg 3 C3+O3+A3+T3	0.75	1.23	0.21	1.54		1.29	1.33	0.49	2.84
Avg 4 C4+O4+A4+T4	1.51	0.91	0.51	2.73		2.66	1.28	1.06	4.23
Avg 5 C3 to T4	1.04	1.11	0.32	1.99		1.81	1.31	0.71	3.37
Ave 6 C2 to T4	0.77	1.18	0.30	1.47		1.23	1.34	0.54	2.79
Avg 7 Coronene Family	0.52	0.78	0.13	0.48		1.08	0.97	0.21	1.24
Avg 8 Ovalene Family	0.63	1.20	0.26	1.14		1.34	1.91	0.84	3.92
Avg 9 Anthracene Family	0.88	1.69	0.54	3.06		1.08	1.47	0.95	5.16
Avg 10 Tetracene Family	1.67	1.83	0.58	3.32		1.56	1.50	0.62	3.02
a $I_{7.6}$ is the sum of all intensity in the range 7.365 to 7.8 $\mu$m, $I_{7.8}$ in the range 7.8 to 8.08 $\mu$m, $I_{8.2}$ in the
range 8.095 to 8.365 $\mu$m, $I_{8.6}$ in the range 8.408 to 8.752 $\mu$m, and $I_{6.2}$ in the range 6.2 to 6.6 $\mu$m. A redshift
of 15 cm-1 is applied before calculating the intensity ratios. See Section 5.1.2 for details.
The observations for NGC 2023 fall in the region 1.49 2.56, 0.35 1.17,
0.13 0.40, and 0.37 0.72 (Fig. 23).

Table 5. Table A.1 Parameters of the fits to the observed SL (top) and SH (bottom) correlations and their Pearson’s correlation coefficients.

PAH intensity ratios1		weighted correlation coefficient			A2			B2			C3
y vs. x		N+S	N	S
6.2/11.2	G8.6/11.2	0.9721	0.9755	0.9749	0.2583	$\pm$	0.0160	1.5154	$\pm$	0.0147	1.7476	$\pm$	0.0045
6.2/11.2	11.0/11.2	0.9376	0.9320	0.9390	0.7682	$\pm$	0.0205	12.8764	$\pm$	0.2044	20.3914	$\pm$	0.1560
6.2/11.2	12.7/11.2	0.8884	0.8462	0.9164	-0.5619	$\pm$	0.0489	7.7560	$\pm$	0.1480	6.0840	$\pm$	0.0230
7.7 (LS)/11.2	6.2/11.2	0.9713	0.9828	0.9729	-0.4170	$\pm$	0.0357	2.0201	$\pm$	0.0191	1.8020	$\pm$	0.0042
7.7 (LS)/11.2	8.6 (LS)/11.2	0.9363	0.9346	0.9392	1.1322	$\pm$	0.0351	4.2934	$\pm$	0.0636	6.2698	$\pm$	0.0389
7.7 (LS)/11.2	11.0/11.2	0.9216	0.9328	0.9208	1.0326	$\pm$	0.0495	26.9987	$\pm$	0.4939	37.0259	$\pm$	0.2652
7.7 (LS)/11.2	12.7/11.2	0.8809	0.8218	0.9139	-1.6614	$\pm$	0.1059	15.9992	$\pm$	0.3206	11.0683	$\pm$	0.0512
G7.6/11.2	G8.6/11.2	0.9872	0.9887	0.9870	0.1574	$\pm$	0.0189	3.0693	$\pm$	0.0174	3.2105	$\pm$	0.0045
G7.6/11.2	11.0/11.2	0.9317	0.9342	0.9322	1.1601	$\pm$	0.0425	26.3498	$\pm$	0.4240	37.6225	$\pm$	0.2621
G7.6/11.2	12.7/11.2	0.9197	0.8647	0.9449	-0.4904	$\pm$	0.0251	4.9616	$\pm$	0.0759	3.5100	$\pm$	0.0143
(G7.6+G7.8)/11.2	6.2/11.2	0.9391	0.9539	0.9447	1.1117	$\pm$	0.0466	1.7208	$\pm$	0.0246	2.3025	$\pm$	0.0077
(G7.6+G7.8)/11.2	G86/11.2	0.9502	0.9450	0.9525	1.5320	$\pm$	0.0369	2.6292	$\pm$	0.0337	4.0092	$\pm$	0.0170
(G7.6+G7.8)/11.2	11.0/11.2	0.8751	0.8676	0.8800	2.2136	$\pm$	0.0578	24.3539	$\pm$	0.5796	46.3829	$\pm$	0.4653
(G7.6+G7.8)/11.2	12.7/11.2	0.8363	0.7307	0.8753	-0.1693	$\pm$	0.1147	14.3072	$\pm$	0.3477	13.8005	$\pm$	0.0536
G8.2/11.2	G7.8/11.2	0.7462	0.8220	0.7519	0.1082	$\pm$	0.0107	0.3508	$\pm$	0.0111	0.4630	$\pm$	0.0022
G8.6/11.2	11.0/11.2	0.9465	0.9507	0.9468	0.3393	$\pm$	0.0125	8.4567	$\pm$	0.1245	11.7402	$\pm$	0.0761
G8.6/11.2	12.7/11.2	0.9197	0.8647	0.9449	-0.4904	$\pm$	0.0251	4.9616	$\pm$	0.0759	3.5100	$\pm$	0.0143
8.6 (LS)/11.2	6.2/11.2	0.9370	0.9369	0.9391	-0.3666	$\pm$	0.0132	0.4735	$\pm$	0.0069	0.2883	$\pm$	0.0020
8.6 (LS)/11.2	11.0/11.2	0.9507	0.9518	0.9519	0.0055	$\pm$	0.0087	5.9965	$\pm$	0.0862	6.0484	$\pm$	0.0290
8.6 (LS)/11.2	12.7/11.2	0.9411	0.9145	0.9588	-0.5709	$\pm$	0.0145	3.4803	$\pm$	0.0439	1.8109	$\pm$	0.0135
12.7/11.2	11.0/11.2	0.9025	0.8622	0.9363	0.1676	$\pm$	0.0031	1.7087	$\pm$	0.0306	3.3512	$\pm$	0.0305
plat.7/11.2	6.2/11.2	0.8702	0.9107	0.8608	2.5193	$\pm$	0.1218	2.6167	$\pm$	0.0638	3.9329	$\pm$	0.0192
plat.7/11.2	7.7 (LS)/11.2	0.8343	0.9153	0.8210	3.1482	$\pm$	0.1156	1.2678	$\pm$	0.0335	2.1812	$\pm$	0.0129
plat.7/11.2	plat11/11.2	-0.1112	0.3343	-0.1735	-89.0386	$\pm$	41.8972	60.5135	$\pm$	26.1939	4.8298	$\pm$	0.0389
7.7 (LS)/6.2	8.6 (LS)/6.2	0.3374	-0.1060	0.3976	0.4450	$\pm$	0.1193	4.8382	$\pm$	0.4239	6.4170	$\pm$	0.0434
(G7.6+G7.8)/6.2	G86/6.2	0.4472	0.5021	0.4378	-1.3032	$\pm$	0.2145	6.2744	$\pm$	0.3716	4.0188	$\pm$	0.0093
G8.2/G8.6	G7.8/G8.6	0.8326	0.9646	0.8165	0.1398	$\pm$	0.0065	0.2944	$\pm$	0.0080	0.4646	$\pm$	0.0021
plat.7/G8.6	7.7 (LS)/G8.6	0.4106	0.6013	0.4023	-40.7359	$\pm$	3.6317	15.0708	$\pm$	1.1572	2.0776	$\pm$	0.0108
11.2/G8.6	12.7/G8.6	0.7307	0.8088	0.7544	-1.9183	$\pm$	0.0959	9.6982	$\pm$	0.3385	3.0312	$\pm$	0.0231
plat. 7/plat.11	6.2/plat.11	0.9469	0.9600	0.9470	1.0956	$\pm$	0.0493	3.0234	$\pm$	0.0456	4.0115	$\pm$	0.0207
plat. 7/plat.11	7.7 (LS)/plat.11	0.9174	0.9521	0.9153	1.2655	$\pm$	0.0566	1.5907	$\pm$	0.0289	2.2246	$\pm$	0.0137
12.7/11.2	11.0/11.2	0.8598	0.9026	0.9142	0.0969	$\pm$	0.0067	1.8465	$\pm$	0.0605	2.7079	$\pm$	0.0279
13.5/11.2	11.0/11.2	0.8129	0.9064	0.7955	0.0179	$\pm$	0.0011	0.1804	$\pm$	0.0100	0.3340	$\pm$	0.0045
13.5/11.2	12.7/11.2	0.7474	0.7725	0.7658	0.0096	$\pm$	0.0018	0.0935	$\pm$	0.0060	0.1250	$\pm$	0.0013
13.5/11.2	16.4/11.2	0.7874	0.8264	0.8071	-0.0017	$\pm$	0.0021	0.5920	$\pm$	0.0309	0.5670	$\pm$	0.0052
13.5/11.2	17.4/11.2	0.5938	0.4945	0.6156	0.0165	$\pm$	0.0015	1.3549	$\pm$	0.0957	2.4152	$\pm$	0.0534
16.4/11.2	11.0/11.2	0.9002	0.9228	0.9534	0.0315	$\pm$	0.0010	0.3191	$\pm$	0.0086	0.5915	$\pm$	0.0062
16.4/11.2	12.7/11.2	0.9133	0.8818	0.9192	0.0147	$\pm$	0.0015	0.1730	$\pm$	0.0050	0.2212	$\pm$	0.0012
17.4/11.2	11.0/11.2	0.9167	0.7917	0.9317	0.0006	$\pm$	0.0005	0.1412	$\pm$	0.0042	0.1464	$\pm$	0.0012
17.4/11.2	12.7/11.2	0.7676	0.7181	0.8211	-0.0043	$\pm$	0.0011	0.0677	$\pm$	0.0037	0.0539	$\pm$	0.0007
17.4/11.2	16.4/11.2	0.7982	0.6805	0.8398	-0.0101	$\pm$	0.0013	0.3918	$\pm$	0.0188	0.2445	$\pm$	0.0034
diffuse 17.4/11.2	11.0/11.2	0.6812	0.8953	0.5791	0.0035	$\pm$	0.0012	0.0621	$\pm$	0.0096	0.0895	$\pm$	0.0021
diffuse 17.4/11.2	12.7/11.2	0.7397	0.9237	0.7000	-0.0052	$\pm$	0.0020	0.0483	$\pm$	0.0060	0.0332	$\pm$	0.0007
diffuse 17.4/11.2	16.4/11.2	0.5825	0.8311	0.5471	-0.0019	$\pm$	0.0023	0.1817	$\pm$	0.0317	0.1560	$\pm$	0.0036
17.8/11.2	13.5/11.2	0.8143	0.0000	0.8143	0.0002	$\pm$	0.0010	0.3425	$\pm$	0.0269	0.3476	$\pm$	0.0048
17.8/11.2	16.4/11.2	0.8210	0.0000	0.8210	-0.0007	$\pm$	0.0009	0.2072	$\pm$	0.0139	0.1967	$\pm$	0.0024
19.0/11.2	11.0/11.2	0.6486	-0.1943	0.6552	-0.0009	$\pm$	0.0022	0.1647	$\pm$	0.0198	0.1569	$\pm$	0.0058
19.0/11.2	12.7/11.2	0.3716	-0.5228	0.4091	0.0009	$\pm$	0.0039	0.0523	$\pm$	0.0128	0.0554	$\pm$	0.0025
19.0/11.2	16.4/11.2	0.5476	-0.2527	0.5702	-0.0117	$\pm$	0.0040	0.4295	$\pm$	0.0594	0.2577	$\pm$	0.0108
11.2/12.7	13.5/12.7	0.3424	0.0000	0.3424	-6.9288	$\pm$	1.5988	82.9106	$\pm$	12.9525	27.3724	$\pm$	0.4232
15.8/12.7	11.2/12.7	0.8878	0.0000	0.8878	-0.0200	$\pm$	0.0039	0.0231	$\pm$	0.0012	0.0173	$\pm$	0.0003
16.4/12.7	11.2/12.7	0.6118	0.6826	0.6001	0.1714	$\pm$	0.0050	0.0151	$\pm$	0.0015	0.0659	$\pm$	0.0008
16.4/12.7	plat.11/12.7	0.7017	0.7446	0.6987	0.1703	$\pm$	0.0038	0.0095	$\pm$	0.0007	0.0401	$\pm$	0.0005
16.4/12.7	plat.17/12.7	0.6980	0.6305	0.7156	0.1445	$\pm$	0.0055	0.0727	$\pm$	0.0052	0.2079	$\pm$	0.0018
plat.11/12.7	11.2/12.7	0.8862	0.8865	0.9012	-0.7051	$\pm$	0.2186	1.8476	$\pm$	0.0703	1.6229	$\pm$	0.0114
plat.11/12.7	plat.17/12.7	0.8794	0.8744	0.8866	-3.2396	$\pm$	0.2889	8.1680	$\pm$	0.2745	5.1377	$\pm$	0.0391
plat.17/12.7	11.2/12.7	0.8785	0.9253	0.8717	0.2831	$\pm$	0.0266	0.2345	$\pm$	0.0081	0.3192	$\pm$	0.0020
12.7/16.4	11.2/16.4	-0.2435	-0.2654	-0.2428	5.1265	$\pm$	0.1473	-0.0468	$\pm$	0.0097	0.3068	$\pm$	0.0041
12.7/15-18plat.	11.2/15-18plat.	0.0481	0.0714	0.0490
16.4/sum(15-20)	12.7/11.2	0.7424	0.6412	0.7663	0.0961	$\pm$	0.0033	0.1955	$\pm$	0.0107	0.5022	$\pm$	0.0047

Table 6. Table A.2 C-H band position maximum (in  μ 𝜇 \mu m), total intensity (in km/mol), and intensity per C-H (in km/mol/H). No redshift is applied.

Molecule	Cation			Neutral			Anion
	$\lambda$     $I$   $I$(CH)			$\lambda$      $I$    $I$(CH)			$\lambda$      $I$   $I$(CH)
C1 C24H12	3.240	1.8	0.1	3.263	8.9	0.7	3.298	20.7	1.7
O1 C32H14	3.242	3.3	0.2	3.263	10.8	0.8	3.291	24.3	1.7
A1 C40H16	3.245	4.9	0.3	3.264	13.0	0.8	3.288	27.4	1.7
T1 C48H18	3.247	6.8	0.4	3.264	14.7	0.8	3.284	30.3	1.7
C2 C54H18	3.249	7.6	0.4	3.264	15.5	0.9	3.284	28.7	1.6
O2 C66H20	3.250	9.6	0.5	3.264	18.2	0.9	3.280	32.5	1.6
A2 C78H22	3.252	11.8	0.5	3.265	21.0	1.0	3.281	36.3	1.7
T2 C90H24	3.254	14.6	0.6	3.266	23.5	1.0	3.279	40.8	1.7
C3 C96H24	3.256	16.0	0.7	3.268	24.9	1.0	3.279	39.5	1.6
O3 C112H26	3.256	18.2	0.7	3.266	27.4	1.1	3.276	43.1	1.7
A3 C128H28	3.258	20.6	0.7	3.268	30.6	1.1	3.277	47.1	1.7
T3 C144H30	3.259	24.5	0.8	3.268	34.3	1.1	3.277	53.6	1.8
C4 C150H30	3.261	27.1	0.9	3.268	34.4	1.1	3.276	52.8	1.8
O4 C170H32	3.261	29.0	0.9	3.272	38.6	1.2	3.276	56.8	1.8
A4 C190H34	3.262	32.2	0.9	3.269	43.1	1.3	3.277	62.4	1.8
T4 C210H36	3.263	37.6	1.0	3.270	47.6	1.3	3.276	69.8	1.9

Full text

1. Introduction

The mid-infrared (IR) spectra of many astronomical objects are dominated by strong emission bands at 3.3, 6.2, 7.7, 8.6, 11.3 and, 12.7 ${\rm\mu m}$, generally attributed to polycyclic aromatic hydrocarbon molecules (PAHs) and related species. Despite the significant progress made in our understanding of these emission bands (further referred to as PAH bands), the specifics of their carrier remain unclear. Indeed, no single PAH molecule or related species has been firmly identified to date. This lack of detailed knowledge hampers our understanding of the PAH emission bands and their use as a diagnostic tool for probing the local physical conditions.

The PAH emission bands show clear variations in peak positions, shapes and (relative) intensities, not only between sources, but also spatially within extended sources (e.g. Hony et al., 2001; Peeters et al., 2002; Smith et al., 2007b; Galliano et al., 2008). In particular, it is well established that the 3.3 and 11.2 ${\rm\mu m}$  PAH bands correlate with each other and that the 6.2, 7.7 and 8.6 ${\rm\mu m}$  PAH bands correlate tightly with each other. Laboratory and theoretical PAH studies have long indicated that the main driver behind these correlations is the PAH charge state: emission from neutral PAHs dominates at 3.3 and 11.2 ${\rm\mu m}$  while ionized PAHs emit strongest in the 5-10 ${\rm\mu m}$  and are responsible for the 11.0 ${\rm\mu m}$  band (e.g. Hudgins et al., 2004, and references therein). However, Whelan et al. (2013) and Stock et al. (2014) recently reported that the strong correlation between the 6.2 and 7.7 ${\rm\mu m}$  PAH bands breaks down on small spatial scales towards the giant star-forming region N66 in the Large Magellanic Cloud and towards the massive Galactic star-forming region W49A. This suggests that in addition to PAH charge, other parameters such as molecular structure influence the astronomical 6.2 and 7.7 ${\rm\mu m}$  PAH emission. A dependence on PAH parameters such as charge, molecular structure, size can be probed by investigating theoretical calculations of PAH spectra (Langhoff, 1996; Bauschlicher & Langhoff, 1997; Ellinger et al., 1998; Pauzat & Ellinger, 2002). Advances in computing power and computational methods now allow spectra for larger PAHs and for a much wider range of PAH species to be calculated and thus for a better and larger scan of the parameter space (Mulas et al., 2006; Malloci et al., 2007; Bauschlicher et al., 2008, 2009; Ricca et al., 2010, 2012; Candian et al., 2014; Candian & Sarre, 2015). Hence, a systematic study of fully resolved PAH emission at high spatial scale over the full mid-IR bandwidth of a large sample of objects in combination with a systematic investigation of theoretical PAH spectra promises to reveal a more detailed and nuanced view of the characteristics of the emitting PAH population.

The high sensitivity and broad wavelength coverage of the Spitzer Space Telescope has provided a unique opportunity to study the variations of the PAH features within extended sources that have been spectrally mapped. In particular, reflection nebulae (RNe) are very interesting in this respect as they are known to be excellent laboratories to determine the characteristics of interstellar dust and gas due to their high surface brightness and relatively simple geometry. Moreover, in contrast to H ii regions, observations of RNe are not confused by diffuse emission due to recombination and cooling of ionized gas, simplifying their analysis greatly. Consequently, the well-known reflection nebulae NGC 7023 and NGC 2023 have been prime targets for PAH studies (e.g. Abergel et al., 2002; Werner et al., 2004b; Sellgren et al., 2007; Compiègne et al., 2008; Fleming et al., 2010; Peeters et al., 2012; Pilleri et al., 2012; Boersma et al., 2013, 2014; Shannon et al., 2015; Stock et al., 2016).

In this paper, we present Spitzer-IRS spectral maps of NGC 2023 in the 5-20 ${\rm\mu m}$  region and theoretical data from oval, compact PAHs in order to investigate and interpret the spatial behaviour of the PAH emission components, in particular in the 6 to 9 ${\rm\mu m}$  region. Section 2 discusses the observations, the data reduction, typical spectra and the continuum and flux determination. The data are analyzed in terms of features’ emission maps and correlations in Section 3 and theoretical data from oval, compact PAHs are presented in Section 4. These are combined and discussed in Section 5. Finally, we end with a short summary in Section 6.

2. The data

2.1. Observations

The observations were taken with the Infrared Spectrograph (IRS, Houck et al., 2004) on board the Spitzer Space Telescope (Werner et al., 2004a) and were part of the Open Time Observations PIDs 20097, 30295 and 50511.

We obtained spectral maps for two positions in the reflection nebula NGC 2023 (see Fig. 1): towards the dense shell south south-west (-11*′′, -78′′) of the exciting star HD 37903 corresponding to the H2 emission peak, and towards a region to the north (+33′′, +105′′) of the exciting star (Burton et al., 1998). The north position is characterized by a much lower density ($\sim$104* cm*-3*) than that at the south position, where densities exceed 105 cm*-3* (Burton et al., 1998; Sheffer et al., 2011). Likewise, the radiation field between 6 and 13.6 eV in the northern region ($\sim$500 G0111G0 is the integrated 6 to 13.6 eV radiation flux in units of the Habing field = 1.6x10*-3* ergs/cm2/s.) is lower than the UV field where the H2 peaks in the southern region ($\sim$104 G0, Burton et al., 1998; Sheffer et al., 2011).

The spectral maps are made with the short-low and short-high modes (respectively SL and SH). The SL mode covers a wavelength range of 5-15 ${\rm\mu m}$  at a resolution ranging from 60 to 128 in two orders (SL1 and SL2) and has a pixel size of 1.8*′′, a slit width of 3.6′′and a slit length of 57′′(resulting in 2 pixels per slit width). The SH mode covers a wavelength range of 10-20 ${\rm\mu m}$  at a resolution of $\sim$ 600 and has a pixel size of 2.3′′, a slit width of 4.7′′and a slit length of 11.3′′*(resulting in 2 pixels per slit width). We obtained background observations in SH for the north position and in SL for the south position. Table 1 gives a detailed overview of the observations. The 15-20 ${\rm\mu m}$  SH spectra have been presented in Peeters et al. (2012), hereafter paper I.

2.2. Reduction

The SL raw data were processed with the S18.7 pipeline version by the Spitzer Science Center. The resulting bcd-products are further processed using cubism (Smith et al., 2007a) available from the SSC website222http://ssc.spitzer.caltech.edu.

As discussed in detail in paper I, the background adds only a small contribution to the on-source PAH flux except for source D from Sellgren (1983). Consequently, we did not apply a background subtraction and exclude source D from the PAH analysis for the remainder of this paper.

We applied cubism’s automatic bad pixel generation with $\sigma_{TRIM}=7$ and Minbad-fraction = 0.50 and 0.75 for the global bad pixels and record bad pixels. We excluded the spurious data at the extremities of the SL slit by applying a wavsamp of 0.05 to 0.95 for the north map and of 0.06 to 0.94 for the south map. Remaining bad pixels were subsequently removed manually.

Spectra were extracted from the spectral maps by moving, in one pixel steps, a spectral aperture of 2x2 pixels in both directions of the maps. This results in overlapping extraction apertures. Slight mismatches in flux level were seen between the different orders of the SL module (SL1 and SL2). These were corrected by scaling the SL2 data to the SL1 data. The applied scaling factor averaged to 5 and 7% for the south and north SL map respectively. Around 8% of the spectra in the south SL map have scaling factors of more than 20%, which are all located at high x-values (i.e. at the south side of the map, below the S ridge). Hence, the applied scaling does not influence our results.

We’ve noticed an apparent small (of the order of a few $\times$ 0.01 ${\rm\mu m}$) wavelength shift in some 2x2 SL spectra, in particular for the south map (see Appendix A). While its effects can be noticed in the top row of the north map and in the bottom two rows of the south map in most feature intensity maps presented in this paper (Figs. 5, 6, 9, 10, A.4, A.5, and A.11), it does not change the conclusions of this paper.

For the reduction of the SH data, we refer to paper I for a detailed discussion. Spectra were extracted from the spectral maps by moving, in one pixel steps, a spectral aperture of 2x2 pixels in both directions of the maps; the same procedure as for the SL data.

2.3. The spectra

Figs. 2 and 3 show typical spectra towards NGC 2023. The complete 5-20 ${\rm\mu m}$  spectra reveal a weakly rising dust continuum, H2 emission lines, C60 emission, and a plethora of PAH emission bands. In addition to the main PAH bands at 6.2, 7.7, 8.6 and 11.2 ${\rm\mu m}$, weaker bands are detected at 5.7, 6.0, 11.0, 12.0, 12.7, 13.5, 14.2, 15.8, 16.4, 17.4 and, 17.8 ${\rm\mu m}$. These PAH bands are perched on top of broad emission plateaus from roughly 5–10, 10–15 and, 15–18 ${\rm\mu m}$. C60 exhibit bands at 7.0, 8.6, 17.4 and 19 ${\rm\mu m}$  (Cami et al., 2010; Sellgren et al., 2010). The latter two are clearly present in these spectra (paper I).

NGC 2023 contains a cluster of young stars (Sellgren, 1983). YSO source D is located in the southern quadrant of the south SL map and just outside the south SH map (see Fig. 1). The spectra of regions close to source D show PAH emission features as well as characteristics typical of YSOs, i.e., a strong dust continuum and a (strong) CO2 ice feature near 15 ${\rm\mu m}$. Furthermore, the background contribution to the 11.2 PAH flux is significant for a large fraction of the spectra containing ice-features. As in paper I, we therefore excluded these spectra in the analysis. The 15-20 ${\rm\mu m}$  spectrum of YSO source C is found in both the SL and SH south map (see Fig. 1). This YSO has the same spectral characteristics as the spectra across NGC 2023 but with enhanced surface brightness. However, its SL2 spectrum suffers from instrumental effects (i.e. a sawtooth pattern likely due to the undersampled IRS PSF). Hence, we excluded this source in the analysis of the SL data but included it for the analysis of the SH data.

2.4. Continuum subtraction and band fluxes

The extinction towards NGC 2023 is estimated to be small: $A_{K}=0.2-0.65$ mag (for a detailed overview, see Sheffer et al., 2011). Moreover, little variation in the extinction is found across the south map except around source C and D (Pilleri et al., 2012; Stock et al., 2016). An extinction $A_{K}=0.4$ mag results in extinction corrections of 11 % for the 8.6 and 11.2 ${\rm\mu m}$  PAH bands (because of its overlap with the silicate band) and between 4-7 % for the other features considered here. In this paper, we determined the PAH fluxes assuming zero extinction.

We applied three different decomposition methods to the SL data. For the first method, we subtract a local spline (LS) continuum from the spectra, consistent with the method of Hony et al. (2001) and Peeters et al. (2002) as shown in Fig. 2. This continuum is determined by using anchor points at roughly 5.4, 5.5, 5.8, 6.6, 7.0, 8.2, 9.0, 9.3, 9.9, 10.2, 10.5, 10.7, 11.7, 12.1, 13.1, 13.9, 14.7, and 15.0 ${\rm\mu m}$. The fluxes of the main PAH bands are then estimated by integrating the continuum subtracted spectra while those of the weaker bands and the H2 lines are measured by fitting a Gaussian profile to the band/line. The central wavelength of the Gaussian used to measure the H2 line fluxes was allowed to vary, to correct for the wavelength shifts. The 6.0, 11.0, and 12.7 ${\rm\mu m}$  features need special treatment because of blending. Since the 6.0 and 6.2 ${\rm\mu m}$  PAH features are blended and given the low resolution of the SL data, we extracted the 6.0 band intensity by fitting two Gaussians with $\lambda$ (FWHM) of 6.026 (0.099) and 6.229 (0.1612) ${\rm\mu m}$  respectively to the data (excluding the red wing of the 6.2 PAH band). These values were obtained by taking the average over all spectra when fitted by two Gaussians having peak positions and FWHM that were not fixed. As can be seen in Fig. 4, this decomposition works reasonably well except, of course, for the red wing of the 6.2 ${\rm\mu m}$  band. The 6.0 ${\rm\mu m}$  PAH flux is then subtracted from the integrated flux of the 6.0+6.2 ${\rm\mu m}$  PAH bands to obtain the 6.2 ${\rm\mu m}$  PAH flux. A similar decomposition method can be applied for the 11.0 and 11.2 ${\rm\mu m}$  bands by fitting two Gaussians with $\lambda$ (FWHM) of 10.99 (0.154) and 11.26 (0.236) ${\rm\mu m}$  respectively (Fig. 4, middle panel). Analogously to the 6 ${\rm\mu m}$  region, this decomposition works remarkably well except for the red wing of the 11.2 ${\rm\mu m}$  band. In order to obtain the fluxes of the 12.7 ${\rm\mu m}$  PAH band and the 12.3 ${\rm\mu m}$  H2 line, we fit a second order polynomial to the local ‘continuum’ (i.e. the blue wing of the 12.7 ${\rm\mu m}$  PAH band) and a Gaussian profile to the H2 line. The flux of the 12.7 ${\rm\mu m}$  PAH band was then obtained by subtracting the H2 flux from that obtained by the (combined) integrated flux of the 12.7 ${\rm\mu m}$  PAH band and the H2 line. We estimated the signal-to-noise ratio of the features as follows: $S/N=F~{}/~{}(rms\times\sqrt{N}\times\Delta\lambda$) with $F$ the feature’s flux [W/m2/sr], $rms$ the rms noise, $N$ the number of flux measurements within the feature and $\Delta\lambda$ the wavelength bin size as determined by the spectral resolution. Although the rms noise is a measure of how accurate the continuum can be determined, this method does not take into account the error in the continuum measurement.

The second method is a modification of the first method. A global spline (GS) continuum is determined by using the same anchor points as in the first method except for the continuum point at roughly 8.2 ${\rm\mu m}$  (Fig. 2). This affects the band profiles and intensities of the 7.7 and 8.6 ${\rm\mu m}$  PAH bands and the underlying plateau. This also creates a new plateau underneath only the 7.7 and 8.6 ${\rm\mu m}$  PAH bands that is then defined by the difference of the local spline continuum and the global continuum and is further referred to as the 8 ${\rm\mu m}$  bump. The plateau continuum (Fig. 2) is obtained by using continuum points at roughly 5.5, 9.9, 10.2, 10.4, and 14.7 ${\rm\mu m}$. The underlying plateaus between 5-10 and 10-15 ${\rm\mu m}$  are then defined by the difference of this continuum with the GS and LS continua respectively. The fluxes are obtained in the same way as discussed above for the LS continuum.

Finally, the third approach employs PAHFIT to analyze the data (Fig. 2, Smith et al., 2007b). The PAHFIT decomposition results in components representing the dust continuum emission (which is a combination of modified blackbodies), the H2 emission and the PAH emission. In particular, the PAH emission is fit by a combination of Drude profiles. A detailed discussion of this decomposition method and its results can be find in Appendix C.

For the SH data, we applied the local spline continuum by using anchor points at roughly 10.2, 10.4, 10.8, 11.8, 12.2, 13.0, 13.9, 14.9, 15.2, 15.5, 16.1, 16.7, 16.9, 17.16/17.19 (for the south and north maps, respectively), 17.6, 18.14/18.17 (for the north and south maps, respectively), 18.5, 19.3, and 19.4 ${\rm\mu m}$  and the plateau continuum by using anchor points at 10.4, 15.2, 18.5 and 19.4 ${\rm\mu m}$. The fluxes of the 11.0, 11.2, 12.7 ${\rm\mu m}$  PAH bands and the 12.3 ${\rm\mu m}$  H2 line are determined in the same way as for the SL data. The parameters for the 11.0 and 11.2 Gaussians are $\lambda$ (FWHM) of 11.018 (0.1205) and 11.243 (0.144) ${\rm\mu m}$  respectively (Fig. 4, right panel). The fluxes of the weaker 12.0 and 13.5 ${\rm\mu m}$  bands are determined by fitting a Gaussian. The 14.2 ${\rm\mu m}$  PAH displays a weaker blue shoulder. These two components were fitted by a Gaussian ($\lambda$ (FWHM) of 13.99 (0.1178) and 14.230 (0.1830) ${\rm\mu m}$  respectively).

The SL and SH fluxes for features in the 10-14 ${\rm\mu m}$  region tend to differ. This difference arises from the following reasons. Firstly, we do not regrid the SL data to the coarser SH grid and do not correct for the different spatial PSFs. However, when the SL data are regridded to the SH grid and the SL data are scaled to the SH data but no correction for the different spatial PSFs are made, it is clear that the features’ strength differ to varying degrees in the SL and SH data (see Fig. 3). Secondly, the continuum is less accurately determined in the SL data due to blending of the emission features (11.2 PAH, 12.0 PAH, H2 and 12.7 PAH) which strongly influences the fluxes of the weaker bands. Hence, we will refrain from comparing the SL and SH fluxes directly with each other and we will not use the weaker 12.0, 13.5 and 14.5 ${\rm\mu m}$  bands in the SL data.

The applied continuum is clearly not unique. Hence neither is the decomposition of the PAH emission, nor the calculated band strengths. Thus, for comparison, we did a full analysis of the SL data using these three methods. In the remainder of the paper, the default applied continuum is the local spline continuum for the intensities of the PAH bands (excluding the 5–10 ${\rm\mu m}$  plateau333Defined as the difference between the plateau continuum and the global spline continuum.) unless stated otherwise (e.g. Section 3.1.2). We will discuss the influence of the continuum and decomposition on the results where necessary.

3. Data analysis

Here we investigate the relationship between individual PAH emission bands, the underlying plateaus, H2 emission, C60 emission, and the dust continuum emission present in the 5-20 ${\rm\mu m}$  region. The north map is characterized by lower flux levels and therefore larger scatter is present in the maps/plots tracing weaker PAH features. Hence, when discussing the spatial distribution of the 14.2, 15.8 and 17.8 ${\rm\mu m}$  bands, we restrict ourselves to the south map.

3.1. SL data

The spatial distribution of the various emission components in the 5-15 ${\rm\mu m}$  SL data are shown in Figs. 5 and 6 (the range in colours is set by the minimum and maximum intensities present in the map, and a local spline continuum is applied except for the 5–10 ${\rm\mu m}$  plateau) and feature correlations in Fig. 7 (with a local spline continua applied except for the 5–10 ${\rm\mu m}$  plateau and the Gaussian components as discussed in Section 3.1.2). To exclude the influence of PAH abundance and column density in the correlation plots, we normalized the band fluxes to that of another PAH band.

3.1.1 Overall appearance

The following trends are derived from the south spectral maps shown in Fig. 5. The 11.2 ${\rm\mu m}$  feature, the 5-10 and 10-15 ${\rm\mu m}$  plateaus and the continuum emission show very similar spatial morphology with distinct peaks at the S and SSE positions. The 8 ${\rm\mu m}$  bump exhibits a similar morphology but shows more enhanced emission at the SE ridge and west of the S’ ridge. The PDR front is well traced by the H2 emission, which also clearly peaks at the S and SSE positions and is heavily concentrated along these two ridges only. In contrast, the distribution of the 6.2, 7.7, 12.7 ${\rm\mu m}$  features are displaced towards the illuminating star and away from the S ridge, with the loss of emission at the S position accompanied by a rise at the S’ and SE ridges. Put another way, the 6.2, 7.7 and 12.7 ${\rm\mu m}$  features show very similar spatial morphology with distinct peaks at the SSE and SE ridges and with weaker peaks at the S and S’ positions. In addition, they show broad, diffuse emission NW of the line connecting the S and SSE ridges. The 8.6 ${\rm\mu m}$  band is further displaced towards the illuminating star: it peaks at the SE and S’ position but does not peak at the S and SSE ridges as do the 6.2, 7.7 and 12.7 ${\rm\mu m}$  emission. In fact, it lacks emission in the S ridge. Similar to the 6.2, 7.7 and 12.7 ${\rm\mu m}$  emission, it does show a broad, diffuse plateau NW of the line connecting the S and SSE ridges. As does the 8.6 ${\rm\mu m}$  PAH emission, the 11.0 ${\rm\mu m}$  PAH emission also peaks at the SE and S’ ridges and lacks emission in the S ridge. But it is also very strong at the broad, diffuse plateau N and NW of the line connecting the S and SSE ridges. The 6.0 ${\rm\mu m}$  PAH emission peaks at the S and SSE ridges but its peaks seem to be displaced towards the east compared to the 11.2 ${\rm\mu m}$  PAH emission. It has weaker emission maxima in the form of an arc south of the S’ ridge and west of the SE ridge. Hence, 6.0 ${\rm\mu m}$  PAH emission seems to be somewhat unique.

The variety in the spatial distribution of the different emission components is even more pronounced in the north map (Fig. 6). The 11.2 ${\rm\mu m}$  PAH emission peaks in the NW ridge. In contrast with the south map, the continuum flux and the 10-15 ${\rm\mu m}$  plateau are displaced from the 11.2 ${\rm\mu m}$  PAH emission and peaks in the N ridge while the 5-10 ${\rm\mu m}$  plateau and the 8 ${\rm\mu m}$  bump trace both the N and (part of) the NW ridge. The 6.2 and 7.7 ${\rm\mu m}$  PAH emission are again very similar and differ from the 11.2 ${\rm\mu m}$  PAH emission. While they also peak in the NW ridge as does the 11.2 ${\rm\mu m}$  PAH feature, they show decreased emission in the northern part of this ridge and extend towards the west in the southern part of this ridge (note the difference in the black and pink contours which trace the 11.2 and 7.7 ${\rm\mu m}$  PAH bands respectively). The morphology of the 12.7 ${\rm\mu m}$  PAH emission is a bit of both that of the 11.2 and 7.7 ${\rm\mu m}$  PAH emission. The 8.6 ${\rm\mu m}$  emission is distinct from the 6.2 and 7.7 ${\rm\mu m}$  PAH emission and peaks slightly west of the southern part of the NW ridge. This is also seen in the spatial distribution of the 11.0 ${\rm\mu m}$  PAH which peaks even further west of the southern part of the NW ridge compared to the 8.6 ${\rm\mu m}$  PAH emission. The 6.0 ${\rm\mu m}$  PAH emission in the north map, as for the south map, has a unique morphology which is highly concentrated and peaks at the intersection of the N and NW ridge. The spatial distribution of the H2 line intensities varies as well: the S(3) [9.7 ${\rm\mu m}$] intensity peaks between the maxima in the 11.2 and 7.7 ${\rm\mu m}$  PAH intensities along NW ridge and the S(2) at 12.3 ${\rm\mu m}$ intensity peaks in N ridge. These S(2) and S(3) distributions are confirmed by the maps obtained with PAHFIT for both the north and south FOVs (see Appendix C).

Fig. 7 shows observed intensity correlations; their fit parameters and correlation coefficients can be found in table D and their line cuts in Fig. A.7. The well known, very tight correlation between the 6.2 and 7.7 ${\rm\mu m}$  PAH bands is also observed within our sample, consistent with their similar spatial morphology. Note that the observed correlation is close to a 1:1 correlation (i.e. through (0,0)). Surprisingly, the 8.6 ${\rm\mu m}$  band correlates with the 6.2 and 7.7 ${\rm\mu m}$  bands despite their differing spatial distributions. However, this correlation is not as tight as that between the 6.2 and 7.7 ${\rm\mu m}$  bands, but overall, in keeping with their close but differing spatial distributions. The enhanced scatter in the correlation of the 8.6 ${\rm\mu m}$  band with the 6.2 and 7.7 ${\rm\mu m}$  bands has regularly been attributed to the influence of extinction and/or the larger uncertainty in determining the 8.6 ${\rm\mu m}$  band intensity. However, we found that larger deviations from the line of best fit are found in locations where the spatial distribution of the bands is different, i.e., the observed scatter originates in the distinct spatial distributions. The 11.0 ${\rm\mu m}$  PAH band correlates best with the 8.6 ${\rm\mu m}$  PAH band although not as tight as the 6.2 with the 7.7 ${\rm\mu m}$  PAH bands. The 11.0 vs. 8.6 correlation exhibits exactly a 1:1 dependence (i.e. the best fit goes through (0,0); see Fig. 7) and has a high correlation coefficient. This may not be immediately clear from their spatial distribution for the south map as shown in Fig. 5 which is attributed to the bottom two rows (in y-direction; see Appendix  B for details). The 11.0 ${\rm\mu m}$  PAH band also correlates with the 6.2 and 7.7 ${\rm\mu m}$  PAH bands but with slightly more scatter. Here as well, the 11.0 ${\rm\mu m}$  PAH emission has a different spatial distribution than the 8.6, 6.2 and 7.7 ${\rm\mu m}$  PAH emission resulting in enhanced scatter in these correlation plots and thus a lower correlation coefficient (ranging from 0.949 to 0.958 compared to 0.978 for the 6.2 vs. 7.7 correlation). The 12.7 ${\rm\mu m}$  PAH band also shows a dependence on the 6.2, 7.7, 8.6 and 11.0 ${\rm\mu m}$  PAH bands but this connection is clearly weaker than those amongst the 6.2, 7.7, 8.6 and 11.0 ${\rm\mu m}$  bands (with correlation coefficients ranging between 0.930 to 0.945).

3.1.2 Decomposition of the 7 to 9 ${\rm\mu m}$  region

The distinct spatial distribution of the 7.7 and 8.6 ${\rm\mu m}$  PAH bands prompts further investigation. As discussed in Sect. 2.4, the chosen local spline continuum clearly influences band intensities. However, if an emission feature is due to a single carrier or distinct subset of the PAH population, the spatial distributions of its sub-components should all be identical, independent of how these sub-components have been defined. Hence, the distinct spatial distribution of the 7.7 and 8.6 ${\rm\mu m}$  PAH emission indicates that they originate in multiple carriers or loosely related PAH subpopulations. In an attempt to resolve the sub-components of the 7 to 9 ${\rm\mu m}$  PAH emission, each related to a single carrier or subset of the PAH population, we subtracted the global spline continuum (GS, see Fig. 2, magenta line) from the spectra and decomposed the remaining PAH emission into four Gaussians with $\lambda$ (FWHM) of 7.59 (0.450), 7.93 (0.300), 8.25 (0.270), and 8.58 (0.344) ${\rm\mu m}$  respectively (see Fig. 8)444Note that the G8.2 component is not present when using the local continuum (LS, see Fig. 2, orange line) as for this continuum, an anchor point at 8.2 ${\rm\mu m}$  is taken and this component thus becomes part of the 8 ${\rm\mu m}$  bump. When using the plateau continuum (see Fig. 2, green line), the G8.2 component sits on top of the plateau emission.. These values were obtained by taking the average over all spectra when fitted by four Gaussians having peak positions and FWHM that were not fixed but were constrained to fit these four bands. These components are further referred to as G7.6, G7.8, G8.2 and G8.6 bands. We chose four components to represent the 8.6 ${\rm\mu m}$  PAH band and the 7.6 and 7.8 ${\rm\mu m}$  subcomponents of the 7.7 ${\rm\mu m}$  complex, whose ratio, 7.6/7.8, exhibit spatial variation with distance from the star (Bregman & Temi, 2005; Rapacioli et al., 2005; Boersma et al., 2014). A fourth component is needed to obtain a good fit in the 7 – 9 ${\rm\mu m}$  region. We did not include a Gaussian at 7.4 ${\rm\mu m}$  because none of the spectra in the map show a ’feature’ near 7.4 ${\rm\mu m}$  like that observed in the ISO-SWS data of NGC 7023 (Moutou et al., 1999). The ‘nominal’ 7.7 ${\rm\mu m}$  PAH complex is dominantly comprised of the G7.6 component and only originates by a relatively small fraction in the G7.8 component (see Appendix E). The G7.6 and G8.6 ${\rm\mu m}$  components exhibit almost identical spatial distributions indicating there is no significant contribution from an additional, spatially distinct component in either of the G7.6 and G8.6 ${\rm\mu m}$  components (Fig.  9). Both components peak at the SE, SSE and S’ ridges in the south map. In the northern FOV, they peak at the southern and centre portion of the NW ridge, with an extension towards the west in the southern part of this ridge. The G7.8 and G8.2 ${\rm\mu m}$  components also have a similar spatial distribution, though discrepancies are present (Fig.  9). They peak at the S and SSE ridges in the south map, similar to the 11.2 ${\rm\mu m}$  emission and the 5-10 and 10-15 ${\rm\mu m}$  plateau emission. In the north map, they trace both the N, centre and northern portion of the NW ridge, which is similar to the 10-15 ${\rm\mu m}$  plateau and the 10.2 ${\rm\mu m}$  continuum emission. As reported for the ratio of the 7.6 and 7.8 ${\rm\mu m}$  subcomponents (Bregman & Temi, 2005; Rapacioli et al., 2005; Boersma et al., 2014), the ratio of G7.6/G7.8 traces the different environments very well: it is clearly smallest where H2 emission peaks and inside the molecular cloud (see Appendix E).

These results are also found in the correlation plots (Fig. 7). The G7.6 and G8.6 ${\rm\mu m}$  components exhibit a very tight correlation (best correlation coefficient of 0.987) which closely resembles a 1:1 relation (i.e. it goes through (0,0)). In contrast, the correlation of the G7.8 and G8.2 ${\rm\mu m}$  Gaussian components shows more scatter resulting in a correlation coefficient of 0.833. Remarkably, this enhanced scatter largely originates in the south map. Indeed, the correlation coefficient for the south map is 0.817 while that for the north map is a whopping 0.965555The same holds when normalized to the 6.2 ${\rm\mu m}$  band with correlation coefficients of 0.823, 0.787 and 0.961 for respectively both the north and south maps. This scatter originates in the slight mismatch of their spatial distributions and either indicates the shortcomings of our decomposition or is due to the fact that they arise in different PAH sub-populations.

Despite the arbitrariness of this decomposition (by assuming a decomposition into four Gaussians), we can conclude that at least 2 spatially distinct components contribute to the PAH emission in the 7 to 9 ${\rm\mu m}$  region. In this respect, it is very enlightening to watch an animation showing the change in spatial distribution, as a function of wavelength, of the continuum-subtracted emission (applying the global spline continuum, Fig. 10; animation is available online). The spatial distribution of the PAH emission in the 7 to 9 ${\rm\mu m}$  region continuously varies between two extremes which are found at around $\sim$7.35 and $\sim$8.1 ${\rm\mu m}$  in both the north and south map666Note that the morphology of the two extremes depends on the chosen continuum (here, the global spline continuum). Using instead the local spline continuum or the plateau continuum will remove the contribution captured by the 8 ${\rm\mu m}$  bump in the former case and add the contribution of the 5-10 ${\rm\mu m}$  plateau in the latter case. These components have a different morphology as that of the 7.35 ${\rm\mu m}$  extreme and are more similar to that of the 8.1 ${\rm\mu m}$  extreme. . For the south map, the extreme $\sim$7.35 ${\rm\mu m}$  emission is spatially very similar to that of the 11.0 ${\rm\mu m}$  PAH emission (Fig. 5): it traces the S’, SE and SSE ridges, the horizontal filaments in the northern region of the map and the broad, diffuse plateau NW of the line connecting the S and SSE ridges. The extreme $\sim$8.1 ${\rm\mu m}$  emission is similar to the G8.2 component (Fig. 9) and peaks at the S and SSE ridges. So does the 11.2 ${\rm\mu m}$  PAH emission (Fig. 5) but the extreme $\sim$8.1 ${\rm\mu m}$  emission has enhanced emission in the horizontal filaments in the northern region of the map and the broad, diffuse plateau compared to the 11.2 ${\rm\mu m}$  PAH emission. For the north map, the extreme distribution at $\sim$7.35 ${\rm\mu m}$  peaks at the southern part of the NW ridge and the extension towards the west of the NW ridge like the 8.6 and 11.0 ${\rm\mu m}$  PAH emission (Fig. 6) and the G8.6 component (Fig. 9). In contrast, the extreme $\sim$8.1 ${\rm\mu m}$  emission peaks in the N ridge as does the 10.2 ${\rm\mu m}$  continuum (Fig. 6) and the G7.8 component (Fig. 9).

3.1.3 Implications for other PAH bands

The morphology of the PAH emission changes continuously with wavelength for all major PAH bands (i.e. the 6.2, 11.2 and, 12.7 ${\rm\mu m}$  bands), as for the south map, albeit to a considerable lesser extent (Fig. 10): none of the wavelengths exhibit a morphology as extreme as that of the $\sim$7.35 ${\rm\mu m}$  or the $\sim$8.1 ${\rm\mu m}$  extremes.

At $\sim$6.0 ${\rm\mu m}$, the morphology is similar to that of 6.0 PAH shown in Figs. 5 and 6 for both maps. For the south map, with increasing wavelength from 6.0 to 6.2 ${\rm\mu m}$, emission increases in the SE ridge, followed by increased diffuse emission NW of the line connecting the SSE and S ridge, and increased emission in the S’ ridge. This is accompanied by decreased emission in the S ridge resulting in a spatial distribution very similar to that of the integrated 6.2 PAH emission map. Subsequently, for longer wavelengths ($>\,\sim$6.25 ${\rm\mu m}$), the diffuse emission and the emission in the S’ ridge decreases and the S ridge becomes a tiny bit stronger while the peak emission is found in the SE and SSE ridge. Considering the north map, at $\sim$6.1 ${\rm\mu m}$, the emission peaks in the NW ridge. With increasing wavelength up to 6.2 ${\rm\mu m}$, emission in the southern part and west of the southern part of the NW ridge increases while the northern part of the NW ridges fades a little, resulting in a morphology somewhat between the 6.2 PAH and the G8.6 emission. Subsequently, the NW ridge becomes a bit stronger and eventually, the emission west of the southern part of the NW ridge fades away.

The most dramatic change in the 11 ${\rm\mu m}$  region is found from 11.07 to 11.14 ${\rm\mu m}$  with the switch from the 11.0 to 11.2 PAH emission. The morphologies near 11.0 and 11.2 ${\rm\mu m}$, as in the 6 ${\rm\mu m}$  region, are well represented by that of the 11.0 PAH emission and 11.2 PAH emission respectively. Once beyond the peak of the 11.2 PAH band, the morphology continues to vary, albeit to a less dramatic extent, and becomes more sharply peaked. In the south map, the emission peaks sharply at the S and SSE ridges (and weaker at S’ ridge), very similarly to the $H_{2}$ emission. In contrast, for the north map, the emission peaks in the northern part of the NW ridge while the remaining of the NW ridge has similar intensity as the N ridge.

The morphology of the integrated 12.7 PAH band closely resembles that at 12.7 ${\rm\mu m}$. Going to shorter wavelengths, emission is enhanced at locations where the ions peak (as traced by the 8.6 and 11.0 ${\rm\mu m}$  emission). This is most noticeable in the south map: emission enhances in the SE ridge, the S’ ridge (in particular the western part) and NW of the line connection the SSE and S ridge (earlier referred to as the diffuse emission). In addition, the local emission peak in the S ridge becomes displaced towards the north to eventually being merged with the diffuse emission NW of the line connection the SSE and S ridge. In the north map, this is seen by a slight enhancement of the emission in the extension to the W of the southern part of the NW ridge (but the peak of the emission remains in the NW ridge).

3.1.4 PAH features versus plateaus

In the north map, the PAH emission and dust continuum emission are displaced from each other (i.e. they peak at different ridges), giving a unique opportunity to further explore the relationship between the PAH emission features and the underlying plateaus. Indeed, in the north map the 10–15 ${\rm\mu m}$  plateau resembles the dust continuum emission (at 10 and 14.7 ${\rm\mu m}$) but not the PAH emission (of any of the PAH bands), the 5–10 ${\rm\mu m}$  plateau emission exhibit a spatial distribution between that of the continuum and that of the PAH emission and the 15-18 ${\rm\mu m}$  plateau emission follows the PAH emission (Sections  3.1.1 and  3.2; paper I). In paper I (Fig. 8), we reported that the 15-18 ${\rm\mu m}$  PAH bands vary independently of the 15-18 ${\rm\mu m}$  plateau. Figures 11 and 12 show a similar exercise for the 5–10 and 10–15 ${\rm\mu m}$  features and their underlying plateaus respectively: we compared the PAH emission in each pixel with the average PAH emission when normalizing the strength of the plateaus. While small variations are present in the shape of the plateaus amongst the different pixels777Specifically, the largest variation in the 5–10 ${\rm\mu m}$  plateau is found between the 6.2 and 7.7 ${\rm\mu m}$  features (Fig. 11, top panel). In case of the 10–15 ${\rm\mu m}$  plateau, the most noticeable change is the ratio of the plateau strength underneath the 11.2 and 12.0 ${\rm\mu m}$  bands versus the plateau strength underneath the 13.5 and 14.2 ${\rm\mu m}$  bands (Fig. 12, top panel)., it is clear in both cases that the features behave independently of their underlying plateau and of each other. After investigating the possible influence of the applied decomposition on this result (see Appendix F), we conclude and confirm earlier reports (Bregman et al., 1989; Roche et al., 1989, paper I) that the plateaus are distinct from the features.

3.2. SH data

The spatial distribution of the various emission components in the 10-20 ${\rm\mu m}$  SH data are shown in Figs. 13 and 14 (the range in colours is set by the minimum and maximum intensities present in the map). Given that paper I gives a detailed discussion of the 15-20 ${\rm\mu m}$  emission, we will only briefly summarize these here (see also Figs. 3 and 4 in paper I). The SH data of the south map only covers the S and SSE ridges and the diffuse emission NW of the line connecting both ridges (see Figs. 1 and 5). The SH data of the north map covers both the N and NW ridge though the northern part of the NW ridge is missing (see Figs. 1 and 6). The feature correlation plots are shown in Fig. 15 and Figs. 5 and 6 in paper I. Their fit parameters and correlation coefficients can be found in table D and their line cuts in Fig. A.7. We normalized the band fluxes to that of another PAH band to exclude the influence of PAH abundance and column density in the correlation plots, as in the SL data analysis.

For the south maps a comparison with the SL maps reveals identical maps (in SL and SH) for the 11.2 PAH feature and the H2 emission bands. They clearly peak at the S and SSE positions with the H2 emission, tracing the PDR front, being more concentrated along these two ridges only. The continuum emission differs in the SL and SH data in terms of the concentration on the S and SSE ridge: the SL 14.7 ${\rm\mu m}$  emission is more extended than the SH 15 ${\rm\mu m}$  emission. The 11.0 PAH feature differs slightly in the SL and SH data: the SH data exhibit peak emission in the SSE ridge while the SL data doesn’t. This may be (partially) due to wavelength shifts occurring in the edges of the SL map (see Section 2.2) and due to different spectral resolution of the SH and SL mode. The morphology of the 12.7 PAH feature in the SH data shows a decreased emission in the S ridge and enhanced diffuse emission NW of the line connecting the S and SSE ridge compared to the 12.7 map based on the SL data, and thus is more similar to e.g. the 7.7 PAH morphology. This discrepancy between the SH and SL data likely originates in the less accurate continuum determination in the SL data due to blending of the emission features (11.2 PAH, 12.0 PAH, H2 and 12.7 PAH). The spatial distribution of the 12.0 ${\rm\mu m}$  PAH emission and the 10-15 ${\rm\mu m}$  and 15-18 ${\rm\mu m}$  plateau emission is very similar to that of the 11.2 ${\rm\mu m}$  PAH emission - although the 12.0 ${\rm\mu m}$  band shows weaker peak emission in the S ridge. The 13.5 ${\rm\mu m}$  PAH emission shows a morphology in between that of the 11.2 and 16.4 ${\rm\mu m}$  PAH emission: compared to the 11.2 ${\rm\mu m}$  PAH map, it exhibits a decreased emission in the S ridge and enhanced emission south of the S and SSE ridge and in the north of the map. In addition, the very weak 14.2 ${\rm\mu m}$  PAH band seems to have a spatial distribution similar to the 11.0 ${\rm\mu m}$  PAH feature. As discussed in paper I, the 16.4 ${\rm\mu m}$  PAH emission exhibits a similar spatial variation as does the 12.7 ${\rm\mu m}$  PAH emission. The weaker 17.8 ${\rm\mu m}$  PAH emission exhibits a morphology between that of the 11.2 PAH and 12.7 PAH emission: its emission peaks in the SSE ridge but it lacks emission in the S ridge and shows slightly enhanced diffuse emission. The 15.8 ${\rm\mu m}$  PAH band also shows similarities with the 11.2 ${\rm\mu m}$  PAH morphology but is more sharply peaked on the S and SSE ridges. Finally, the 17.4 ${\rm\mu m}$  emission (due to both PAH and C60) resembles the 11.0 emission most closely. It shows enhanced emission in the north corner of the map but this is mostly due to C60 as discussed in paper I. The diffuse emission is dominated by PAH emission and is identical to that seen for the 11.0 emission except that it exhibits decreased emission in the SSE ridge (hence overall, it is identical to the SL 11.0 emission).

For the north map, the continuum emission, the 11.0 and 11.2 ${\rm\mu m}$  PAH bands are the same in both SH and SL data. The SH 12.7 ${\rm\mu m}$  emission north map shows small deviations from the SL map, as for the south map. Specifically, for the SH map, the south section of the NW ridge seems to be brighter while the middle section seems to be fainter. In contrast with the south map, the spatial distribution of the 12.0 and 13.5 ${\rm\mu m}$  emission deviates from that of the 11.2 ${\rm\mu m}$  emission. Specifically, both bands exhibit a lack of emission in the N ridge while the emission in the NW ridge is more centred towards the lower half of the ridge for the 12.0 ${\rm\mu m}$  band and is restricted to the southern part of the ridge only for the 13.5 ${\rm\mu m}$  band. The latter is very similar to the 11.0 emission. As discussed in paper I, the 16.4 and 12.7 ${\rm\mu m}$  PAH bands exhibit similar morphologies. The 17.4 ${\rm\mu m}$  emission peaks in the NW ridge, however, due to its weakness, it’s hard to further distinguish its morphology but it seems like it encompasses the entire NW ridge. Finally, the morphology of H2 S(2) emission is the same in both the SH and SL data, though the H2 S(2) SH map is clearly of superb quality with respect to the SL map. Combining the H2 results of both SH and SL observations of the north FOV, we note that the H2 emission peaks in the N ridge for the S(2) and S(1) transitions (as does the continuum emission) while it peaks in the NW ridge for the S(3) and S(5)888 as measured using PAHFIT, see Appendix C. transitions (as does the PAH emission). In contrast, the morphology of all H2 lines in the south map is similar. A full analysis of this interesting result is beyond the scope of this paper, however the difference might be attributed to different excitation processes for low J versus high J lines. Sheffer et al. (2011) found that in the south map collisions dominate the excitation to the J = 7 level, however in the north the line maps can look quite different if the low J lines are collisionally excited while the high J lines are UV pumped. The higher J lines, whether they are UV pumped or collisionally excited, should arise from regions that are closer to the surface of the PDR where the UV field and temperatures are higher (see e.g. Sheffer et al., 2011).

The correlation plots reveal less tight correlations amongst the SH bands and hence the obtained correlation coefficients are lower compared to those obtained with the SL data. However, the correlation coefficient for the 11.0/11.2 vs. 12.7/11.2 data is similar for both the SL and SH data (see Table D). Note however that an offset is present between the SH and SL data which is likely due to the less accurate continuum determination in SL resulting in a large underestimation of the weaker band intensities. In addition, an offset is present between the north and south map in both the SL and SH data of which the origin remains unclear. If much weaker emission in the N map were the culprit and thus the accuracy of the 11.0 ${\rm\mu m}$  flux, it would effect all correlations involving the 11.0 ${\rm\mu m}$  PAH which is not the case. Nevertheless, the fact that we see similar correlation coefficients in the SL and SH data for the 11.0/11.2 vs. 12.7/11.2 data indicates that the PAH features at the longer wavelengths are intrinsically less connected with each other. Keeping this in mind, we can further investigate the behaviour of the longer wavelength bands.

In paper I, we found correlations between the 11.2 PAH, 15.8 PAH and the 15-18 ${\rm\mu m}$  plateau (correlation coefficients of 0.88–0.89; group 1) and correlations between the 16.4 and 12.7 ${\rm\mu m}$  PAH emission (correlation coefficient of 0.91, group 2). In addition to these, we find that the 11.0 ${\rm\mu m}$  PAH emission correlates well with the 12.7 and 16.4 PAH emission (correlation coefficients of 0.90–0.95; amongst highest in the SH dataset) and that the 10-15${\rm\mu m}$  plateau correlates well with the features in group 1. Also, as revealed in the south maps, the 17.4 ${\rm\mu m}$  emission correlates well with the 11.0 (correlation coefficient of 0.93). Note that this comprises not only PAH emission but also a small contribution of C60 emission, in particular in the north corner of the map (see paper I). To a lower degree, the 17.4 ${\rm\mu m}$  emission also correlates with the 12.7 and 16.4 ${\rm\mu m}$  emission having correlation coefficients of 0.82 and 0.84. Finally, the weak 13.5 ${\rm\mu m}$  emission band exhibits weaker correlations with the 11.0 and 16.4 ${\rm\mu m}$  PAH emission while the weak 12.0 ${\rm\mu m}$  PAH emission correlates well with the 11.2 ${\rm\mu m}$  PAH emission.

3.3. Spatial sequence

The spatial distributions of the individual PAH band intensities (and continuum emission) show great diversity (Sections 3.1 and 3.2, and Figs.  5, 6, 13 and 14). They reveal that the (peak) emission for these bands occurs at different distances from the illuminating source. To further exemplify this, Fig. 16 presents normalized emission profiles for PAH features, H2 lines and the continuum component for both maps. The following spatial sequence emerges (from furthest to the source to nearest away from the source):

• •

Group 1: 8.1 ${\rm\mu m}$  extreme, continuum emission, 10-15 ${\rm\mu m}$  plateau

• •

Group 2: G7.8, G8.2, 5-10 ${\rm\mu m}$  plateau

• •

Group 3: 11.2, 15.8, 15-18 ${\rm\mu m}$  plateau, 12.0

• •

Group 4: 12.7 (SL)999A slightly different morphology is present in the SL and SH data (see Section 3.2)., 17.8

• •

Group 5: 6.2, 7.7, 16.4, 12.7 (SH)15

• •

Group 6: G7.6, G8.6, 8.6 (GS)101010For clarity, we make a difference between the 8.6 ${\rm\mu m}$  strength when applying a global spline continuum (GS) or local spline continuum (LS).

• •

Group 7: 11.0, 17.4 PAH, 8.6 (LS)16

• •

Group 8: 7.35 ${\rm\mu m}$  extreme

The distinction between groups 1, 2, and 3 is made based on the N map where these features peak in either the N or NW ridge or both. And while these features exhibit a slightly different line profile in the S map, they all peak in the S ridge, as does the H2 emission (see Fig. 16). Hence, the differentiation between groups 1, 2, and 3 is not clearly present in the south map. Note that the different lines of H2 emission in the N map show different morphologies (see Section 3.2). The spatial behaviour of the 15-20 ${\rm\mu m}$  bands has already been reported in paper I and Shannon et al. (2015) with respect to the 11.0, 11.2 and 12.7 ${\rm\mu m}$  features. Specifically, groups 1 and 2 presented in paper I are part of group 3 and 5 respectively in this spatial sequence. A similar spatial behaviour is found in NGC 7023 for the 15-20 ${\rm\mu m}$  bands (Shannon et al., 2015).

The spatial morphology of individual features, and thus the spatial sequence as given above, is influenced by the applied decomposition of the spectra into individual components. Therefore, we also analyzed the spectral maps with PAHFIT and a detailed comparison with the spline method is given in Appendix C. The PAHFIT decomposition also reveals distinct spatial morphology for different components but the (detailed) spatial sequence is different. This originates in the fact that some individual features within PAHFIT trace different emission components than in the spline method: for instance, several PAHFIT features include a fraction of the underlying plateaus (as defined by the spline method). As these plateaus are independently of the features at similar wavelengths, this influences the obtained spatial sequence.

When considering a specific decomposition method however, rather than being exact, this spatial sequence is indicative of the overall, gradual and continuous (rather than discrete) change in PAH population with distance from the star. Indeed, features in consecutive bins may switch bins and/or consecutive bins may be merged. For example, some features (e.g. G8.2 and continuum emission) exhibit the same spatial distribution in the south map while being distinct in the north map. This reflects the local geometries and different conditions of the maps (the FUV radiation field is $\sim$500 G0 and $\sim$104 G0 for respectively the N and S position and the density is $\sim$104 cm*-3* and $>$105 cm*-3* respectively; Burton et al., 1998; Sheffer et al., 2011). However, the fact that they exhibit different spatial distributions in the north map suggests they are intrinsically different. Hence, towards many other objects, such a detailed sequence is likely not spatially resolved.

4. The infrared emission features and PAHs

While these IR emission features at 3.3, 6.2, 7.7, 8.6, 11.3 and 12.7 ${\rm\mu m}$  are generally assigned to a family of vibrationally excited PAHs, the identification of the specific molecules has yet to be made. Nevertheless, comparison of the observations with experimental and theoretical studies have put strong constraints on the emitting PAH population.

First, experimental and theoretical studies of PAHs clearly demonstrate the importance of PAH charge (Szczepanski & Vala, 1993; Hudgins et al., 1994; Langhoff, 1996; Kim et al., 2001; Hudgins et al., 2004, and ref. therein). The observed correlations between the main PAH bands (i.e. between the 3.3 and 11.2 ${\rm\mu m}$  bands and between the 6.2, 7.7 and 8.6 ${\rm\mu m}$  bands) are then interpreted to reflect the charge dominance of the bands: the CH stretching mode at 3.3 ${\rm\mu m}$  and the solo CH out-of-plane bending mode at 11.2 ${\rm\mu m}$  are dominated by neutral PAHs while the CC stretching modes at 6.2 and 7.7 ${\rm\mu m}$  and the CH in-plane bending mode at 8.6 ${\rm\mu m}$  are governed by ionized PAHs (Joblin et al., 1996; Hony et al., 2001; Galliano et al., 2008). At longer wavelengths, the solo CH out-of-plane bending mode at 11.0 ${\rm\mu m}$  is assigned to ionized PAHs (Hudgins & Allamandola, 1999; Hony et al., 2001) while the 15.8 ${\rm\mu m}$  band and the broad 15-18 ${\rm\mu m}$  plateau (similar to the broad 17 ${\rm\mu m}$  band as defined by PAHFIT, Smith et al., 2007b) are assigned to neutrals (Peeters et al., 2012; Shannon et al., 2015).

A secondary effect influencing the intrinsic PAH spectra is the PAH structure. The 16.4 ${\rm\mu m}$  band seems to be systematically present in PAHs exhibiting pendent rings (Moutou et al., 2000; Van Kerckhoven et al., 2000; Peeters et al., 2004; Boersma et al., 2010) and large PAHs with pointy edges (Ricca et al., 2012). However, it is also found to correlate with the 6.2 and 12.7 ${\rm\mu m}$  bands (Boersma et al., 2010; Peeters et al., 2012; Shannon et al., 2015). Hence, both charge and molecular structure likely influences the observed 16.4 ${\rm\mu m}$  band (for a detailed discussion, see Peeters et al., 2012). Likewise, the 12.7 ${\rm\mu m}$  band is assigned to duo and trio CH out-of-plane bending vibrations of large PAHs and hence depends on the PAH edge structure (Hony et al., 2001; Bauschlicher et al., 2008, 2009). It may arise from both neutral and ionized PAHs despite its correlation with the 6.2 and 7.7 ${\rm\mu m}$  bands (for a detailed discussion see Hony et al., 2001; Boersma et al., 2010; Peeters et al., 2012). In addition, recent theoretical studies have shown that PAH structure also influences the emission at shorter wavelengths, in particular for the 8.6 ${\rm\mu m}$  band. Indeed, symmetric PAHs containing approximately 100 carbon atoms111111Smaller PAHs (less than about 96 carbon atoms) do not show a clear 8.6 ${\rm\mu m}$  band (Bauschlicher et al., 2008, 2009) while larger PAHs (more than about 150 carbon atoms up to 384 carbon atoms) still had band positions consistent with observations but the relative intensities of the 6.2, 7.7 and 8.6 $\mu$m bands are not consistent with observations (Ricca et al., 2012). were required to obtain a distinct, sizeable band at 8.6 $\mu$m (Bauschlicher et al., 2008, 2009). Lowering the symmetry by adding more rings or adding irregular edge structures dramatically reduced the size and distinctiveness of the 8.6 $\mu$m band (Bauschlicher et al., 2009).

The reduction in the 8.6 $\mu$m band with lower PAH molecular symmetry and irregular edge structures suggests that the majority of the emitting PAHs have smooth edge structure. Since circular compact PAHs are inherently more thermodynamically stable than non-compact, irregular PAHs, it is not surprising that they dominate the astronomical PAH population. However, comparison of the computed and observed relative intensities of the 6.2, 7.7, and 8.6 $\mu$m bands suggests the emitting population does not contain only circular, compact species (Ricca et al., 2012). Oval, compact PAHs are also exceptionally stable (Ricca et al., 2012) and, until now, have been underrepresented in the collection of computed spectra available to compare with observations. Here, we address this deficiency and present the spectra of compact oval PAHs ranging in size from C66 to C210. The discussion on their astronomical implications will be focussed on the 6 to 9 ${\rm\mu m}$  PAH bands and applied to the presented analysis of NGC 2023 (Sect. 5).

4.1. Model and Methods

The species studied in this work are shown in Figure 17. The bottom line in the figure shows the hexagonal cores around which the rings are systematically added. They are the carbon skeletons of the molecules benzene, naphthalene, anthracene, and tetracene. We do not consider these molecules in this work as they have been discussed earlier (Langhoff, 1996), but start our study with the second row. The structures were fully optimized and the harmonic frequencies computed using density functional theory (DFT). We use the hybrid B3LYP (Becke, 1993; Stephens et al., 1994) functional in conjunction with the 4-31G basis set (Frisch et al., 1984). All of the DFT calculations are performed using Gaussian 09 (Frisch et al., 2009). The interactive molecular graphics tool MOLEKEL (Flükiger et al., 2000) is used to aid the analysis of the vibrational modes.

Previous work (Bauschlicher & Langhoff, 1997) has shown that the computed B3LYP/4-31G harmonic frequencies scaled by a single scale factor of 0.958 are in excellent agreement with the matrix isolation mid-IR fundamental frequencies of the PAH molecules.

The computed intensities obtained using 14 different functionals, including both pure and hybrids, are in good agreement with each other and with experiment for neutral naphthalene (Bauschlicher & Ricca, 2010; Szczepanski & Vala, 1993). As previously noted, the band positions obtained using these functionals also agree well with experiment. While the same 14 functionals have band positions in good agreement with experiment for naphthalene cation, their intensities are not in particularly good agreement with the two experimental results (Szczepanski et al., 1992; Hudgins et al., 1994; Bauschlicher & Ricca, 2010). The computed intensities of the ions obtained using the hybrid functionals are in good agreement, as are those obtained using the pure functionals. However, the pure functional results are somewhat smaller than those obtained using the hybrid functionals. While there is a systematic difference for the absolute intensities, the relative intensities are more similar for the hybrids and pure functionals. An inspection of the results in the PAH database shows that the good agreement between theory and experiment for band intensities of naphthalene is quite common for neutrals, as is the qualitative agreement between theory and experiment for the intensities of the cations. Given the agreement between the different functionals, the agreement of the theory with the experimental data currently available and, the limited number of PAHs that can be studied experimentally, we assume that the theory forms a more consistent set of data than the experiment, and we use the theoretical results to analyze the observed intensities (see Sect. 5.1.2). However, we need to stress that while the computed band positions have been consistently shown to be accurate when compared to available experimental data, the reliability of the computed intensities of ionized PAHs is not yet known. We suspect that changes in the relative band intensities with shape and size are more reliable than absolute intensities, so we should correctly identify trends in intensity ratios. Finally we note that any uncertainty in the intensities does not detract from our identification of molecular characteristics that lead to shifts in the band positions (see Sect. 5.1.1).

A linewidth of 30 cm*-1* is taken for the bands shortward of 9 $\mu$m, a linewidth of 10 cm*-1* is used for the bands between 10 and 15 $\mu$m and a linewidth of 5 cm*-1* is used between 15-20 $\mu$m; values consistent with current observational and theoretical understanding (see discussion in Ricca et al., 2012). For the 9 to 10 $\mu$m region, the FWHM is scaled in a linear fashion (in wavenumber space) from 30 to 10 cm*-1*. In addition to ignoring any further variations of linewidth as a function of mode, Fermi resonances as well as overtone and combination bands are not taken into account. Despite these limitations, these idealized spectra have proven very useful in better understanding the astronomical spectra.

The computational studies yield integrated band intensities in km/mol, which we broaden in wavenumber space because it is linear in energy. Thus the units of our synthetic spectra are cm*-1* for the $x$ axis and km/(mol cm*-1*) for the $y$ axis. The latter is converted to units of a cross section (given in 10${}^{5}~{}$cm2 mol*-1*). The $x$ axis is converted to $\mu$m to compare with observational results, which are commonly reported in $\mu$m.

Astronomical PAHs are typically observed as the emission from highly vibrationally excited molecules. Hence when comparing with observations, our computed 0 K absorption spectra should be shifted to the red to account for the difference between absorption and emission from vibrationally excited molecules. The size of this shift depends on many factors such as the size of the molecule, the anharmonicity of modes, and temperature of the emitting species. In the past we have redshifted the computed spectra by 15 cm*-1* to compare with observational spectra (see the discussion in Bauschlicher et al., 2009). For the presentation of the theoretical data of oval compact PAHs (Sect. 4.2), we do not apply any shift to the synthetic spectra. However, when comparing to the astronomical observations (Sect. 5.1), we include a 15 cm-1 redshift.

4.2. IR spectroscopy of the oval, compact PAHs

The presented theoretical data are part of our ongoing investigations of IR properties of large PAHs (Bauschlicher et al., 2008, 2009; Ricca et al., 2012). In particular, Ricca et al. (2012) reported a detailed description of the spectroscopic properties of the coronene family. Here, we extend this to the ovaline, anthracene, tethracene family.

The number of vibrational bands in these large molecules is so large that it is not practical to publish them in this manuscript and only plots of the spectra will be given (except for the 3 ${\rm\mu m}$  region). The full list of the band positions and intensities are available on-line in the NASA Ames PAH IR database (Bauschlicher et al., 2010; Boersma et al., 2014, www.astrochem.org /pahdb).

4.2.1 The CH stretching Vibrations (3.2 - 3.3 $\mu$m)

The C-H stretching region is very similar for all compact families (coronene, ovalene, anthracene and tethracene), with the aspect ratio having very little effect on this region of the spectra. This thus extends our conclusions for the compact circular PAHs (Ricca et al., 2012) to the compact oval PAHs. A detailed description can be found in Appendix G.

4.2.2 The CC/CH in-plane stretching Vibrations (6 - 9 $\mu$m)

The 6-9 $\mu$m region of the neutral, cation, and anion spectra are shown in Figure 18. As expected, the intensities of the neutral species are at least 10 times smaller than those of the ions and we do not consider them in detail. Inspection of the cation spectra show that 6.3, 7.7, and 8.6 $\mu$m bands increase in intensity with increasing size in each family. However, there are similarities and differences between the families. The increase in the 7.7 $\mu$m band with size is similar for the different families, while the 6.3 and 8.6 $\mu$m bands in the coronene family grow more rapidly with size than in the others. The 6.3 and 8.6 $\mu$m band intensities seem to follow a decreasing trend in family from coronene to tetracene, ovalene and finally anthracene. The 6.3 and 8.6 $\mu$m bands appear to be somewhat coupled, while the 7.7 $\mu$m band seems to vary more independently. For the anions, the coronene, ovalene, and anthracene families follow similar trends as the cation, while the tetracene family is clearly different for the anions than the cations, with the intensities of the tetracene anion being the largest of the families considered.

4.2.3 The CH out-of-plane vibrations (9 - 15 $\mu$m)

The 9-15 $\mu$m region of the neutral, cation, and anion spectra are shown in Figure 19. Excluding the smallest member of each family, where the number of duo hydrogens exceeds the number of solo hydrogens, as expected, the neutral and cation spectra look very similar as a function of size and family; the most notable change is the increase in the intensity of the solo band with size and across the families because of the increasing number of solo hydrogens. The anions show some differences with family, namely the coalescence of two bands near 11 $\mu$m occurs more slowly for the coronene family than the other families, probably because there are fewer solo hydrogens. In addition, the tetracene family shows the growth of a band near 9 $\mu$m with increasing size. It is interesting to note that while the tetracene anion seems to have more intensity in the 6-9 $\mu$m region, it seems more consistent with the other families in the 9-15 $\mu$m range.

4.2.4 The transition to skeletal vibrations (15 - 20 $\mu$m)

The 15-20 $\mu$m region of the neutral, cation, and anion spectra are shown in Figure 20. There tends to be an increase in the number of bands with increasing size, both within a family and across families. As in our previous studies (e.g. Boersma et al., 2010; Ricca et al., 2012), some bands appear to line up with some of the observed bands, however one is unable to make any assignments.

5. Discussion

5.1. Theory versus observations

The spectroscopic properties of the oval compact PAHs presented in Section 4.2 are now combined with earlier studies and analyses of the NASA Ames PAH database (PAHdb, Bauschlicher et al., 2010; Boersma et al., 2014), and applied to the PAH characteristics of NGC 2023 (as discussed in Section 3). In particular, in Sections 5.1.1 and 5.1.2, we will focus on the new decomposition of the 7 to 9 ${\rm\mu m}$  PAH emission into the four Gaussian components (G7.6, G7.8, G8.2 and G8.6) in terms of band assignments and relative intensities. We conclude with a discussion on the CH modes in Section 5.1.3.

5.1.1 Band assignments

To summarize the ensuing discussion, we ascribe the G7.6 ${\rm\mu m}$  emission primarily to compact, positively charged PAHs with sizes in the range of 50 to 100 C-atoms, the G7.8 ${\rm\mu m}$  emission to predominantly neutral, very large PAHs (100 $\leq\#$ C $<$ 150) / PAH clusters with bay regions or modified duo CH groups (like adding or removing hydrogens or substituting a N for a CH group), the G8.2 ${\rm\mu m}$  emission to PAHs/PAH clusters with multiple bay regions (e.g. PAHs with very irregular structures or corners), and the G8.6 ${\rm\mu m}$  emission to very large compact, symmetric PAHs (96 $\leq\#$C $<$ 150). The interested reader can find a detailed discussion of these assignments in the following paragraphs and of the possible cluster assignment in Section 5.2.3. Within these assignments and following the discussion in Andrews et al. (2015), the spatial morphology of the different features and their spatial sequence (as detailed in Section 3.3) then reveals the photochemical evolution of the interstellar PAH family as they are increasingly exposed to the strong radiation field of the central star in the evaporative flows associated with the PDRs in NGC 2023. As depicted in Figure 21 and summarized in Table 3, with decreasing distance from the source of radiation, the PAH family thus evolves evolves from VSG and PAH clusters comprised of PAHs with irregular edge structures and multiple bay regions to very large PAHs with multiple bay regions, irregular edges and modified duo CH groups (cf., G7.8 & G8.2 components), which are then broken down to a mixture of neutral and positively charged, compact, highly symmetric neutral PAHs (e.g., the circumcoronene family; e.g., 11.2 ${\rm\mu m}$  band). Subsequently, these too are photochemically broken down to smaller structures and in this process “pass” again through the stage of PAHs with irregular edge structures and/or corners (i.e., PAHs with duo/trio CH groups and strong 12.7 ${\rm\mu m}$  band; see Section 5.1.3). Eventually, only the most stable, compact, and highly symmetric, large ($\geq$ 70 C-atoms) PAH cations and C60 fullerenes survive in the region closest to the star (e.g., G7.6 and G8.6 components, and 11.0 ${\rm\mu m}$  band). For reference, we note that the grandPAHs as described in Andrews et al. (2015) are located in the PDR front, where the IR continuum and H2 emission peaks.

The G7.6 and G7.8 ${\rm\mu m}$  components

All PAHs with more than $\sim$ 20 carbon atoms have C-C stretching modes in the 7.6-7.8 ${\rm\mu m}$  range. The frequency of the most intense bands varies with charge, size and structure. Indeed, Bauschlicher et al. (2008, 2009) reported that PAH anions emit at slightly longer wavelengths than do the PAH cations. From their comparison of observations to theory, these authors furthermore conclude that the 7.7 ${\rm\mu m}$  complex is comprised of a mixture of small and large PAH cations and anions with the “small” species contributing more to the 7.6 ${\rm\mu m}$  component and the large species (#C $\geq$ 100) more to the 7.8 ${\rm\mu m}$  component. In addition, an upper limit to the PAH size of 150 carbon atoms is put forward by Ricca et al. (2012) based on the fact that very large compact PAHs (#C $\geq$ 150) exhibit broad complex emission between 6 to 7 ${\rm\mu m}$, which is unlike the astronomical observations. Similarly, Bauschlicher et al. (2009) reported that the 7.7 ${\rm\mu m}$  complex in irregular PAHs is broader than that of the compact PAHs and is merging with the 8.6 ${\rm\mu m}$  band (see their Fig. 8). This suggests that also the molecular edge structure determines the frequency of the most intense bands. Indeed, a further detailed investigation of the PAHdb suggests that the presence of bays (see Fig. 22, top, for examples of bay regions) and modifications to the duo CH groups, like adding or removing hydrogens or substituting a N for a CH group, tend to shift the intensity from 7.6 to 7.8 ${\rm\mu m}$. We also find that the broadening and merging of the 7.7 ${\rm\mu m}$  complex and the 8.6 ${\rm\mu m}$  bands for very large irregular PAHs as reported by Bauschlicher et al. (2009) holds in general for PAHs with irregular edge structures. On the other hand, if the molecule is more compact and the shift is caused by adding H or removing H then there is less broadening and hence less merging with the 8.6 ${\rm\mu m}$  band. However, not all such changes in structure shift the intensity in a significant manner. Hence, multiple factors are likely contributing to the intensity distribution in this range, and we may not have identified and/or quantified all of them.

The G8.2 ${\rm\mu m}$  component

No band assignments exist for the 8.2 ${\rm\mu m}$  emission as it had previously not been identified as an individual feature in the PAH band family but rather was included in the plateau emission (i.e. 8 ${\rm\mu m}$  bump) underlying the 6.2 and 7.7 ${\rm\mu m}$  PAH bands when the local spline decomposition is used or in the 7.7 and 8.6 ${\rm\mu m}$  features in case of the global spline continuum and Lorentzian decomposition. We therefore investigated all PAHs present in the PAHdb for 8.2 ${\rm\mu m}$  emission and found that it seems to originate in C-H in-plane bending modes at bay sites. We further explored this possible origin based on the molecule C150H30. We have gradually modified the edge structure of this parent molecule to enhance the number of bay regions (see Fig. 22, top) and computed the cation spectrum for these new structures. The results are shown in Fig. 22121212 These will be added to the PAHdb in the next update.. The parent structure exhibits two strong emission bands near $\sim$7.8 and $\sim$8.6 ${\rm\mu m}$. By increasing the number of bay regions, the 7.8 ${\rm\mu m}$  band shifts towards 8.2 ${\rm\mu m}$  while the 8.6 ${\rm\mu m}$  component gets weaker and shifts to shorter wavelengths. Based on this analysis of the PAHdb, we therefore attribute the 8.2 ${\rm\mu m}$  emission to C-H in-plane (ip) bending modes in PAHs with multiple bay regions, as for example found in PAHs with very irregular structures or corners.131313Note that this also holds for neutral PAHs. However, neutral PAHs typically have very weak emission in this wavelength range and hence do not significantly contribute to the emission here. We do not find a dependence on PAH size: all PAHs with multiple bay regions regardless of size have emission at 8.2 ${\rm\mu m}$.

The G8.6 ${\rm\mu m}$  component

Bauschlicher et al. (2008, 2009) reported that the 8.6 ${\rm\mu m}$  PAH emission is due to the C-H in-plane bending mode produced by large compact, symmetric PAHs. Indeed, perusal of the PAHdb indicates that almost any change in symmetry reduces the intensity of the 8.6 ${\rm\mu m}$  emission. Typically, PAHs of 96 or more carbon atoms show a significant 8.6 ${\rm\mu m}$  band but – as noted by Ricca et al. (2012) – PAHs larger than 150 C-atoms show a broad complex of emission between 6 and 7 ${\rm\mu m}$  and hence they are not important in the ISM. For completion, we list here the small and medium sized PAHs in the PAHdb that do exhibit emission at 8.6 ${\rm\mu m}$: some PAH ions with 20 to 22 carbon atoms, C54H18 (0, 2+, 3+), and C66H19/20 anions. Keep in mind that PAHs with #C-atoms $<$ 30 are easily destroyed and not expected to be prevalent in the ISM (Micelotta et al., 2010; Berné & Tielens, 2012).

5.1.2 Relative intensities

The astronomical 8.6/6.2 and 7.7/6.2 intensity ratios are not well reproduced by the PAH intrinsic relative intensities for a collection of PAHs (Ricca et al., 2012). This discrepancy depends on the PAH structure: more elongated PAHs (e.g. oval PAHs) seem to fall closer to the observed values than circular PAHs (i.e. the coronene family). Here we find that the 7-9 ${\rm\mu m}$  PAH emission is due to at least two PAH subpopulations. These PAH subpopulations are not traced with the nominal main PAH bands which may explain the lack of overlap between the astronomical and intrinsic intensity ratios. Therefore, we calculated the intrinsic intensities in the range of 7.365 to 7.8 ${\rm\mu m}$  for the G7.6 component, of 7.8 to 8.08 ${\rm\mu m}$  for the G7.8 component, of 8.095 to 8.365 ${\rm\mu m}$  for the G8.2 component and of 8.408 to 8.752 ${\rm\mu m}$  for the G8.6 component, applying a 15 cm*-1* redshift and an excitation of 8 eV. This corresponds to the typical photon energy for the illuminating star. We normalized these intensities to the 6.2 ${\rm\mu m}$  band intensity. Note that changing the integration range for the 6.2 ${\rm\mu m}$  PAH band will change the calculated ratios. We therefore set this integration range to 6.2–6.6 ${\rm\mu m}$, corresponding to the range of frequencies of the most intense band in this wavelength range for all PAH ions with #C $\geq$ 20 (a total of 456 molecules, applying a 15 cm*-1* redshift). The results for the oval PAH ions are given in Table 4 and the results for all PAH ions in the PAHdb with #C $\geq$20 (including the oval PAHs presented here) are shown with the observations in Fig. 23. The spatial distribution and the line cuts of these intensity ratios are also shown in Figs. A.7 and A.8.

When comparing the theoretical results for the oval PAH molecules discussed in Section 4.2 (Table 4), we note that there is no clear pattern. The ratios tend to be larger for the anions than the cations. Within a single family, the $I_{G7.8}/I_{6.2}$ and the $I_{G8.6}/I_{6.2}$ intensity ratios as well as the $I_{G8.2}/I_{6.2}$ anionic intensity ratio do increase with size. A few molecules as well as a few averages are consistent with the observations but do not represent the bulk of the ratios found in NGC 2023. Considering the entire database, a large range of values is present for all four intensity ratios, which can not be attributed solely to the smallest PAHs in the PAHdb. The observed range of each intensity ratio individually is well covered by the intrinsic PAH ratios (Fig. 23). However, the intrinsic values do not cover well both the range in the observed G7.6/6.2 and in the observed G8.6/6.2 at the same time: indeed, only a handful of molecules (all anions) do exhibit both observed intensity ratios. The match is better for the G7.8/6.2 versus the G8.2/6.2 intensity ratios where several molecules coincide with the observed ratios (including anions, cations and dications). Note that no large PAHs showing strong 8.6 ${\rm\mu m}$  emission (#C $\geq$ 96, see previous section) matches the G7.6/6.2 versus G8.6/6.2 intensity ratios in NGC 2023. The C66H19/20 anions as well as C20H12/14 ions do agree141414These include C66H${}_{19}^{-}$, PAHdb uid of 712 and 711; C66H${}_{20}^{-}$, PAHdb uid of 605 and 116; C20H${}_{12}^{-}$, PAHdb uid of 398; and, to a lesser extent, C20H${}_{14}^{+}$, PAHdb uid of 370.. Indeed, these large PAHs (#C $\geq$ 96) exhibit emission near 7.8 ${\rm\mu m}$  instead of 7.6 ${\rm\mu m}$  in contrast with the observed trend between the 7.6 and 8.6 ${\rm\mu m}$  PAH bands. Such an issue does not arise with bay PAH cations responsible for the 8.2 component as they also emit at 7.8 ${\rm\mu m}$. Here, the complexity occurs for the CH bending modes: the frequency of the solo CH out-of-plane bending mode is longwards of 11.2 ${\rm\mu m}$, well in the realm of the traditional solo CH out-of-plane bending mode for neutral PAHs.

Since there are a few matches between the PAHdb and the observations, we can investigate if the observed trend in the G7.6/6.2 versus G8.6/6.2 and G7.8/6.2 versus the G8.2/6.2 intensity ratios (i.e. the slope in the observations) depends on PAH parameters such as molecular structure and PAH excitation. To explore possible dependence on PAH structure, we selected a PAH (C36H${}_{20}^{+}$, PAHdb uid of 542) that matches well with the G7.8/6.2 versus the G8.2/6.2 observational trend and checked the intensity ratios of related PAHs, i.e. PAHs of the same size but different structure (i.e. C36H${}_{16}^{+}$, PAHdb uid of 155 and C36H${}_{16}^{+}$, PAHdb uid of 129).151515While these are not relevant structures as discussed in Section 5.1.1, we were constraint by the limited match between the observations and pure PAH cations. These however do not match the observations suggesting (although based on a very small sample) that molecular structure is not the dominant driver of the observed trend. Likewise, we explored the effect of different excitation for both trends. We selected C66H${}_{20}^{-}$ (PAHdb uid of 116) and C66H${}_{20}^{+}$ (PAHdb uid of 595) for respectively the G7.6/6.2 versus G8.6/6.2 and G7.8/6.2 versus the G8.2/6.2 observational trend and calculated their intensity ratios for excitation energies of 4, 6, 8 and 10 eV. An increase in excitation energy results in decreased band ratios. While a change in excitation energy mimics both observed trends well (Fig. 23), the implied excitation energy variations are much larger than expected in the bright zones of a PDR. Thus, while variations in the excitation energy could be one of the parameters contributing to the observed trends, it is unlikely to be the dominant one.

Despite the fact that the PAHdb covers a wide range in relative intensities, only a small number of PAHs have relative intensities that match those observed in NGC 2023. This can be due to (a combination of) various reasons including:

1. (i)

only a small fractions of PAHs are responsible for these 4 Gaussian components, and thus the 7.7 and 8.6 ${\rm\mu m}$  PAH emission, and/or 2. (ii)

the PAHdb is incomplete and lacks the necessary PAH molecules to match the observations, and/or 3. (iii)

the uncertainty on the relative intensities is considerable and inhibits a meaningful comparison with the observations (see Section 4.1 for a detailed discussion), and/or 4. (iv)

instead of single PAH molecule obeying the observed trends, a selection of PAHs fulfills the observational requirement as a group.

Point (i) is consistent with the grandPAH hypothesis: only a set of the most stable PAH species can survive in the harsh conditions of the ISM and are responsible for the observed PAH emission (e.g. Tielens, 2013; Andrews et al., 2015). This also implies however that these (or some of these) grandPAHs are currently not present in the PAHdb as it does not provide a good match with the observations, in particular for the G7.6 and G8.6 subcomponents. However, we note that the PAHdb now does comprise all large, compact, symmetric PAHs (i.e. circular and oval PAHs with 96 $<$ #C-atoms $<$ 150). These are also the most stable ones, an essential aspect to the grandPAH hypothesis. The lack of a good match then may pose a challenge to this hypothesis.

Point (ii) is not that surprising as the PAHdb is known to be incomplete and, although biased towards smaller PAHs, is still improving (Bauschlicher et al., 2010; Boersma et al., 2014). The current incompleteness and bias may prove to be problematic for the 8.6 ${\rm\mu m}$  emission and it’s strong connection to the 7.6 ${\rm\mu m}$  emission. Indeed, the large, compact, symmetric PAHs that exhibit a strong 8.6 ${\rm\mu m}$  emission also have emission at 7.8 instead of 7.6 ${\rm\mu m}$. In addition, any reduction of symmetry decreases the 8.6 ${\rm\mu m}$  emission. This may then indicate that they are not responsible for the 8.6 ${\rm\mu m}$  emission but instead (a few) smaller PAHs are (as discussed in Section 5.1.1). Within this hypothesis, these structure as well should have been missed in the PAHdb. Indeed, only C66H19/20 and C20H12 anions have G7.6 and G8.6 relative intensities well within the realm of the observations while two PAH cations (C20H${}_{12}^{+}$; C20H${}_{14}^{+}$) straddle the observed boundaries for these relative intensities. It’s interesting though that some of these smaller sized species do emit at 7.6 ${\rm\mu m}$.

Finally, we note that excellent fits to astronomical data with the PAHdb have already been obtained (e.g. Cami, 2011; Boersma et al., 2013; Andrews et al., 2015). This may imply that the PAHs may not satisfy the observed trends individually but do so rather as a group. This would however require that all members of the ‘group’ behave spatially the same under a variety of physical conditions as e.g., the G7.6 and G8.6 exhibit the best correlations found between the PAH emission components. As a corollary, photochemistry is very sensitive to excitation and hence size, implying that these PAHs must all have very similar sizes. However, this seemingly contradicts with the assignment of the 7.6 ${\rm\mu m}$  emission to PAHs with $\#$ C-atoms $<$ 100 and the 8.6 ${\rm\mu m}$  emission to PAHs with $\#$ C-atoms $\geq$ 96.

Hence, it is clear that while none of these arguments undermines the PAH and/or grandPAH hypotheses, it does indicate that further investigations and possibly fine-tuning of the hypotheses are warranted.

5.1.3 CH modes

The spatial distribution of the 8.6, 11.0 and 17.4 ${\rm\mu m}$  PAH emission is remarkably similar but quite distinct from that of, for example, the 6.2 and 11.2 ${\rm\mu m}$  PAH emission (Section 3). Within the framework of the photochemical evolution of PAHs presented in Section 5.1.1, this then implies they arise from similar PAH subpopulations. However, experimental and theoretical studies of PAHs assign these bands to different PAH subpopulations. Indeed, the 8.6 ${\rm\mu m}$  PAH emission is due to CHip bending modes in large symmetric PAH ions and/or some smaller PAH ions (Section 5.1.1, Bauschlicher et al., 2008, 2009). The 11.0 ${\rm\mu m}$  PAH emission generally is attributed to solo CH out-of-plane (oop) bending mode in (any) large PAH cation (Hony et al., 2001; Bauschlicher et al., 2008, 2009). Candian & Sarre (2015) argue instead for a possible origin in neutral acenes. For NGC 2023, Shannon et al. (2016) report a neutral contribution of roughly 8–16% to the 11.0 ${\rm\mu m}$  band limiting the possible acene contribution in this object. Finally, the 17.4 ${\rm\mu m}$  PAH emission is systematically seen in large PAHs of all charges (Boersma et al., 2010; Ricca et al., 2010). These three “attributions” represent quite different PAH subpopulations: going to longer wavelengths, the size of the PAH population that can contribute to these bands increases significantly. Yet, they have identical spatial morphologies. We can therefore put further constraints on their assignments based on their observational properties. As argued by Peeters et al. (2012) and Shannon et al. (2015), the observational characteristics of the 17.4 ${\rm\mu m}$  PAH emission points towards a carrier of PAH ions rather than neutrals. In addition, based on their chemical model calculation, Boschman et al. (2015) report that PAH cations are largely dehydrogenated. In this case, these effects combined reduces its PAH sub-population similar to that for the 11.0 ${\rm\mu m}$  PAH band. On the other hand, rather few PAHs (present in the PAHdb) have strong 8.6 ${\rm\mu m}$  emission and thus the 8.6 ${\rm\mu m}$  emission possibly imposes very stringent restrictions to the CH ionic PAH population (and thus the PAH population overall, point (ii) above). Alternatively, while the photochemical evolution of the PAH family in PDRs as described in Section 5.1.1 is consistent with the data, other interpretations are possible as well. Indeed, further investigations are needed as the simple assumption that a small collection of 8.6 ${\rm\mu m}$  carriers implies a limited set of carriers for the 11.0 and 17.4 ${\rm\mu m}$  band as well may not hold. For example, consider large symmetric PAHs that emit strongly at 8.6 ${\rm\mu m}$. These PAHs are dominated by solo CH groups with a few duo CH groups present. In addition, their solo CHoop intensity is even stronger than one would expect since the solo CH-groups steal intensity from the duo CH-groups (Bauschlicher et al., 2008; Ricca et al., 2012). In this case, the 11.0 ${\rm\mu m}$  emission is dominated by these large symmetric PAHs minimizing the effect of the larger PAH subpopulation contributing to the 11.0 ${\rm\mu m}$  band compared to the PAH subpopulation contributing to the 8.6 ${\rm\mu m}$  band. It is yet unknown whether a similar effect is present for the 17.4 ${\rm\mu m}$  band.

Another relatively strong CH band is found at 12.7 ${\rm\mu m}$. Generally, this band has been ascribed to CHoop bending vibration of duo and trio CH groups in (very) large PAHs (Hony et al., 2001; Bauschlicher et al., 2008, 2009; Ricca et al., 2012). In addition, the 12.7 ${\rm\mu m}$  intensity is further enhanced by the coupling between the CHoop mode and the C–C ring deformation mode in (elongated) armchair PAHs (Candian et al., 2014). Thus the assignment for the 12.7 ${\rm\mu m}$  band, as well as for the G7.8 and G8.2 components, indicate the importance of bay regions, yet their spatial morphologies are very distinct. Indeed, we attribute the G7.8/G8.2 components to large PAHs with bay regions. As the smaller ones of these ($\leq$ 70 C-atoms) are broken down, the largest and most stable ones remain abundant the longest (e.g., closest to the star). Theoretical studies have shown that carbon loss from compact PAHs lead to the formation of armchair structures (Bauschlicher & Ricca, 2014). These surviving species are then the carriers of the 12.7 ${\rm\mu m}$  band, where we note that, for armchair PAHs, the CHoop mode couples with the C–C skeletal modes, resulting in an enhanced intrinsic strength of the 12.7 ${\rm\mu m}$  band Candian et al. (2014). Clearly, further laboratory studies on the photochemical evolution of PAHs are warranted.

5.2. Clusters and VSGs

5.2.1 PAH features versus plateaus

We conclude and confirm earlier reports (Bregman et al., 1989; Roche et al., 1989, paper I) that the plateaus are distinct from the features. As a consequence, we need to fine-tune applied spectral decompositions to obtain the PAH components. Several decomposition methods exist in the literature and they differ greatly with respect to the treatment of the broad underlying plateaus. These plateaus are either incorporated in the PAH bands themselves through the wings of Drude or Lorentzian profiles (e.g. Boulanger et al., 1998; Smith et al., 2007b; Galliano et al., 2008) or treated separately (e.g. Hony et al., 2001; Peeters et al., 2002; Galliano et al., 2008, this paper). Our findings, along with those of Bregman et al. (1989), Roche et al. (1989) and paper I, strongly imply the plateaus need to be treated independently from the features in spectral decompositions and argue against using Lorentzian and Drude profiles for the PAH features (see also Tielens, 2008).

5.2.2 Carrier of the plateaus

The spatial characteristics of the plateaus’ emission also reveal further information about their origin. The morphology of the 10–15 ${\rm\mu m}$  plateau in the north map suggests that its carrier may be more closely related to the carrier of the MIR dust continuum than to the PAH molecules responsible for the PAH features. Likewise, the 5–10 ${\rm\mu m}$  plateau emission has a different carrier from that of the 10–15 ${\rm\mu m}$  plateau, one that is more closely related to the PAH molecules, for instance by having a size between that of the carriers of the 10–15 ${\rm\mu m}$  plateau and that of the PAH bands. In contrast to the 5–10 and 10–15 ${\rm\mu m}$  plateau, the morphology of the 15–18 ${\rm\mu m}$  plateau resembles more that of the PAH emission. Paper I and Shannon et al. (2015) reported that it likely originates in large neutral PAHs. Indeed, while it has been found to vary independently of the 15–18 ${\rm\mu m}$  features overall, it does correlate well with the very weak 15.8 ${\rm\mu m}$  band. This is also the only band in this wavelength range that is attributed solely to neutral PAHs because of its correlation with other neutral charge proxies, while the other features have been assigned to either PAH ions or a mixture of neutral and charged PAHs (Paper I; Shannon et al., 2015).

A larger sized carrier for the 5–10 and 10–15 ${\rm\mu m}$  plateaus is consistent with previous suggestions. For instance, Allamandola et al. (1989) and Bregman et al. (1989) concluded that these broad components likely arise from a semi-continuum produced by a mixture of larger PAHs, PAH clusters and small amorphous carbon particles containing on the order of $\sim$ 400-500 carbon atoms. More recently, blind signal separation (BSS) analysis of spectral maps revealed three spatially distinct components contributing to the ‘nominal’ PAH emission, which are attributed to neutral PAHs, ionized PAHs and evaporating very small grains (PAH0, PAH+ and eVSGs respectively; Boissel et al., 2001; Rapacioli et al., 2005; Berné et al., 2007). Instead of applying a BSS analysis to our spectral maps, we have used its derivative PAHTAT (Pilleri et al., 2012) for comparison with our results (for a detailed discussion on PAHTAT and its application to our data, see Appendix H). The spatial morphology of the eVSG component and the plateau emission is very similar though small differences are present (Figures A.11). In hindsight, this is not so surprising as, spectroscopically, the eVSG component is quite distinct and considerably narrower and smaller compared to the plateau emission as defined in this paper (Figure A.10).

The similarity between eVSG and 8.1 ${\rm\mu m}$  extreme argues for a stronger contribution of the VSG to the 7-9 ${\rm\mu m}$  emission than what is found by the PAHTAT analysis and thus the BSS analysis (i.e. the fraction of the eVSG component to the total 7–9 ${\rm\mu m}$  emission, in particular near 8.1 ${\rm\mu m}$). Combined with the fact that we found that the PAH features behave independent of the underlying plateau emission, this suggests that the PAHTAT and BSS “decomposition” is in need of more components to decompose the PAH emission into its distinct PAH sub-populations.

5.2.3 Carrier of the G7.8 and G8.2 component

As noted, the eVSG emission is spatially most similar to that of the 8.1 ${\rm\mu m}$  extreme in the north map (Figure A.11). This suggests that the eVSG component is related to the G8.2 component despite the fact they are spectroscopically two distinct components (see Fig. A.10). In addition, the spatial distribution seen at 8.1 ${\rm\mu m}$, the G8.2 component and of the eVSGs is not seen at other wavelengths in the movies regarding the 6.2, 11.2 and 12.7 ${\rm\mu m}$  bands (Figure 10) nor is it seen for the G7.6 and G8.6 components (Figure 9). This is therefore inconsistent with a carrier like the typical PAHs that are responsible for the PAH emission bands.

Recent results on small neutral PAH clusters (Roser & Ricca, 2015) indicate that they have very similar infrared spectral characteristics in the 6-9 ${\rm\mu m}$  region (peak position) to those of their respective neutral PAH components. If we assume that this finding also holds for clusters of larger size and/or different charge state, we can extend the earlier proposed assignment for the G8.2 emission based on the PAHdb (see Section 5.1.1) to PAH clusters: PAHs and PAH clusters with multiple bay regions exhibit emission near 8.2 ${\rm\mu m}$, independent of their charge. In addition, nanograins and/or very small amorphous carbon particles could also possibly contribute to the G8.2 component. However, this cannot be currently confirmed due to the limited amount of infrared data available. An origin in larger sized species (PAH clusters, nanograins and/or very small amorphous carbon particles) is also more consistent with its spatial distribution being more similar to that of the eVSG and continuum emission than that of typical PAH bands. This can be further extended to the G7.8 component which has a similar morphology. While spectroscopically, the eVSG component obtained with the BSS method is quite distinct and considerably broader and weaker compared to the G7.8 emission as defined in this paper (Figure A.10), it’s worth noting that Rapacioli et al. (2005); Berné et al. (2007) attributed the 7.8 ${\rm\mu m}$  subcomponent to eVSGs. Note that the ‘traditional’ 7.7 ${\rm\mu m}$  component has contributions of both the G7.6 and G7.8 component as defined in this paper and thus certainly has a contribution of PAH molecules.

5.3. Profile classes

Peeters et al. (2002) classified the profiles of the ‘nominal’ 7.7 and 8.6 ${\rm\mu m}$  PAH bands in class A, B and C. We applied our 4 Gaussian decomposition to the same sample of ISO-SWS observations as these authors. The class A profiles are well fitted (see also Stock & Peeters, 2016) as well as class B profiles which are not redshifted. However, our decomposition fails to fit the majority of the class B profiles: indeed many class B profiles are redshifted compared to class A profiles and our decomposition (with fixed peak positions and FWHM) does not allow for a red-shifting of the entire band. If we let the peak position of the 4 Gaussian components vary slightly while keeping the same FWHM, we can reproduce all class B profiles. As expected, class B profiles have a lower G7.6/G7.8 intensity ratio compared to class A profiles. Moreover, the G7.6/G7.8 intensity ratio decreases with peak position of the G7.8 component (when allowing for slight red-shifting) albeit with large scatter. Given the assignments discussed earlier, the band profiles thus reflect different PAH populations and size distributions with class B sources.

5.4. Tracer for PAH charge

The 6.2, 7.7, 8.6 and 11.0 ${\rm\mu m}$  PAH bands have traditionally been used as a tracer for PAH ions and quantitative relations have been deduced to relate the 6.2/11.2, 7.7/11.2 or 8.6/11.2 intensity ratio with the ionization parameter and thus the physical conditions (i.e. the radiation field, G0, the electron density, ne, and the gas temperature, Tgas; e.g. Galliano et al., 2008; Fleming et al., 2010; Rosenberg et al., 2011; Boersma et al., 2015). Based on a BSS analysis, Rosenberg et al. (2011) argue that the 11.0 ${\rm\mu m}$  band may be better suited as a tracer for ions as it is a largely cationic band while the 6.2, 7.7, and 8.6 ${\rm\mu m}$  features include a mixture of PAH+ and PAH0 components. The results presented in this paper are partially consistent with this. Clearly, the ‘nominal’ 7.7 ${\rm\mu m}$  band traces at least two PAH subpopulations, one of which being PAH cations. By virtue of its strong correlation and spatial resemblance with the ‘nominal’ 7.7 ${\rm\mu m}$  band, the 6.2 ${\rm\mu m}$  band likely is also contaminated albeit to a lesser degree. To the contrary, the 8.6 and 11.0 ${\rm\mu m}$  bands behave similarly on a spatial scale and are located closest to the illuminating star indicating that both bands would work well as a proxy for PAH charge (we present the spatial distribution of the 8.6/11.2 in Appendix E). The reliability of the 6.2, 7.7 and 8.6 ${\rm\mu m}$  bands as a charge proxy was recently investigated by Boersma et al. (2015). The fitting approach of these authors recognizes that both neutral and charged PAHs contribute to these bands. Nevertheless, these authors find that the 6.2 and 8.6 ${\rm\mu m}$  band are indeed a good tracer of PAH charge while the 7.7 ${\rm\mu m}$  band is not. This has been attributed to a rapid growth of the 7.8 ${\rm\mu m}$  component that is coincident with a stronger red wing for the 11.2 ${\rm\mu m}$  band (traditionally assigned to neutral PAHs in the denser regions of the PDR). Comparing the applied decomposition methods, it is clear that their 7.7 ${\rm\mu m}$  band includes the G7.6, G7.8, G8.2 and a fraction of the G8.6 components and thus contains a larger contribution of the second PAH subpopulation responsible for the G7.8 and G8.2 components than a ‘nominal’ 7.7 ${\rm\mu m}$  band161616These authors extracted the 7.7 ${\rm\mu m}$  band in the following way (Boersma et al., 2014): a continuum similar to our global spline continuum is fitted (which we used for the 4 Gaussian decomposition) but the 7.7 ${\rm\mu m}$  band intensity includes emission underneath the 8.6 ${\rm\mu m}$  band represented, in our decomposition, by the difference between the local and global spline continuum (see Figure 2).. It is therefore not surprising that their 7.7 ${\rm\mu m}$  band intensity does not live up to expectations as a charge proxy. It’s noteworthy to re-iterate that the 7.7 ${\rm\mu m}$  band fails as a tracer for charge in regions that are also characterized by a stronger red wing for the 11.2 ${\rm\mu m}$  band. Rosenberg et al. (2011) have attributed this component to eVSGs, another indirect link that PAH clusters and/or VSGs may be the carrier of the G8.2 component.

5.5. Synopsis

The results presented in this paper, amongst others, indicate that we are in need to further fine-tune our ‘simple’ picture of the PAH population and its emission. While the charge balance of the PAH population is certainly the major factor driving the observed spatial variations of the major features, other parameters such as chemical structure and size drive more subtle, perhaps second-order, variability that we are now beginning to probe. Hence, our ‘picture’ of the PAH population in space should then reflect the variability of interstellar PAH structure and size, resulting in a spatially and continuously changing PAH population. Indeed, at any given location in the nebula, the ionic and neutral PAHs have the same chemical structure because the timescale for ionization and recombination is fast, and charge balance changes on the order of months while the timescale for chemical modifications or processing is much larger (105 yrs; Le Page et al., 2003; Berné & Tielens, 2012). However, with for example distance from the illuminating source, not only will the charge balance vary but also the chemical structure of the PAH population. Although the latter may happen on a larger spatial scale than the variation in charge balance. Likewise, the size of the PAHs present have been proven to vary with distance from the illuminating source (Croiset et al., 2015).

Clearly, this ever changing PAH population is reflected in the observed PAH emission inviting us to go beyond the traditional approach of investigating the ‘nominal’ PAH features. As shown, even within individual PAH emission bands, the spatial morphology varies with wavelength (Section 3.1.2). One should not think of a PAH band as a single emission feature. Instead, each PAH band is a combination of different components that each respond to the local physical conditions in their own way; they thus each (or their sub-population) have their own spatial distribution. This of course holds for the G7.6 and G8.6 Gaussian components which show a 1:1 correlation and have the highest correlation coefficient (0.987). The animation indicates that even at wavelengths where they are the only Gaussian components contributing (from the four used for the decomposition), the morphology changes continuously. Given that the PAH emission features are due to a collection of PAH molecules, this should not be surprising. In addition, these observations clearly indicate that there are several PAH subpopulations beyond the neutral and charged PAHs. This then softens the grand-PAH hypothesis in which a few extremely stable PAH molecules make up the PAH emission in space. It would at least imply that instead of one grandPAH set, at each location with unique physical conditions and possibly unique history, a distinct grandPAH collection is present.

6. Conclusions

Two Spitzer-IRS spectral maps from 5–15 ${\rm\mu m}$  of the northern and southern part of the reflection nebula NGC 2023 are presented. The spatial distribution of the individual emission components including the dust continuum emission, the H2 emission and the PAH emission, and relative PAH intensity correlations are investigated. These observational analyses resulted in the following conclusions:

• •

PAH emission bands exhibit a variety of different spatial morphologies: their (peak) emission occurs at different distances from the illuminating source revealing a ‘spatial sequence’ within the PAH emission.

• •

Multiple correlations are present between individual PAH bands, including the well-known correlation between the 6.2, 7.7 and 8.6 ${\rm\mu m}$  bands.

• •

Despite being well correlated, some PAH bands exhibit distinct spatial morphologies, in particular the traditional charge proxy bands at 6.2, 7.7 and 8.6 ${\rm\mu m}$. Hence, spectral maps reveal very important subtleties missed by correlation plots.

• •

The 7–9 ${\rm\mu m}$  PAH emission is decomposed into 4 Gaussian components (G7.6, G7.8, G8.2 and G8.6). The G7.6 and G8.6 components exhibit the same spatial distribution and this closely resembles that of the 11.0 ${\rm\mu m}$  PAH emission. These are the strongest correlated of all PAH components. The G7.8 and G8.2 components have spatial distributions similar to that of the dust continuum emission and not the PAH features emission (in case of the north map) and are not as tightly correlated. Thus, at least two PAH subpopulations with different spatial distributions contribute to the 7–9 ${\rm\mu m}$  PAH emission.

• •

The morphology of all the major PAH bands varies with wavelength indicating that multiple components contribute to a single feature. In particular, we find that the spatial distribution of the 7 to 9 ${\rm\mu m}$  PAH emission continuously vary between two extremes which are bound between $\sim$7.35 and $\sim$8.1 ${\rm\mu m}$.

• •

The PAH features also behave independently from the underlying plateau emission. Existing spectral decompositions clearly need further fine-tuning to treat these components separately.

We compared our observational results with those of other analysis methods such as PAHFIT and PAHTAT. We present spectra of compact oval PAHs ranging in size from C66 to C210, determined computationally using density functional theory, and investigate trends in the band positions and relative intensities as a function of PAH size, charge and geometry within this family.

Based on the entire NASA Ames database (PAHdb), we attribute the 7.6 ${\rm\mu m}$  emission to compact PAHs with sizes in the range of 50 to 100 C-atoms, the 7.8 ${\rm\mu m}$  emission to very large PAHs (100 $\leq\#$ C $<$ 150) with bay regions or modified duo CH groups (like adding or removing hydrogens or substituting a N for a CH group), the 8.2 ${\rm\mu m}$  emission to PAHs/PAH clusters with multiple bay regions (e.g. PAHs with very irregular structures or corners), and the 8.6 ${\rm\mu m}$  emission to very large compact, symmetric PAHs (96 $\leq\#$C $<$ 150). Within these framework, the observed spatial sequence then reveals the photochemical evolution of the interstellar PAH family with distance from the illuminating star.

We compared observed relative intensities with intrinsic relative intensities (as present in the PAHdb) in the 7 to 9 ${\rm\mu m}$. The astronomical range in individual relative intensities is well represented by PAHs in the PAHdb but very few are consistent with the observed trends (correlation of G7.6 with G8.6 and of G7.8 with G8.2). We highlight possible reasons for this apparent inconsistency.

We discuss the assignments of the CH bands within the framework of the presented photochemical evolution of interstellar PAHs. Based on the observed morphologies, we argue that the plateau emission and likely also the G8.2 ${\rm\mu m}$  component originate from larger sized species such as larger-than-typical PAHs, PAH clusters and/or very small grains. Finally, we discuss the astronomical implications of our results with regards to PAH charge proxies and the PAH population in space.

The authors thank the referee Kris Sellgren for the careful reading of the manuscript and very useful feedback. The authors are extremely grateful to Henry Leparskas for creating Figure 21. EP thanks Dr. D. Stock for fruitful discussions and feedback on the manuscript. We gratefully acknowledge support from the NASA Spitzer Space Telescope General Observer Program. EP gratefully acknowledges sustained support from the Natural Sciences and Engineering Research Council of Canada (NSERC: Discovery grant and Accelerator grant). LJA gratefully acknowledges support from NASA’s Astrobiology Program and Laboratory Astrophysics Carbon in the Galaxy consortium. AR thanks NASA’s Laboratory Astrophysics program grant NNX11AK09A for its generous support of this work. Studies of interstellar chemistry at Leiden Observatory are supported through advanced-ERC grant 246976 from the European Research Council, through a grant by the Dutch Science Agency, NWO, as part of the Dutch Astrochemistry Network and through the Spinoza premie from the Dutch Science Agency, NWO.

Appendix A Observed wavelength shifts

We’ve noticed an apparent small (of the order of a few $\times$ 0.01 ${\rm\mu m}$) wavelength shift in some 2x2 SL spectra, in particular for the south map. For the north map, only single-pixel spectra from the majority of pixels (not all) in the top row exhibit a wavelength shift. This affects the 2x2 spectra incorporating pixels of the top rows in two ways: 1) these 2x2 spectra will exhibit a shift and 2) the intensity at each wavelength is influenced in particular in wavelength regions with a large intensity gradient. In contrast, for the south map, the single-pixel spectra shift from row to row. This is known to occur for maps made with cross-slit step sizes larger than 1 pixel (as is the south SL map) and is a result of an under-sampling of the PSF (Spitzer helpdesk & Dr. Bregman, private communication). This effect is generally mitigated by using 2x2 extraction windows as the amplitude of the shift varies approximately between two values (one for even rows and one for uneven rows). However, because the shifts for even and uneven rows are approximately the same and not equal, the 2x2 spectra exhibit residual shifts (of a few $\times$ 0.01 ${\rm\mu m}$). Moreover, spectra in the second row from the bottom exhibit an opposite shift compared to those in other even rows. This will influence the two bottom rows in the 2x2 spectral maps. The effect of these wavelength shifts can clearly be seen in the peak position of the H2 emission lines (Fig. A.1). Moreover, it reveals an additional dependence with column171717The red column around x=50 in the left panel of Fig. A.1 is due to a spurious ’spike’ in the data on top of the H2 line which we were unable to remove. This features influences the peak position and derived strength of the H2 line (see also Fig. 5) and leads to a 12.7 ${\rm\mu m}$  detection less than 2 sigma (Fig. 5).. This issue may possibly be resolved by correcting the spectra for these wavelength shifts based on the H2 emission lines on a single pixel scale. However, not all SL2 spectra show H2 emission lines. In addition, the wavelength shift in the SL2 spectra seems to depend on wavelength since all features at the shorter wavelength (including the 6.2 ${\rm\mu m}$  PAH band) seem to shift while the wing of the 7.7 ${\rm\mu m}$  PAH band does not (although this is tricky to assess as the PAH bands are known to exhibit variations). Hence, we cannot correct this wavelength shift, which is negligible with respect to the width of the PAH features. While its effects can be noticed in the top row of the north map and in the bottom two rows of the south map in several feature intensity maps presented in this paper, it does not change the conclusions of this paper.

Appendix B The spatial distribution of the 8.6 and 11.0 ${\rm\mu m}$  emission.

For both the 8.6 and 11.0 ${\rm\mu m}$  maps, the maximum value for the colour table is determined by the peak emission in the SE ridge, which is located in the bottom two rows (in y-direction) and thus affected by the observed wavelength shift in the data (see discussion in Section 2.2 and Appendix A). This in turn affects the visual appearance of these maps. When the maximum of the intensity range for the colour table is determined from the S’ peak emission region (well represented by the pink contours in Fig A.2) both PAH maps exhibit an almost identical spatial distribution (see Fig A.2), given rise to their high correlation coefficient and 1:1 dependence.

Appendix C Results of PAHFIT decomposition

To investigate the influence of the decomposition method on the results, we applied PAHFIT (Smith et al., 2007b) to the SL data. PAHFIT models the spectra assuming a combination of starlight continuum, featureless thermal dust continuum, H2 emission lines, fine-structure lines, and dust emission features while taking into account dust extinction. The starlight continuum is represented by blackbody emission, the thermal dust continuum by a set of modified blackbodies, the H2 and fine-structure lines by Gaussian profiles and the dust emission features by a combination of Drude profiles. An example of the resulting fit and its separate components is shown in Fig. 2. Clearly, this decomposition is quite different from the spline decomposition and hence a direct comparison of the intensities and profiles of the individual emission components derived by the different methods is not possible. A major difference in these two decomposition methods is the fact that the underlying emission plateaus (as defined in the spline decomposition) are accounted for by (the broad wings of) the Drude profiles in the PAHFIT decomposition and hence are part of the individual emission features. Also, several PAH features/complexes are fitted by two or more Drude profiles. In particular, within PAHFIT, the ”7.7 ${\rm\mu m}$  Complex” is defined as the blend of the Drude profiles with $\lambda$ (fractional FWHM) of 7.42 (0.126), 7.60 (0.044) and 7.85 (0.053) ${\rm\mu m}$, the ”11.2 ${\rm\mu m}$  Complex” as the blend of Drude profiles with $\lambda$ (fractional FWHM) of 11.23 (0.012) and 11.33 (0.032) ${\rm\mu m}$, and the ”12.6 ${\rm\mu m}$  Complex” as the blend of Drude profiles with $\lambda$ (fractional FWHM) of 12.62 (0.042) and 12.69 (0.013) ${\rm\mu m}$. Finally, the dust continuum emission generally has a lower intensity in the PAHFIT decomposition compared to spline decomposition, influencing the resulting dust features’ intensities. Despite these significant differences, it has been shown that the overall conclusions on PAH intensity correlations for a large sample of objects are independent of the chosen decomposition approach (e.g. Uchida et al., 2000; Smith et al., 2007b; Galliano et al., 2008). Here, we investigate if this also holds when considering spectral maps of individual sources.

We set the following conditions when applying PAHFIT. First, we assumed no extinction is present, similar as for spline decomposition (see Section 2.4 for a discussion on extinction). Second, we set the contribution of star light to zero. The effective temperature of the illuminating star is 22,000 K and the resulting fits do not include starlight continuum emission when a contribution is allowed in the fitting routine. Third, we fixed the wavelengths and FWHMs for the dust features and lines (i.e. the default in PAHFIT). Fourth, we set the intensity of the H2 S(7), S(6) and S(4) to zero (at 5.511, 6.109 and 8.026 ${\rm\mu m}$  respectively). These H2 lines are not discerned in the spectra. Unfortunately, the H2 S(4) line is blended with the 7.7 ${\rm\mu m}$  PAH complex. For the south map, the results from PAHFIT includes emission from this H2 S(4) line when a contribution is allowed in the fitting routine. However, as mentioned, we do not discern its contribution in the spectra, even in the spatial positions where this lines is strongest. Furthermore, its spatial distribution shows emission in the S ridge (as do the H2 S(5) and S(3) lines) but it exhibits additional strong emission in the north-west region of the map where no emission is seen from the H2 S(5) and S(3) lines. Hence, we conclude that this is an artefact of the fitting routine and thus do not allow a contribution of the H2 S(7), S(6) and S(4) lines.

The dust continuum emission is fitted by only two and exceptionally three of the maximum allowed eight thermal dust components (with T of 135, 200 and 300 K). Typically, either the 200K or 300K component is used in combination with the 135K component. The 300K component is used towards the south ridge and sources C and D while the remaining positions require a 200K component. This results however in quite distinct dust continuum intensities at the shorter wavelengths (see Fig.  A.3): the continua for the spectra fitted with a 300K dust component are very close to the observed intensities at for example 5.5 ${\rm\mu m}$. In contrast, the continua for the spectra fitted with a 200K dust component have very little intensity at the shorter wavelengths. This distinction in continuum intensity largely reflects the constraints of the fitting routine rather than a physical origin. However, a smaller dust continuum intensity requires a larger contribution of the (wings of the) Drude profiles in order to fit the data. Hence, the strength of the PAH components is influenced by the fitted dust continuum, particularly for the weaker bands thus inhibiting any useful interpretation of these weaker bands. The spatial distribution and line cut of the fitting components as well as these PAH complexes are shown in Figs. A.4, A.5, and A.6. A comparison with the results of the spline decomposition (see Sect. 3.1 and Figs. 5, 6, 9, 10, and 16) reveals the following:

• •

The spatial distribution of the 6.2 ${\rm\mu m}$  PAH emission is very similar for both decompositions. Their small discrepancies between each other can be understood by the different spatial distribution of the LS 6.2 ${\rm\mu m}$  PAH band and the LS 5-10 ${\rm\mu m}$  plateau which both contribute (with different degrees) to the PAHFIT 6.2 ${\rm\mu m}$  PAH band.

• •

The morphology of the PAHFIT 7.7 ${\rm\mu m}$  Complex is slightly distinct from that of the LS 7.7 ${\rm\mu m}$  PAH band arising from the (partial) contribution of the GS 5-10 ${\rm\mu m}$  plateau emission to the PAHFIT 7.7 ${\rm\mu m}$  Complex and its components. Specifically, the PAHFIT 7.7 ${\rm\mu m}$  Complex shows a decreased relative intensity in the horizontal filaments in the northern area of the south map and the broad, diffuse plateau north of the S - SSE ridges. For the north map, it shows an extension towards the N ridge but a decreased extension towards the west of the NW ridge. Likewise, the PAHFIT 7.6, 7.8 and 8.3 ${\rm\mu m}$  components resemble the G7.6, G7.8 and G8.2 components respectively. There is slightly more deviation in the south map for the 7.6 and 7.8 ${\rm\mu m}$  components: the PAHFIT 7.6 component shows decreased relative intensity in the SSE and S ridges compared to the G7.6 ${\rm\mu m}$  emission and the PAHFIT 7.8 ${\rm\mu m}$  emission shows a very small increased relative intensity in the S’ ridge. The PAHFIT 7.4 ${\rm\mu m}$  component follows the PAHFIT 8.3 morphology in the south map. In contrast, in the north map, it shows a unique morphology.

• •

The PAHFIT 8.6 ${\rm\mu m}$  emission spatially resembles the G8.6 ${\rm\mu m}$  emission (GS). Their discrepancies arise from the (partial) contribution of the GS 5-10 ${\rm\mu m}$  plateau emission. Nevertheless, the PAHFIT 8.6 ${\rm\mu m}$  morphology shows the same trends as the LS 8.6 and GS G8.6 ${\rm\mu m}$  bands: 1) lack of intensity in locations where the 11.2 ${\rm\mu m}$  emission peaks (S ridge and center/north of NW ridge) and 2) it is distinct from the morphology of the 7.8 and 8.2/8.3 ${\rm\mu m}$  components while being similar to the spatial distribution of the 7.6 ${\rm\mu m}$  component.

• •

Similar to the 6.2 ${\rm\mu m}$  PAH band, the LS 11.2 ${\rm\mu m}$  band and the PAHFIT 11.2 ${\rm\mu m}$  Complex exhibit a similar spatial distribution. The PAHFIT 11.2 ${\rm\mu m}$  Complex includes part of the LS 10-15 ${\rm\mu m}$  plateau. The latter’s spatial morphology is distinct from that of the LS 11.2 ${\rm\mu m}$  PAH band in the north map resulting in small discrepancies between the LS and PAHFIT 11.2 ${\rm\mu m}$  band/complex. Within PAHFIT, the 11.2 ${\rm\mu m}$  Complex is composed of bands positioned at 11.2 and 11.3 ${\rm\mu m}$  which exhibit distinct spatial morphology. Specifically, the 11.3 ${\rm\mu m}$  band emission is most similar to the 11.2 ${\rm\mu m}$  Complex emission (it contributes the largest fraction to the complex). On the other hand, the morphology of the 11.2 ${\rm\mu m}$  band deviates from that of the 11.2 ${\rm\mu m}$  Complex and is an intermediate step in a progression towards that of the 7.6 and 8.6 ${\rm\mu m}$  bands.

• •

The morphology of the PAHFIT 12.7 ${\rm\mu m}$  Complex emission is similar to that of the PAHFIT 11.2 ${\rm\mu m}$  Complex emission and that of the LS 11.2 ${\rm\mu m}$  emission while it is distinct from that of the LS 12.7 ${\rm\mu m}$  emission. This distinction arises due to the significant contribution of the LS 10-15 ${\rm\mu m}$  plateau to the PAHFIT 12.7 ${\rm\mu m}$  Complex. Similar to the 11.2 ${\rm\mu m}$  Complex, the two components (positioned at 12.62 and 12.69 ${\rm\mu m}$) accounting for the 12.7 ${\rm\mu m}$  Complex have distinct spatial distributions. Given that the 12.62 ${\rm\mu m}$  component contributes on average 11 $\pm$ 6% and 16 $\pm$ 2% to the 12.7 ${\rm\mu m}$  Complex for the south and north map respectively, the 12.7 ${\rm\mu m}$  Complex emission resembles closely the 12.62 ${\rm\mu m}$  emission. In contrast, the 12.69 ${\rm\mu m}$  emission shows a spatial distribution between that of the LS 12.7 ${\rm\mu m}$  emission and that of the LS 11.0 ${\rm\mu m}$  emission.

Appendix D Fitting results

The parameters of the fits to the observed correlations and their Pearson’s correlation coefficients are listed in Table A.1.

Appendix E PAH ratio maps

Line cuts of PAH ratios shown in Figs. 7, 15, and 23 are presented in Fig.A.7. A selected number of PAH ratio maps are presented in Fig. A.8.

The first set shows the fraction of the 7.7 ${\rm\mu m}$  feature due to either the G7.6 or G7.8 component. As we use the global spline (GS) continuum before applying the four Gaussian decomposition, we also use this continuum for the 7.7 ${\rm\mu m}$  feature in these maps. Clearly, the G7.6 component makes up the dominant fraction of the 7.7 ${\rm\mu m}$  complex. The spatial distribution of the G7.6/7.7(GS) is closest to, but not similar, that of the 11.0 ${\rm\mu m}$  emission.

Boersma et al. (2014) showed that peak emission of the ratio of the 7.6 and 7.8 ${\rm\mu m}$  subcomponents is parallel with the PDR front and is located between the PDR front and the star in NGC 7023. We report a similar result for NGC 2023: G7.6/G7.8 clearly peaks in regions where ionized PAHs reign (as traced by the 8.6 and 11.0 ${\rm\mu m}$  PAH bands) and shows minimal emission on the PDR front (as traced by H2) and in the molecular cloud.

Given the results of this paper, the ion-to-neutral PAH fraction is best traced by the 8.6/11.2 or 11.0/11.2 PAH ratio. Since both ratios have similar morphologies, Fig. A.8 only shows that of the 8.6/11.2 PAH ratio. This ion-to-neutral PAH fraction is smallest in all ridges in both the N and S maps and thus also in the PDR fronts, and peaks closest to the illuminating star.

Finally, we show the spatial distributions of the ratios presented in Fig. 23: the four Gaussian components normalized on the 6.2 ${\rm\mu m}$  PAH band. Given the charge attribution of the G7.6, G8.6 and 6.2 PAH features, the G7.6/6.2 and G8.6/6.2 shows little variation across the FOV. This is less clear in the N map, originating in the fact that the range of the color table is set by the edge columns (in x coordinates). Nevertheless, the weak variation in pixels not located at the edges is present in both maps. Note the clear decrease in emission on the PDR fronts of the S and SSE ridges in the south map, another indication that the 6.2 ${\rm\mu m}$  PAH band is not solely due to ions.

Appendix F Influence of the decomposition on the plateaus’ behavior

While the applied decomposition may well be flawed to some extent, the observed spatial differences of the plateaus and the PAH features nevertheless imply that the bulk of the emission in the plateaus behaves independently from that of the PAH features (Section 5.2). Indeed, we applied a linear fit to the dust continuum (cf. plateau continuum) which is consistent with the observed dust continuum emission in PDRs (e.g. Berné et al., 2007; Compiègne et al., 2008). The uncertainty on the dust continuum emission is therefore such that it will only have a relatively small influence on the strength of the plateau emission and hence does not change its spatial distribution significantly. The disentangling of the plateau emission and the features is more challenging. In this regard, the ‘in-between’ morphology of the 5–10 ${\rm\mu m}$  plateau may possibly indicate that we have not selected the ‘true’ plateau emission, though there is no way to tell based on this dataset. But we can turn the problem around: if both the features and the underlying plateau arise from the same PAH population, a similar spatial behavior would occur, independently on how the division is made between these individual components. Clearly, within the PAH population several subpopulations exist (the most obvious one is neutral vs charged PAHs) and so spatial differences are present within the PAH features themselves (as discussed in this paper and many others, e.g. Joblin et al., 1996; Sloan et al., 1999; Hony et al., 2001; Galliano et al., 2008). However, none of the ‘nominal’ PAH bands have their strongest emission in the NW ridge in the north map as do the 5–10 and 10–15 ${\rm\mu m}$  plateaus, indicating that the plateaus have a different carrier. A caveat is in place here for the 5–10 ${\rm\mu m}$  plateau. In this paper, we showed that the ‘nominal’ 7–9 ${\rm\mu m}$  PAH emission is comprised of at least two sub-populations, one of which does show strong emission in the NW ridge (as the continuum emission; cf. G7.8 and G8.2, Fig. 9). Even more so, the morphology of the 7–9 ${\rm\mu m}$  PAH emission changes with each wavelength and varies between two extremes, one peaking in the NW ridge and one peaking in the N ridge (as the ‘nominal’ PAH emission; Fig. 10). Fig. 11 clearly shows though that even in the $\sim$ 7.7–8.3 ${\rm\mu m}$  wavelength range, the features vary independently of the plateau emission. This ‘complication’ does not hold for the other ’nominal’ PAH bands (e.g. 6.2 and 11.2 ${\rm\mu m}$  PAH bands) and thus does not affect the 10–15 ${\rm\mu m}$  plateau. We thus conclude and confirm earlier results (Bregman et al., 1989; Roche et al., 1989, paper I) concluding that the plateaus are distinct from the features.

Appendix G The CH stretching Vibrations (3.2 - 3.3 $\mu$m) of oval, compact PAHs

The band position maxima for the C-H stretch of the molecules shown in Figure 17 is summarized in Table A.2. We first note that the trends for all four families are similar and follow those noted previously (Bauschlicher et al., 2008, 2009; Ricca et al., 2012). In summary, there is a slow increase in band position with increasing size for the neutral and cation and a slow decrease for the anions. For all sizes, the band position of the cations is always at the shortest wavelengths, followed by that of the neutrals and subsequently that of the anions. As with other large PAHs, the intensity of the cations is slightly weaker than that of the neutrals while the anion intensities are about twice as large. The intensity per H increases with size for all charge states, with the increase being the smallest for the anion and largest for the cation. As a result, the cation intensity becomes rather sizeable for the largest species.

The C-H stretching bands, plotted using a FWHM of 15 and 30 cm*-1*, are shown in Figure A.9. All of the species studied have 12 duo hydrogens, with the number of solos varying with size and with family class (see Table 2). Since the duo C-H stretches and the solo C-H stretches fall at slightly different wavelengths, the constancy of the C-H stretching frequency, when plotted with a 30 cm*-1* bandwidth, might appear to be a bit odd. To investigate this in more detail, we also plot the C-H stretching region using a FWHM of 15 cm*-1*. The important factor to note is that while the number of duo hydrogens is constant, the intensity of the shorter wavelength duo C-H stretch increases as the intensity of the longer wavelength solo C-H component increases. That is, there is a coupling of the solo and duo C-H stretching modes so that the intensities of both increase with increasing molecular size. This results in a much smaller change in the 30 cm*-1* FWHM band maximum than would have been expected if there was no coupling of the solo and duo C-H stretches.

Appendix H PAHTAT results

The BSS results (Boissel et al., 2001; Rapacioli et al., 2005; Berné et al., 2007) have been employed to develop PAHTAT, a fitting program that allows the extraction of the distinct components found by the BSS analysis when no spectral map is available (Pilleri et al., 2012). Specifically, PAHTAT fits individual spectra assuming a combination of four templates (PAH0, PAH+, PAHx and eVSGs), three minor PAH bands, gas lines and a featureless dust continuum, while taking into account dust extinction. As the PAHx template was only invoked to better fit highly excited regions (Joblin et al., 2008), we excluded this component from our fit. We obtain very good fits by using the PAHTAT tool (Fig. A.10). The largest discrepancy is found for the 6.2 ${\rm\mu m}$  band and is not resolved when including the PAHx template.

Fig. A.10 shows a comparison of the PAHTAT and the four Gaussian components for a typical spectrum. Spectroscopically, the eVSG component is quite distinct from both the plateau emission and any of the four Gaussian components as defined in this paper (Figure A.10). In contrast, the G7.6 and the G8.6 component is more similar to respectively the 7.7 and 8.6 component in the PAH0 and PAH+ templates.

By definition, the spatial morphology of the PAHTAT components (PAH0, PAH+, eVSGs) is unique (Figures A.11 and A.12) but they, of course, do resemble morphologies and line cuts presented in this paper. As expected, the morphology (and thus line cut) of the PAH0 component mimics the 11.2 ${\rm\mu m}$  PAH morphology, that of the PAH+ resembles the spatial distribution of the 8.6 and 11.0 ${\rm\mu m}$  PAH emission and that of the continuum component mimics the spatial distribution of the continuum emission. Interesting is the eVSG component. In the south map, its morphology is similar to that of the H2 emission, very sharply peaked on the S ridge unlike the plateau emission. In the north map, it is most similar to that of the 8.1 ${\rm\mu m}$  extreme. Compared to the plateau emission, it is more narrowly peaked in the N ridge and has considerably less emission in the NW ridge.