Size distribution of polycyclic aromatic hydrocarbons in space: an old new light on the 11.2/3.3 μm intensity ratio

Abstract

The intensity ratio of the 11.2/3.3 μm emission bands is considered to be a reliable tracer of the size distribution of polycyclic aromatic hydrocarbons (PAHs) in the interstellar medium (ISM). This paper describes the validation of the calculated intrinsic infrared (IR) spectra of PAHs that underlie the interpretation of the observed ratio. The comparison of harmonic calculations from the NASA Ames PAH IR spectroscopic database to gas-phase experimental absorption IR spectra reveals a consistent underestimation of the 11.2/3.3 μm intensity ratio by 34%. IR spectra based on higher level anharmonic calculations, on the other hand, are in very good agreement with the experiments. While there are indications that the 11.2/3.3 μm ratio increases systematically for PAHs in the relevant size range when using a larger basis set, it is unfortunately not yet possible to reliably calculate anharmonic spectra for large PAHs. Based on these considerations, we have adjusted the intrinsic ratio of these modes and incorporated this in an interstellar PAH emission model. This corrected model implies that typical PAH sizes in reflection nebulae such as NGC 7023 – previously inferred to be in the range of 50 to 70 carbon atoms per PAH are actually in the range of 40 to 55 carbon atoms. The higher limit of this range is close to the size of the C₆₀ fullerene (also detected in reflection nebulae), which would be in line with the hypothesis that, under appropriate conditions, large PAHs are converted into the more stable fullerenes in the ISM.

---

1. Introduction

The most widely accepted explanation for the emission of distinct infrared (IR) features, observed towards many different phases of the interstellar medium (ISM), is the PAH hypothesis.^1–3 From different regions of the ISM a small, but clear variation of spectral profiles is observed, which is attributed to a variation in the population of PAHs in terms of size, structure and charge state.^4,5 This variation is thought to be controlled by the local physical conditions.^6,7 Since the introduction of the PAH hypothesis, much effort has gone into modeling the mechanisms of molecular IR emission aiming to place better constraints on the molecular inventory of the ISM and to learn how they are affected by their environment. Apart from the charge state, abundance and structure, one of the constraints that receives a lot of attention is the size of PAHs, mostly in terms of the number of carbon atoms.

The IR emission by isolated PAHs starts with the absorption of a UV photon. Because under the dilute conditions of the ISM the collisional cooling rate is much smaller than the spontaneous emission rate, most energy is lost by vibrational relaxation which is accompanied by the fluorescence of IR photons. The observed IR spectrum is an average of various PAH species relaxing and emitting IR photons via different relaxation cascades.⁸ Interstellar PAHs are exposed to a spectrum of UV photon energies, but in first approximation and because the UV photon energy does not vary a lot within objects the UV photon energy is often taken as a single value.⁶

Currently, the commonly accepted approach for deriving the size distribution of PAHs involves the intensity ratio of the 3.3 μm band to that of a band at a longer wavelength, typically the 11.2 μm band.^2,9–16 The 3.3 μm peak corresponds to CH stretch vibrations and the relative intensity of this band was found to have the strongest (inverse) dependence on the number of carbon atoms within a PAH.⁹ The 11.2 μm peak is due to the CH out-of-plane (OOP) mode from solo Hs. The complete CH OOP region, which includes also the duo, trio and quartet H OOP modes, has been nicely decomposed and analyzed by Bauschlicher et al. and ranges from about 650–950 cm⁻¹ (15.4 to 10.5 μm).¹⁷

The main reason behind the sensitivity of the 11.2/3.3 μm intensity ratio to molecular size is that smaller PAHs reach a higher temperature upon absorption of the same UV photon due to their smaller heat capacity, and therefore emit more photons at higher frequencies than larger PAHs.⁴ The heat capacity, and thereby the cascade starting temperature that results from the absorption of a particular UV photon, can be described in terms of isolated harmonic oscillators or alternatively be calculated using more general PAH properties (approximation by Stein and Dwek et al.).^18–20 Additional advantages of using the 11.2/3.3 μm intensity ratio as a size tracer are that both bands originate mainly from CH bonds in neutral species and that they are spectrally well separated.⁶ Moreover, the influence of the number of hydrogen atoms or of the PAH molecular structure on the ratio is minimal.¹²

The NASA Ames PAH IR spectroscopic database (PAHdb) supports a number of analyses of the observed 11.2/3.3 μm ratio.²¹ It has been used to determine the PAH sizes in several astrophysical objects^12,22–25 and quite recently Maragkoudakis et al. explored extensively how the 11.2/3.3 μm ratio relates to the PAH size distribution.⁹ The calculated 0 K absorption intensities (hereafter called intrinsic strengths) that underlie the emission models not only affect the emission ratio, but also the band profile. By artificially adjusting the intrinsic ratio it was found that a large 11.2/3.3 μm ratio results in a red-wing at the 11.2 μm feature which fades away at smaller ratios.²⁶

The IR intensities in the PAHdb that are employed in the emission models are calculated at a relatively low level of theory (B3LYP/4-31G) and are based on the harmonic approximation. Because of the significant effect of intrinsic strengths on the outcome of the emission models in terms of ratio and peak profile, it is of key importance to validate these calculations with suitable experiments. Even more so because it was noticed already in the 90s by Langhoff and later by others that the intensity ratios derived at different levels of theory can vary significantly.^27–30 First, the use of a small basis set results in an overprediction of CH stretch intensities by about 20% when compared to matrix-isolation absorption spectra.²⁷ Secondly, it was recently shown that anharmonicity leading to phenomena such as combination bands and resonances has a dominant effect on the CH stretch region of the IR spectrum. These anharmonic contributions can carry significant intensity.^28,31

Moreover, also from an experimental perspective it has been concluded that one needs to be cautious. Compared to the gas phase, the CH stretch intensity in IR absorption measurements in the solid phase is about a factor 3 lower.³² In matrix-isolation IR absorption measurements, on the other hand, the CH OOP intensity is about a factor 5 lower. Both experimental conditions will thus have a large effect on the 11.2/3.3 μm ratio.³²

The aim of this study is to compare theory that is commonly used to interpret astronomical observations to experimental data in terms of intensity, utilizing direct absorption techniques for reliable intensities over a wide IR wavelength range. The theory at a (relatively) low vibrational temperature should be validated before extrapolating the models to describe the emission of highly excited PAHs and using these models to put constraints on the size distribution of PAHs. Gas-phase FTIR spectroscopy provides a good benchmark for oscillator strengths, because at the relatively low experimental temperatures there is no significant occupation of higher vibrational levels. In first approximation, a higher temperature will thus only lead to band broadening.^32,33 Moreover, with FTIR spectroscopy the data are collected simultaneously over a wide spectral range. We will show that the present analysis leads to the conclusion that previously employed emission models need to be reconsidered, and, importantly, that the application of modified models predicts notably smaller PAH sizes than before.

2. Methods

We compared (relative) intensities from the IR spectra of PAHs from four different experimental and theoretical sources, namely the NASA Ames PAH IR spectroscopic database (PAHdb, version 3.20),²¹ anharmonic calculations extracted from the work by Mackie et al.³⁴ or performed here using the Gaussian 16 software package, experimental gas-phase FT-IR absorption spectra from the NIST database³⁵ and experimental matrix isolation IR spectra extracted from the PAHdb.³⁶ To study the relation between molecular size and the intrinsic 3.3/11.2 μm ratio, we extracted theoretical PAH spectra from the PAHdb of neutral PAHs that do not contain any oxygen or nitrogen atoms (Fig. 1). The intensity is integrated between 650–950 cm⁻¹ (15.4 and 10.5 μm) and will be indicated in the following as ‘oop’. The integrated range for the CH stretch region, indicated in the following as ‘stretch’, spans from 2750–3250 cm⁻¹ (3.6–3.1 μm). For the ratio in emission (Fig. 3) we selected neutral PAHs that contain more than one hydrogen atom, no aliphatic groups and did not include Fe, Si, Mg, N, or O substituted species. Finally, we only considered PAHs with solo (more than 4) and duo hydrogens, as they are generally more compact and thus more stable. One can reasonably expect that this restriction will only affect to a minor extent the conclusions as the structure of the PAH should not have a large influence on the oop/stretch emission ratio.^27,32 Calculations referred to as anharmonic are either described previously and extracted from ref. 34, or they have been performed at the anharmonic^37,38 B3LYP³⁹/N07D^40,41 level of theory using Gaussian 16 (ref. 42) as indicated.

3. Results and discussion

3.1 Theory and gas-phase experiments

In order to investigate the intrinsic intensity ratio as a function of PAH size, Fig. 1a shows the band strength ratio extracted from the (harmonically calculated) PAH database for the complete available size range. Fig. 1a shows that this ratio varies by a factor of about 3 between the smallest and largest PAHs considered here. Interestingly, a large variation of the band strength ratio is associated with small PAHs containing up to about 35 carbon atoms, which results in a larger uncertainty in the calculated emission ratio over this size range. The color scale of Fig. 1a corresponds to the number of trio hydrogens in the PAH. A clear separation of PAHs that contain more than 4 trio hydrogens is visible from PAHs with 4 or less trio hydrogens. The presence of trio hydrogens increases the oop/stretch ratio by suppressing the CH stretch band intensity.

Fig. 1 (a) The intrinsic oop/stretch intensity ratio of neutral PAHs in the PAHdb as a function of the number of carbon atoms. The color scale corresponds to the number of trio hydrogens in the PAH. (b) Same plot as (a), but leaving out PAHs that contain more than 4 trio hydrogens (light red). The dark red symbols represent the average for fixed sizes (size tolerance is 5 carbon atoms). The black line is an exponential decay fit to the latter.

In the most stable PAH species, one does not expect many trio hydrogens as they are in general associated with an irregular shape rather than a compact one.⁴³ To further investigate the relation between the intrinsic oop/stretch ratio and PAH size, PAHs with more than 4 trio hydrogens have therefore been omitted in Fig. 1b. Moreover, to remove the bias towards the more numerous smaller PAHs an average has been taken for fixed sizes (tolerance 5 carbon atoms). Quite interestingly, we then find that the oop/stretch ratio and the number of carbon atoms in a PAH exhibit a clear drop with size that converges to an intrinsic oop/stretch ratio of 0.35 for large PAHs. This drop can be well represented by an exponential function given by

With R_oop/stretch the intrinsic oop/stretch intensity ratio and N_C the number of carbon atoms.

The large number of theoretical IR spectra in the PAHdb have provided extensive insight into the interstellar PAH population and spectroscopic properties of PAHs. However, considering the intensities, the harmonic, relatively low-level calculations are known to consistently overestimate the CH stretch band intensity.²⁹ To corroborate these findings, in Fig. 2 we compare the calculated ratio from the PAHdb with the ratio determined from gas-phase direct absorption measurements from the NIST database.³⁵ As discussed in the introduction, anharmonic calculations lead to a significant improvement in predicting IR spectra. The improvement mainly lies in the addition of combination bands and resonances that carry a significant amount of intensity, especially in the CH stretch region. Fig. 2 thus reports on a comparison in terms of intensities resulting from such calculations. Finally, a comparison is made between the gas-phase absorption measurements to cold, matrix isolation absorption measurements.

Fig. 2 Histogram of the oop/stretch intensity ratio determined using different theoretical and experimental methods divided by the gas-phase direct absorption measurements. The bin size is 0.15. The intensities were integrated between 650–950 cm⁻¹ (15.4–10.5 μm) and 2750–3250 cm⁻¹ (3.6–3.1 μm). The PAH species that are taken into consideration are listed in the ESI.† The mean and deviation of the distributions are indicated in the plot.

The oop/stretch intensity ratio determined using the theoretical spectra of the PAHdb is smaller than the experimental gas-phase ratio by a factor of 0.66. This corroborates the previous findings that the CH stretch intensities are overestimated, although this overestimation appears to be even slightly larger than reported previously.²⁷ Interestingly, the ratio is quite consistently off, which supports the approach to scale the band strengths computed at the B3LYP/4-31G level to obtain more accurate ratios. This approach was used before to obtain band strengths for methylated PAHs that were as accurate as those calculated at far more expensive levels (MP2/6-311+G(3df,3pd)).⁴⁴ We compare the DFT calculated ratios for the coronene family for different basis sets and find that upon increasing the basis set the oop/stretch intensity ratio increases (see Table S2†).

The ratios that result from anharmonic calculations show on average a remarkably good agreement with the gas-phase absorption experiments in terms of accuracy. It should be noted that a larger basis set (N07D, which is a modification of 6-31G*) is used in these anharmonic calculations than has been used for the PAHdb and such a larger basis set is expected to lead to a more accurate ratio, mainly resulting from a smaller stretch intensity.²⁹ Nevertheless, examination of the effects of the inclusion of anharmonicity on the ratio is of great interest as a significant amount of intensity is involved in the combination bands and resonances. The ratio determined from matrix isolation spectroscopy measurements shows the largest spread with respect to gas-phase absorption experiments. The comparison of the ratios of individual PAHs can be found in Fig. S1 in the ESI.†

3.2 Correcting for systematic deviation

In view of the above observations, it becomes quite interesting to determine how the underestimation of the oop/stretch intensity ratio affects conclusions on the PAH size distributions. To this purpose, we used the harmonically calculated intrinsic band strengths and applied a full temperature cascade emission model from the PAHdb at 6.5 eV on a large selection of PAHs (for details, see the Methods section). The results are shown in Fig. 3 (grey squares). In order to correct for the underestimation of the intrinsic oop/stretch intensity ratio by low-level calculations, a coarse modification has been made by scaling the emission ratio with the factor by which the low-level theory calculations underestimate the intrinsic ratio (1/0.66, light-green squares in Fig. 3). Such an approach was used before for methylated PAHs, where the band strengths of inexpensive calculations were scaled to match the more expensive predictions.⁴⁴ This figure also displays the results of a more accurate approach in which the intrinsic ratio is scaled and afterwards the full temperature cascade emission model is applied (dark-green squares). Aiming to facilitate a readout of these plots and serving as a guide to the eye, the data sets have been fitted with an exponential function.

Fig. 3 Grey squares indicate the oop/stretch intensity ratio as a function of the number of carbon atoms determined using the calculated intrinsic band strengths and a full temperature cascade emission model from the PAHdb on a large selection of PAHs (for details, see the Methods section). The green squares are scaled by a factor of 1/0.66 (see text) to correct for the systematic underestimation of the intrinsic oop/stretch intensity ratio by low-level calculations. Intensities were integrated between 830–985 cm⁻¹ (12.05–10.15 μm) and 2950–3115 cm⁻¹ (3.4–3.2 μm).

Since the average absorbed UV photon energy varies between astronomical objects, we have calculated the emission ratio as a function of molecular size for cascade starting energies of 6.5 and 8 eV, these energies being representative for the average photon energy of a reflection nebula such as NGC 7023 (ref. 12) and for a Photon Dominated Region such as the Orion Bar,²⁴ respectively. We assume one-photon absorption processes as cooling rates are high compared to absorption rates in the objects studied here.⁶ Multiphoton events would result in a decrease of the oop/stretch ratio. The results of these calculations are shown in Fig. 4 with the ratio of correcting both before and after applying the emission model plotted as transparent and opaque color, respectively. While correcting the intensity ratio before applying the emission model is most appropriate, we also include the results derived by applying the correction factor after applying the emission models as a convenient way to interpret the results of the current PAHdb. An exponential function to facilitate readout has been fitted to all plots and is provided in the ESI.†

Fig. 4 The oop/stretch intensity ratio as a function of the number of carbon atoms determined using the calculated band strengths and a full temperature cascade emission model from the PAHdb on a large selection of PAHs. The transparent squares correspond to the ratio scaled after applying the emission model, the opaque squares to scaling the intrinsic strengths before applying the emission model, both cases aiming to correct for the systematic underestimation of the intrinsic oop/stretch intensity ratio by low-level calculations. The solid and dotted black horizontal lines are examples of lower and higher limits of the observed intensity ratio in NGC 7023 and the Orion Bar, respectively.^12,24 The exponential fits indicated with the colored solid lines are phenomenological and meant to facilitate readout (see Table S1† for fit parameters).

With both the ‘original’ model and the corrected model, the 11.2/3.3 μm intensity ratio can now be translated into a PAH size distribution. The two horizontal solid black lines in Fig. 4 are typical lower and higher limits of observed intensity ratios in NGC 7023.¹² The two horizontal dotted black lines indicate lower and higher limits of the intensity ratios in the Orion Bar.²⁴ For NGC 7023, using the ‘original’ model and taking the indicated upper and lower observed ratio limits, one would conclude that typical PAH sizes are in the 50–70 carbon atoms range. Importantly, however, we find that when using the corrected model which takes the systematic underestimation of the intensity ratio into account, this range is reduced to 40–55 carbon atoms. For the Orion Bar region, the PAH size range is determined to be between 37–53 carbon atoms.

The results of the present model are quite remarkable but need to be substantiated further with more extensive analyses that incorporate accurate high-level anharmonic calculations on PAHs with a wide range of sizes in the emission model. Such analyses are now indeed high on our priority list. However, as for larger PAHs these calculations are not available yet, this model awaits completion. Nevertheless, our results do indicate that the underestimation of the oop/stretch intensity ratio by low-level calculations translates into a significant overestimation of the PAH size distribution. Equally important to notice is that our results are still in line with the hypothesis that an interconversion between large PAHs and the even more stable C₆₀ fullerene takes place at a certain level of the UV radiation field.⁷ From the higher limit of about 55 carbon atoms that our results imply, it appears that PAHs containing more than 60 carbon atoms should indeed be readily converted into fullerene species.

4. Conclusions

By comparing calculations from the NASA Ames PAH IR spectroscopic database (PAHdb) and state-of-the-art anharmonic calculations with gas-phase direct absorption IR spectra from the NIST chemistry webbook, this study validates the predicted intrinsic oop/stretch intensity ratio of PAHs. Currently, this ratio is recognized as the most reliable tracer for determining the size distribution of PAHs in the interstellar medium. Our studies show that the relatively low-level calculations that have been employed so far to interpret astronomical observations consistently underestimate the intrinsic intensity ratio by about 34%. Predictions made with a larger basis-set and taking anharmonic effects such as combination bands and resonances into account are, in contrast, in good agreement with the experiment. Unfortunately, the high-level anharmonic calculations are not available yet for large PAHs. Nevertheless, the consistency of the underestimation of the intensity ratio allows for making an ad hoc, coarse correction to the emission model based on the PAHdb. Such a corrected ratio–PAH size relation suggests that the size distribution, which up till now has been taken to be in the 50–70 carbon atoms range, could very well need to be adjusted to the 40–55 range. It is quite satisfying to notice that the higher limit of the PAH size conforms to the hypothesis that an interconversion takes place between larger sized PAHs and fullerenes in the interstellar medium.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We dedicate this article to Tim Lee who passed away November 3, 2022. His contributions have been key to advancing quantum chemistry in general and astrochemistry in particular. He will be sorely missed. We would like to thank Dr Christof Iserlohe for his computational support. This work is part of the Dutch Astrochemistry Network (DAN) II (project no. 648.000.029).

References

1. A. G. G. M. Tielens, Annu. Rev. Astron. Astrophys., 2008, 46, 289–337.
2. L. J. Allamandola, A. G. G. M. Tielens and J. R. Barker, Astrophys. J., Suppl. Ser., 1989, 71, 733–775.
3. A. Leger and J. L. Puget, Astron. Astrophys., 1984, 137, L5–L8.
4. E. Peeters, C. Mackie, A. Candian and A. G. G. M. Tielens, Acc. Chem. Res., 2021, 54, 1921–1933.
5. F. Galliano, S. C. Madden, A. G. G. M. Tielens, E. Peeters and A. Jones, Astrophys. J., 2008, 679, 310–345.
6. A. G. G. M. Tielens, Molecular Astrophysics, Cambridge University Press, Cambridge, 2021.
7. O. Berné and A. G. G. M. Tielens, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 401–406.
8. C. J. Mackie, T. Chen, A. Candian, T. J. Lee and A. G. G. M. Tielens, J. Chem. Phys., 2018, 149, 134302.
9. A. Maragkoudakis, E. Peeters and A. Ricca, Mon. Not. R. Astron. Soc., 2020, 494, 642–664.
10. W. A. Schutte, A. G. G. M. Tielens and L. J. Allamandola, Astrophys. J., 1993, 415, 397.
11. A. Ricca, C. W. Bauschlicher, C. Boersma, A. G. G. M. Tielens and L. J. Allamandola, Astrophys. J., 2012, 754, 75.
12. B. A. Croiset, A. Candian, O. Berné and A. G. G. M. Tielens, Astron. Astrophys., 2016, 590, A26.
13. M. Jourdain de Muizon, P. Cox and J. Lequeux, Astron. Astrophys., Suppl. Ser., 1990, 83, 337–355.
14. T. I. Mori, I. Sakon, T. Onaka, H. Kaneda, H. Umehata and R. Ohsawa, Astrophys. J., 2012, 744, 68.
15. C. Pech, C. Joblin and P. Boissel, Astron. Astrophys., 2002, 388, 639–651.
16. H. Andrews, C. Boersma, M. W. Werner, J. Livingston, L. J. Allamandola and A. G. G. M. Tielens, Astrophys. J., 2015, 807, 24.
17. C. W. Bauschlicher, E. Peeters and L. J. Allamandola, Astrophys. J., 2009, 697, 311–327.
18. S. E. Stein, J. N. Spencer, J. R. Seigart, E. Brown, R. L. Sensing and S. E. Stein, J. Phys. Chem., 1978, 82, 566–571.
19. E. Dwek, R. G. Arendt, D. J. Fixsen, T. J. Sodroski, N. Odegard, J. L. Weiland, W. T. Reach, M. G. Hauser, T. Kelsall, S. H. Moseley, R. F. Silverberg, R. A. Shafer, J. Ballester, D. Bazell and R. Isaacman, Astrophys. J., 1997, 475, 565–579.
20. C. W. Bauschlicher, C. Boersma, A. Ricca, A. L. Mattioda, J. Cami, E. Peeters, D. Armas, G. P. Saborido, D. M. Hudgins and L. J. Allamandola, Astrophys. J., Suppl. Ser., 2011, 189, 341–351.
21. C. W. Bauschlicher, A. Ricca, C. Boersma and L. J. Allamandola, Astrophys. J., Suppl. Ser., 2018, 234, 32.
22. C. Boersma, J. Bregman and L. J. Allamandola, Astrophys. J., 2015, 806, 121.
23. B. T. Draine, A. Li, B. S. Hensley, L. K. Hunt, K. Sandstrom and J.-D. T. Smith, Astrophys. J., 2021, 917, 3.
24. C. Knight, E. Peeters and A. G. G. M. Tielens, Mon. Not. R. Astron. Soc., 2022, 509, 3523–3546.
25. C. Knight, E. Peeters, D. J. Stock, W. D. Vacca and A. G. G. M. Tielens, Astrophys. J., 2021, 918, 8.
26. C. J. Mackie, A. Candian, T. J. Lee and A. G. G. M. Tielens, Theor. Chem. Acc., 2021, 140, 124.
27. S. R. Langhoff, J. Phys. Chem., 1996, 100, 2819–2841.
28. C. J. Mackie, A. Candian, X. Huang, E. Maltseva, A. Petrignani, J. Oomens, W. J. Buma, T. J. Lee and A. G. G. M. Tielens, J. Chem. Phys., 2015, 143, 224314.
29. H. Geindre, A. R. Allouche and D. Peláez, J. Comput. Chem., 2021, 42, 1018–1027.
30. C. J. Mackie, A. Candian, T. J. Lee and A. G. G. M. Tielens, J. Phys. Chem. A, 2022, 126, 3198–3209.
31. A. K. Lemmens, D. B. Rap, J. M. M. Thunnissen, C. J. Mackie, A. Candian, A. G. G. M. Tielens, A. M. Rijs and W. J. Buma, Astron. Astrophys., 2019, 628, A130.
32. C. Joblin, L. D'Hendecourt, A. Leger and D. Defourneau, Astron. Astrophys., 1994, 281, 923–936.
33. C. Joblin, P. Boissel, A. Léger, L. D'Hendecourt and D. Défourneau, Astron. Astrophys., 1995, 299, 835–846.
34. C. J. Mackie, The Anharmonic Infrared Spectra of Polycyclic Aromatic Hydrocarbons, 2018.
35. NIST Mass Spectrometry Data Center, W. E. Wallace, director, in NIST Chemistry WebBook, NIST Standard Reference Database Number 69, ed. P. J. Linstrom and W. G. Mallard, National Institute of Standards and Technology, Gaithersburg, MD, Infrared Spectra, accessed Jan 2022.
36. A. L. Mattioda, D. M. Hudgins, C. Boersma, A. Ricca, E. Peeters, J. Cami, F. Sánchez de Armas, G. Puerta Saborido, C. W. Bauschlicher and L. J. Allamandola, Astrophys. J., 2020, 251, 22.
37. V. Barone, J. Chem. Phys., 2005, 122, 014108.
38. V. Barone, M. Biczysko and J. Bloino, Phys. Chem. Chem. Phys., 2014, 16, 1759–1787.
39. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
40. V. Barone, P. Cimino and E. Stendardo, J. Chem. Theory Comput., 2008, 4, 751–764.
41. C. J. Mackie, A. Candian, X. Huang, E. Maltseva, A. Petrignani, J. Oomens, W. J. Buma, T. J. Lee and A. G. G. M. Tielens, Phys. Chem. Chem. Phys., 2018, 20, 1189–1197.
42. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16, Revision A.03, Gaussian, Inc., Wallingford CT, 2016.
43. J. Bouwman, P. Castellanos, M. Bulak, J. Terwisscha van Scheltinga, J. Cami, H. Linnartz and A. G. G. M. Tielens, Astron. Astrophys., 2019, 621, A80.
44. X. J. Yang, A. Li, R. Glaser and J. X. Zhong, Astrophys. J., 2017, 837, 171.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2fd00180b

---