Anharmonic Infrared Spectra of Thermally Excited Pyrene(C_{16}H_{10}): the combined view of DFT-based GVPT2 with AnharmonicCaOs and approximate DFT molecular dynamics with DemonNano
Shubhadip Chakraborty, Giacomo Mulas, Mathias Rapacioli, Christine Joblin
MModeling Anharmonic Infrared Spectra of Thermally Excited Pyrene(C H ): the combined view of DFT AnharmonicCaOs and approximateDFT molecular dynamics Shubhadip Chakraborty a , Giacomo Mulas a,b, ∗ , Mathias Rapacioli c , Christine Joblin a a Institut de Recherche en Astrophysique et Plan´etologie, Universit´e de Toulouse (UPS), CNRS, CNES, 9 Av. du ColonelRoche, 31028 Toulouse Cedex 4, France. b Istituto Nazionale di Astrofisica (INAF), Osservatorio Astronomico di Cagliari, 09047 Selargius (CA), Italy c Laboratoire de Chimie et Physique Quantiques (LCPQ/IRSAMC), Universit´e de Toulouse (UPS),CNRS, 118 Route deNarbonne, 31062 Toulouse, France
Abstract
Aromatic Infrared Bands (AIBs) are a set of bright and ubiquitous emission bands, observed in regionsilluminated by stellar ultraviolet photons, from our galaxy all the way out to cosmological distances. Theforthcoming James Webb Space Telescope will unveil unprecedented spatial and spectral details in the AIBspectrum; significant advancement is thus necessary now to model the infrared emission of polycyclic aro-matic hydrocarbons, their presumed carriers, with enough detail to exploit the information content of theAIBs. This requires including anharmonicity in such models, and to do so systematically for all speciesincluded, requiring a difficult compromise between accuracy and efficiency.We propose a new recipe using minimal assumptions on the general behaviour of band positions and widthswith temperature, which can be defined by a small number of empirical parameters. We explore here theperformances of a full quantum method, AnharmoniCaOs, relying on an ab initio potential, and MolecularDynamics simulations using a Density Functional based Tight Binding potential to determine these param-eters for the case of pyrene for which high temperature gas-phase data are available. The first one is veryaccurate and detailed, but it becomes computationally very expensive for increasing T; the second tradessome accuracy for speed, making it suitable to provide approximate, general trends at high temperatures.We propose to use, for each species and band, the best available empirical parameters for a fast, yet suffi-ciently accurate spectral model of PAH emission properly including anharmonicity. Modelling accuracy willdepend critically on these empirical parameters, allowing for an incremental improvement in model results,as better estimates become gradually available.
Keywords:
Anharmonic Infrared Spectroscopy, Polycyclic Aromatic Hydrocarbons, QuantumChemistry, Molecular dynamics, Astrochemistry
1. Introduction
The study of the Aromatic Infrared Bands (AIBs)in astronomical environment has opened interestingspectroscopic questions on the effect of anharmonic-ity on the infrared (IR) spectrum of hot polycyclicaromatic hydrocarbons (PAHs) and related speciesin an isolated environment. The AIBs are a familyof IR bands whose main components fall at ∼ ∗ Corresponding author µ m. They are ob-served in astronomical environments in which theircarriers can be excited by the absorption of ultra-violet (UV) photons. PAH species constitute thebest carriers since they have IR bands which corre-spond at least to first order to the AIBs and theycan reach high temperatures (1000 K or more) uponthe absorption of a single UV photon due to theirlow heat capacity. Because of the extreme isolationin space, the hot molecules can then relax by slowinfrared emission [1].The coming James Webb Space Telescope (JWST) Preprint submitted to Journal of Molecular Spectroscopy February 15, 2021 a r X i v : . [ a s t r o - ph . GA ] F e b ill provide a wealth of new data on the AIBs andtheir spatial variation, gathering a mine of infor-mation about the chemical identity of the emittingPAHs and their chemical evolution in various en-vironments. However access to this mine will beconditioned by our capacity to produce relevantsynthetic spectra with a model that can describethe photophysics of a given PAH molecule in agiven UV-visible astrophysical radiation field [2].For a comparison with astronomical spectra, such amodel needs to take into account anharmonic effectswhile simulating the IR emission spectra. Howeverthe effect of anharmonicity is usually included ina simplistic way using an average band shift (typ-ically ∼
10 cm − ) and broadening (typically ∼ − ) (cf. for instance the tool available in theAmesPAHdbIDLSuite in order to calculate an emis-sion spectrum from a theoretical IR spectrum asavailable in the NASA Ames PAH IR SpectroscopicDatabase [3, 4]). So far, only a couple of modelshave tried to go beyond these approximations whilesimulating the cooling of hot PAHs to match astro-nomical observations [5, 6, 7, 8]. The first threemodels were based on empirical parameters thatquantify the linear slopes derived from the evolutionof band positions with temperature in the spectra ofthermally excited PAHs [9]. The last model by Mu-las et al. [8] was based only on theory and focusedon the case of the small molecules, namely naphtha-lene (C H ) and anthracene (C H ), for whichthe matrix of anharmonic parameters of the Dun-ham expansion expressing vibrational energy couldbe calculated, and included in the Monte Carlo sim-ulation. With the coming JWST, it appears verytimely to consider these effects in a more system-atic way.In the laboratory, several techniques were used toquantify anharmonicity. The evolution of the IRspectra of neutral gas-phase PAHs at high temper-atures was studied using thermal excitation in gascells [10, 9]. In the case of PAH ions, other tech-niques have to be used and Oomens et al. [11] dis-cussed the potential of IR multiple photon dissocia-tion (IRMPD) action spectroscopy of trapped ions.Although IRMPD spectra contain non-linearitiesand cannot be considered as representative of IRabsorption or emission spectra, they carry infor-mation on the impact of anharmonicity on bandpositions and widths. The technique can be con-veniently used for large PAH ions of astrophysicalinterest [12] to provide at a first order IR spectraof PAHs containing an internal energy compara- ble to that of astro-PAHs absorbing UV photons.Experiments recording the IR emission spectra ofsmall UV-excited PAHs were performed in the early1990ies [13, 14, 15]. This work performed on the3.3 µ m band of small PAHs (up to pyrene in thecase of Shan et al. [14]) was extended to the fullmid-IR range thanks to the IR photon-counting ex-periment developed in the Saykally’s group [16, 17].These difficult experiments demonstrated the valid-ity of the PAH emission model in astrophysical envi-ronments. They illustrated that anharmonic effectsredshift band positions and lead to broadened andasymmetric band profiles. They constitute a veryvaluable dataset if one could precisely quantify theexcitation and de-excitation conditions (role of col-lisions?) in the experiments as discussed in Cooket al. [17].So far, the experimental measurements performedon thermally excited gas-phase PAHs have providedthe easiest way to quantify anharmonicity for astro-physical models [5, 7]. This is due to the fact thatthe evolution with temperature of the band posi-tions and widths was found to be linear in the stud-ied range ( ∼ ab ini-tio engine (e.g. DFT or DFTB) for on the fly eval-uation of the potential energy and electric dipole(and/or higher electric/magnetic multipole) hyper-surfaces. This method treats ionic motion clas-sically, thereby losing information on the discretestructure of individual quantised vibrational states.However, it does not make any a priori assumptionon selection rules, merely following the time evo-2ution of physical quantities and obtaining normalfrequencies, transition energies and intensities fromthe numerical autocorrelation function of the ap-propriate function. It therefore naturally obtainsalso combination/overtone bands, provided one ac-cumulates enough data through many, long enoughsimulations. Achieving ergodicity at relatively lowenergies can be challenging though, making thismethod better suited for T &
500 K. This methodwas used with some success in the context of theinvestigation of interstellar PAHs spectral proper-ties [22, 23, 24, 25, 26, 27].The second approach is that of a full quantum cal-culation, the most successful of which, for PAHs,was the Van Vleck formalism of 2nd order pertur-bation theory, with explicit variational treatmentof resonances, based on a 4th order Taylor expan-sion of the PES at a stable geometry and a 2nd or-der Taylor expansion of the electric dipole (and/orhigher electric/magnetic multipole) moment [28].The Taylor expansions are obtained via ab initio calculations, e. g. with DFT. This representation ofthe PES and multipole moments makes this methodbest suited for fairly rigid molecules and/or not toohigh vibrational excitation, so that these approxi-mations remain good for the vast majority of theaccessible phase space. This method is thus verywell suited for PAHs at not too high temperatures.The main disadvantage is that sampling adequatelyphase space with this method becomes computa-tionally very expensive at moderately high T, dueto the huge number of resonating states, resultingin huge effective Hamiltonians to be built and diag-onalised, and in a correspondingly huge number ofdipole matrix elements between their eigenvectorsto be computed. This method was successfully ap-plied to naphthalene [29, 30] and, more recently toslightly larger PAHs [see e. g. 31, 32, 33].We developed a quantum chemistry code, Anhar-moniCaOs, to study the detailed anharmonic spec-tra of PAHs with a fully quantum approach [34]. Afirst milestone for this code has been to model theeffects of anharmonicity at 0 K, including the effecton the position of the fundamental bands but alsothe contribution of the overtone, combination, anddifference bands to the spectra. The results of thecode were compared to available experimental dataon pyrene (C H ) and coronene (C H ). Whilenew experimental data in the gas-phase and at lowtemperatures are becoming available [35, 36, 31],this code can support the interpretation of suchexperimental data and, once calibrated and vali- dated against them, extend calculations to spectrawhich are more difficult to measure in the labora-tory, e.g. vibrational spectra of ions and/or rad-icals. In this article the focus is to evaluate theability of this code in modelling spectra at highertemperatures. We expect that at temperatures wellabove ∼
500 K, the system will spend an increas-ingly larger part of its time in regions of the po-tential energy surface far from the optimised ge-ometry, which cannot be anymore described by aquartic scheme. In order to access anharmonicityfactors in such cases, we will use the alternative ap-proach given by molecular dynamics, in which theIR spectrum is extracted from a Fourier transformof the dipole autocorrelation function. As millionsof energy and forces calculations are mandatory forsuch scheme, the potential energy will be computedfrom a parameterized method, namely the DensityFunctional based Tight Binding (DFTB), an ap-proximated DFT scheme with a much lower com-putational cost. This approach has already beenshown to provide reasonable qualitative and quan-titative trends of anharmonicity factors for PAHbased systems [24, 25].This article compiles the results which have beenobtained on the evolution of the IR band posi-tions and widths for pyrene as derived from twodifferent modelling approaches, the ab initio An-harmoniCaOs code for the low temperature range(up to ∼
500 K) and the DFTB-based MD simula-tions for the higher temperatures. The obtained re-sults are also compared with available experimentaldata. This is the first study of the kind dedicatedto anharmonic effects in the IR spectra of PAHs inwhich the performances of two modelling methodsare discussed in relation with known experimentaldata.
2. Methodology
To obtain anharmonic vibrational spectra ofPyrene at non-zero temperature, we made use ofthe AnharmoniCaOs code, developed by some ofthe authors as described in Mulas et al. [34]. Inbrief, AnharmoniCaOs takes in input the harmonicfrequencies of all normal modes, the third andfourth derivatives of the potential energy with re-spect to normal coordinates, as well as the firstand second derivatives of the electric dipole mo-ment, and obtains anharmonic states and electric3ipole-permitted transitions. It does so by ap-proximating the potential with its Taylor expan-sion around its minimum (the stable geometry ofthe molecule), truncated to fourth order, and theelectric dipole with its Taylor expansion around thesame point, truncated to second order. These Tay-lor expansions for pyrene were taken from Mulaset al. [34]. AnharmoniCaOs then proceeds with amixed variational-perturbative approach using theVan Vleck formalism: resonant terms are consid-ered explicitly, building effective Hamiltonian ma-trices and numerically solving them to obtain an-harmonic eigenstates, while non-resonant terms aredealt with by perturbation theory. Anharmoniceigenstates are represented as linear combinationsof harmonic states. Matrix elements of the elec-tric dipole operator (approximated by its 2 nd orderTaylor expansion) are then obtained between theanharmonic eigenstates obtained from the diago-nalization of the effective Hamiltonians.The set of harmonic states included in each calcula-tion is built around a given harmonic “start state”,which results in the formation of truncated polyadsof states connected by resonances. For more de-tails on the algorithms used to build and truncatepolyads, see Mulas et al. [34]. As a result of the waypolyads are built and truncated, one finds the an-harmonic states with a significant component alongthe initial “starting state”. States with the largestprojection along the “starting state” are numeri-cally most accurate, while accuracy degrades as an-harmonic states get farther and farther “away” (interms of projection) from the “starting state”. Foreach of these calculations, therefore, only the an-harmonic states closest to the “starting state” areconsidered in building the spectra, with a weightproportional to the square module of their projec-tion along the “starting state”.If one were to repeat the calculation for all possible“starting states”, one would obtain all anharmonicstates and transitions. However, the density of har-monic states, for a polyatomic molecule, is an ex-tremely steep function of vibrational energy, mak-ing complete sampling of harmonic states unfeasibleabove a relatively low energy. We therefore enumer-ated harmonic states, in order of harmonic energy,till we reached a predetermined maximum number.From that harmonic energy on ( ∼ − ), weresorted to statistical sampling. We used the Wang-Landau method [37, 38] to perform a random walkon harmonic states, obtaining a nearly uniformsampling in harmonic energy space of the “start- ing states”. This does not result in a completelyuniform sampling in the anharmonic energy spacefor the states obtained by AnharmoniCaOs; still,this produces a very slowly varying sampling den-sity in anharmonic energy space, ensuring that allenergy ranges randomly explored are evenly sam-pled. This sampling method was successfully used,e. g., by Basire et al. [30] to estimate the anhar-monic density of states of naphthalene.The collection of these AnharmoniCaOs runs pro-duces anharmonic vibrational spectra as a functionof vibrational energy: the vibrational energy is dis-cretised in a number of bins, each covering an in-terval of 5 cm − ; individual spectra consisting oftransitions from an anharmonic state in that in-terval are averaged, to estimate the spectrum ofmolecules in thermal equilibrium at the correspond-ing vibrational energy, in the microcanonical en-semble; these spectra can then undergo a numeri-cal Laplace transform to finally obtain anharmonicvibrational spectra as a function of temperature inthe canonical ensemble, directly comparable withlaboratory experiments. We have also computed the anharmonic spec-tra for pyrene from a dynamical exploration of thePES at various temperatures. Such dynamics re-quires the computation of millions of single pointenergies and gradients and could hardly be per-formed at the ab initio level. A good compromiserelies on the use of Density Functional based TightBinding (DFTB)[39, 40], which presents a muchcheaper computational cost than DFT while pre-serving a quantum description of the electronic sys-tem and allowing for the calculation of the molec-ular dipole ~µ ( t ) on the fly. In practice, we usedthe DFTB scheme in its second order formulation[41] implemented in deMonNano [42] with the MATset of parameters [43] combined with a correctionfor atomic charges [44] that already showed its ef-ficiency in the context of IR spectra calculations[24, 23, 26, 27]. The IR spectra were obtained fromthe Fourier transform of the auto-correlation dipolemoment using the formulation presented by [45] : α ( ω ) ∝ ω Z + ∞−∞ d t h ~µ (0) · ~µ ( t ) i e iωt (1)where <> indicates a statistical average to removedependency on the initial conditions.4s the equation is intended to be used to com-pute the spectum at a given temperature, the sim-ulations should be done in the canonical ensemble,which in practice means using a thermostat. Thelatter might however introduce errors as the pro-cess of energy exchange between the system and thethermostat involves frequencies which might pol-lute the computed spectrum. To reduce such arti-facts, we first perform, for each of the 6 investigatedtemperatures between 600 and 1600 K, a canoni-cal dynamics of 50 ps using a Nose Hoover chainof thermostat [46, 47, 48] (5 thermostats with en-ergy exchange frequency of 800 cm − ). About 50snapshots are selected along this dynamics and thecorresponding geometries and velocities are used togenerate the starting conditions for microcanoni-cal simulations of 5ps each with timestep of 0.1 fs.The IR spectrum is computed and averaged fromthese last simulations without thermostats, the ini-tial conditions sampling however the canonical dis-tribution at a given temperature. One of the aims of this work was to obtain,from more or less complete theoretical calculations,the evolution of band positions and bandwidthsas a function of temperature, for use in simplifiedmodelling of AIB emission in astronomical environ-ments, comparing results from different theoreticalmethods and experimental ones. This poses a fun-damental problem: how to define in an unambigu-ous way band positions and widths, reducing to aminimum their dependence on subjective choices ofthe scientist deriving them from the data. This isparticularly difficult for bands that result from theblending of several components, which sometimescan be resolved only at low T. We therefore decidedto define the centroid ¯ ν i of the band i as the math-ematical average of its position, weighted by bandintensity I ( ν ), over a fiducial interval [ ν min i , ν max i ]defined by the scientist for that band, i. e.¯ ν i = R ν max i ν min i d ν ν I ( ν ) R ν max i ν min i d ν I ( ν ) = 1 I i Z ν max i ν min i d ν ν I ( ν ) , (2)where I i is the integrated intensity of band i . Simi-larly, we define the bandwidth σ i as the square rootof σ i = 1 I i Z ν max i ν min i d ν ( ν − ¯ ν i ) I ( ν ) . (3) These definitions are reminiscent of the definitionof the first two central momenta of a statistical dis-tribution function. They are identically applicableto DFT-AnharmoniCaOs spectra, DFTB/MD, andexperimental ones, making it possible to meaning-fully compare positions and widths, thus defined,obtained from all three sets of data. The intervals[ ν min i , ν max i ] for each band were defined by visualinspection of the spectra, and are the only remain-ing subjective choice. The intervals we adopted aregiven in Tables 2 and 3 in the supplementary ma-terial.
3. Results
The raw spectra resulting from AnharmoniCaOsare “stick” spectra, each “stick” resulting from anindividual anharmonic transition and arbitrarily (inprinciple infinitely) sharp, since AnharmoniCaOsdoes not include rotational structure. For the re-sults presented here, the resolution of the calcu-lations yields a “stick” width of 0.2 cm − . Fig-ure 1 gives an overall comparison with available ex-perimental spectra, both in solid state and in gasphase. For viewing convenience, and easier com-parison with experimental data, in which bands arebroadened either by unresolved rotational structure(in gas phase) or by solid state effects, we con-volved AnharmoniCaOs spectra with Gaussian pro-files with a FWHM of ∼ − . The theoreticalspectra clearly reproduces fairly closely the exper-imental ones, not only in terms of fundamentalsbut also including combination and overtone bands,permitting their individual identification with fairlyhigh accuracy.Figure 2 zooms in on the spectral region around theband at ∼ − and shows the Anharmoni-CaOs spectrum computed at 14 K and 523 K. Be-side the fundamental band near 1583 cm − , evenat low temperature one can clearly see many ad-ditional combination bands. Relevant individualbands in this spectral range were identified in Mulaset al. [34]. They can be rather accurately predictedthanks to the inclusion in AnharmoniCaOs of thesecond derivatives of the dipole moment, and tothe detailed treatment of resonances. On the otherhand, this detailed treatment of resonances, whichin case of several terms of the dipole moment ex-pansion contributing to the same transition causesthe code to merge several polyads in a single, larger5ne, renders the calculation increasingly heavy athigher excitation energies, making it computation-ally expensive, while still feasible, at very high tem-peratures. With rising temperature, a huge numberof additional hot bands progressively appear, cor-responding to transitions arising from increasinglyexcited states. These cluster mostly around the in-dividual bands already present at low T, with vari-able shifts around them, resulting in a clear broad-ening and, in some cases, shift of band centroids.Bands which are close at low T tend to blend athigh T, making it difficult, or even impossible, to re-solve them beyond some temperature. In addition,a multitude of very weak, resonance-activated fea-tures appear, with increasing T, even farther awayfrom the main bands, “borrowing” intensity fromthem. This results, qualitatively, in main bandslosing a fraction of their intensity at high T, whichgoes into a featureless, broad plateau below them.In theoretical spectra, this intensity is still presentin the spectrum, spread over a wider spectral range;in experimental ones, even if present (and there isno reason why it should not be), it is further spreadby rotational structure or solid-state effects, it thusmerges with the continuum, cannot be clearly dis-tinguished from it, and is thus easily partially orcompletely subtracted away with it, being lost tomeasurement.Figure 3 zooms in on the band at ∼
844 cm − , foran example of a spectral region with a well-resolved,single fundamental band. With the exception of theblending of initially resolved features, the same be-haviour described above is observed, with hot bandsappearing all around the fundamental, producing ashift of the centroid of the resulting band, as well asa broadening. In this case, as well, we also observethe appearance, at the highest T, of a “carpet” ofvery weak bands, “stealing” some of the intensityof the fundamental and spreading it over a broadplateau.Figure 4 gives another example to illustrate the re-markably large variability of bandwidths with tem-perature. In both cases, one band remains rela-tively narrow with increasing T, while neighbouringones show a much larger broadening.Figure 5 shows the evolution of the band profile forsome out of plane CH bending modes, as predictedby AnharmoniCaOs for T ranging from 14 K to523 K, compared with the gas-phase spectrum atthe closest temperature, i. e. 573 K. As for Fig-ure 1, the AnharmoniCaOs spectra were convolvedwith a Gaussian profile with ∼ − FWHM, to ease the comparison. The band position pre-dicted by AnharmoniCaOs very clearly tend to theones in the gas-phase spectrum, which is broadenedby rotational structure non included in Anharmon-iCaOs. The peak intensity predicted for the bandat ∼
741 cm − clearly decreases with temperature,as it is distributed over many hot bands around itand many more, weakly resonant features furtheraway.Figure 6 shows some more bands simulated byAnharmoniCaOs for the same temperature range.Many of them present significant resolved structureat low temperature, which is progressively lost as Tincreases. All of them show the general behaviourwe described above.Analysing the AnharmoniCaOs spectra as de-scribed in Sect. 2.3, we identified bands, or bandgroups, which could be followed for the whole tem-perature range spanned by AnharmoniCaOs cal-culations, and determined their T-dependent cen-troids and widths. Figure 7 reports, as black dots,the band positions for some bands, the others aregiven in the Supplementary Information. Similarly,Figure 8 shows the bandwidths for the same bands,the others being given in the Supplementary Infor-mation. In all cases, one can see that both positionsand widths follow different trends in different tem-perature ranges. This is in agreement with whatwas observed in laboratory experiments [18]. Weobtained the linear slopes by fitting data points intemperature ranges in which the behaviour is ap-proximately linear, which are reported in Table 3. The DFTB spectrum computed at the harmoniclevel at 0 K is a stick spectrum whereas the spec-tra obtained at finite temperature from DFTB-MDsimulations present broadenings, shifts and merg-ing of bands. The latters appear naturally, incor-porating the anharmonic effects of the full DFTBPES without truncation of the latter at a given or-der. Thermal effects result from MD explorations,and not from an a posteriori convolution of thesespectra with gaussian or lorentzian functions. Oneshould however remember that the frequency reso-lution scales as the reciprocal of the total durationof the MD segments used to numerically evaluateequation 1. In this work the frequency resolution is1.7 cm − .As an example, the spectrum computed at 600 Kis shown on Figure 9 together with the highest6 AnharmoniCaOsSolid-phase
Wavenumber (cm ) R e l a t i v e i n t e n s i t y AnharmoniCaOsGas-phase
Wavenumber (cm ) R e l a t i v e i n t e n s i t y (a) (b) Figure 1: (a)Theoretical infrared spectrum of pyrene from 3200 to 650 cm − at 300 K, simulated using DFT-AnharmoniCaOstogether with the experimental [18] condensed phase spectrum at 300K, (b)Theoretical infrared spectrum of pyrene from 3200to 650 cm − at 523 K, simulated using DFT-AnharmoniCaOs together with the experimental [9] gas phase spectrum at 573 KFigure 2: Stick spectra of pyrene from 1620 to 1530 cm − at 14, 300, 423 and 523 K computed using DFT-AnharmoniCaOs. - 1 )0100200300400500600 I n t e n s i t y ( k m m o l - c m ) - 1 )020406080 I n t e n s i t y ( k m m o l - c m ) - 1 )0102030405060 I n t e n s i t y ( k m m o l - c m ) - 1 )0510152025 I n t e n s i t y ( k m m o l - c m ) T = 14K T = 300KT = 423K T = 523K
Figure 3: Stick spectra of pyrene from 870 to 820 cm − at 14, 300, 423 and 523 K computed using DFT-AnharmoniCaOs. - 1 )051015 I n t e n s i t y ( k m m o l - c m ) - 1 )01234 I n t e n s i t y ( k m m o l - c m ) - 1 )0.00.51.01.52.0 I n t e n s i t y ( k m m o l - c m ) - 1 )0.00.51.01.5 I n t e n s i t y ( k m m o l - c m ) T =14K T = 300KT = 423K T = 523K
Figure 4: Stick spectra of pyrene from 470 to 560 cm − at 14, 300, 423 and 523 K computed using DFT-AnharmoniCaOs. µ m) Integration range (cm − ) Intensities (km · mol − )Gas 570 K[49] AnharmoniCaOs0 K[34] 300 K 523 K3.3 [2850-3250] 140( ±
10) 99 78.0 62.05.2 [1900-1950] 9.5 ( ± ± ± ± ± ± ± ± ± ± ± ± ±
6) 101 103.6 77.913.4 [730-750] 20.8 ( ± ±
1) 45.3 40.7 36.718.5 [526-550] 2.3 2.3 1.420.0 [487-501] 2.7 3.0 1.620.5 [474-487] 1.9 1.5 1.929.0 [335-354] 1.4 1.6 1.550.0 [185-212] 9.7 10.9 5.7
Table 1: Integrated intensities of the bands computed using AnharmoniCaOs at 0, 300 and 523 K together with the gas-phaseexperimental band intensities Wavenumber (cm ) I n t e n s i t y Figure 5: Temperature evolution of the band profiles for two bands of pyrene at 741.1 and 711.6 cm − simulated using DFT-AnharmoniCaOs together with gas-phase data at 573 K[49]. Wavenumber (cm ) R e l a t i v e i n t e n s i t y Wavenumber (cm ) R e l a t i v e i n t e n s i t y Wavenumber (cm ) R e l a t i v e i n t e n s i t y Wavenumber (cm ) R e l a t i v e i n t e n s i t y Wavenumber (cm ) R e l a t i v e i n t e n s i t y Wavenumber (cm ) R e l a t i v e i n t e n s i t y Figure 6: Temperature evolution of the simulated band profiles for a selection of bands of pyrene. The spectra were simulatedusing DFT-AnharmoniCaOs − ; two well defined bands below1000 cm − and a forest of bands between 1000 and2000 cm − . The main differences rely on the bandintensities, in particular for the DFTB bands at3000 and 720 cm − which have, respectively, higherand lower intensities when compared to their equiv-alent in the AnharmoniCaOs model. We remarkhere that, with respect to bands potitions, rela-tive bands intensities are much more challengingto reproduce by an approximated model and thisremains true for ab initio schemes. As a conse-quence, it is more difficult to reproduce the spec-tral evolution of a specific band located in a regionof high vibrational states density. In such a region,when the temperature increases, many bands willblend in a single, unresolved spectral feature. Thelarge errors, due to DFTB, in the estimated intensi-ties of individual, unresolved components will thenseverely affect the predicted evolution of the profileand position of the resulting blended feature.DFT and DFTB spectra computed at the har-monic levels allowed us to identify the nature ofthe bands on the basis of vibrational modes visu-alisation, the DFT and DFTB normal modes be-ing similar. Table 1 in Supplementary informationlists all harmonic normal modes resulting from theDFT vs DFTB calculations, matching all of themafter visual inspection of the corresponding ionicmotions. We will focus on the evolution of band po-sition and width for four relevant bands : bands 1and 11 which are the two most intense ones; band14 which is a satellite of band 11 whose evolutionwith temperature allowed for unambiguous deriva-tion of band position and width and, finally, band 2as representative of a more complex spectral region.The band positions for these 4 bands at 0 K (har-monic) and 600 K are reported in table 2. Whencompared to experimental values and keeping inmind the level of theory, the MD-DFTB band po-sitions at 600 K are in relatively good agreementwith experimental results for bands 1, 3 and 14,presenting relative errors of 5.1 , 3.1 and 2.6 %, alarger error of 15 % being observed for band 2. Itcan also be seen that the differences between DFTand DFTB band positions for temperatures around600 K are already present at the harmonic level:bands 1, 11 and 13 being redshifted in the DFTBspectra with respect to DFT spectra, the oppositebeing observed for band 2. Figure 10 shows the evo- lution with temperature of the most intense band,namely band 11, at the DFT-AnharmoniCaOs andMD-DFTB levels. It shows that, despite the factthat the absolute band positions differ between thetwo models, the trend, i.e. a redshift and broad-ening when increasing the temperature, seems tobe consistent, encouraging the use of the DFTB-MD model to compute the bands evolution at hightemperatures rather than their absolute positions.
4. Discussion
From the data we showed in Section 3, andtheir comparison with available experimental data,we can draw some general considerations. Inthe temperature interval explored by our DFT-AnharmoniCaOs calculations, we cannot give a“rule of the thumb” as to average band shift andbroadening as a function of temperature, valid forall bands in a given energy range: there are e. g.some notable cases of bands that remain muchsharper (more than an order of magnitude) thanneighbouring ones. One such example is apparentin Figure 4. Clearly, in these cases, none of the vi-brational modes that can be populated at ∼
523 Khave an important anharmonic coupling with thatband. Our current calculations cannot rule out thatat much higher temperatures these bands may showa more significant shift and/or broadening. Check-ing this, while possible, would require pushing ourAnharmoniCaOs calculations much further, imply-ing a quite substantial computational cost, way be-yond the scope of the present work.Table 1 shows the evolution of band intensities, asobtained from DFT-AnharmoniCaOs. We can seethat a bit less than half of the reported bands ap-pear to lose intensity with increasing temperature,in the calculated spectra. This is due to the effect,that we described in Section 3, of some of the inten-sity of main bands being “spread” over an increas-ingly large number of very weak bands, forming abroad plateau below them which, at 523 K, mergesin a sort of featureless continuum. The weakest ofthese hot bands are not even recorded by the code,falling below the intensity threshold to save them.This artifact results in the total integrated intensityover the whole vibrational spectrum to decrease byalmost 20% between 14 K and 523 K.The behaviour of band position and width (Fig-ures 7 and 8) appears not linear at low tempera-tures in the spectra obtained via a fully quantum13nharmonic calculation, in agreement with avail-able laboratory data [18]. This can be rationalisedby considering that, at any given temperature T,vibrational states whose frequency is much largerthan kT /h are not populated, and thus do not par-ticipate at all in the structure of hot bands collec-tively producing band broadening and shift. If avibrational mode has significant anharmonic cou-pling with a given band, it will produce a ratherabrupt change in the evolution of the band whenthe temperature becomes sufficient to populate it.In a temperature range in which all modes withlarge anharmonic coupling with a band are alreadypopulated, and no new important modes come intoplay, band position and width instead evolve lin-early with temperature. This accounts for the ob-served behaviour in Figures 7, 8, and figures in theSupplementary Information, which appears to re-semble broken lines. It also explains why, at highenough temperatures, all trends appear linear, asobserved in experimental data [9, 18].Classical molecular dynamics cannot account forsome modes being completely depopulated at lowtemperature, since they do not consider quantisedstates and include vibrations of arbitrarily smallamplitude, forbidden by quantum mechanics. Asa result, these simulations always predict a linearbehaviour of band positions and widths, regard-less of temperature. This behaviour is expected toapproach the quantum one when vibrational ener-gies are large enough. Figures 7, 8, and 2 and 4in the Supplementary Information show band po-sitions and widths obtained from DFTB-MD re-sults. As already mentioned in Section 3, abso-lute band positions from DFTB-MD are much lessaccurate compared to experiment. Still, their vari-ation with temperature seems qualitatively consis-tent with DFT-AnharmoniCaOs results, when theycan be compared. Some bands, in DFTB-MD spec-tra, become so shallow to be hardly distinguish-able from the continuum, making the measurementof their position and width unreliable. These areomitted our analysis.
5. The combined view for astrophysical ap-plications
Addressing the vibrational spectral evolutionover a wide temperature range, extending up to ∼ evolution with temperature of band posi-tions and widths, at very high T values, regardlessof their absolute positions and intensities.Figure 10 shows a zoom-in on the out of plane CHbend mode, as computed by DFT-AnharmoniCaOsfrom 14 K to 523 K, and by DFTB-MD from 600 Kto 1600 K. One can clearly see that while the abso-lute position is shifted, the slope of the trend givenby DFTB-MD seems to follow smoothly the one ofDFT-AnharmoniCaOs. Figure 7 reports togetherthe DFT-AnharmoniCaOs and DFTB-MD posi-tions for four specific bands. DFT-AnharmoniCaOs14nd DFTB-MD positions for all other bands aregiven in Figures 1 and 2 in the Supplementary Infor-mation. DFTB-MD positions in Figure 7 are rigidlyshifted to smoothly join DFT-AnharmoniCaOs val-ues at 523 K. Three panels in the same figurealso report data from available laboratory measure-ments [9]. Table 3 reports the slopes χ of the evolu-tion of band positions vs. temperature for T rangesin which they behave linearly. At frequencies higherthan 1300 cm − , the values of χ from DFTB-MDappear to be systematically a factor ∼ − . On the opposite, except for a couple ofbands, DFT-AnharmoniCaOs somewhat underesti-mates the χ values.A similar behaviour is observed when consideringbandwidths. Figure 8 shows bandwidths deter-mined from DFT-AnharmoniCaOs and DFTB-MD,including laboratory measurements for three bands[9]. Since the latter were performed in gas phasein thermal equilibrium, they also include a signifi-cant contribution from rotational structure, whichis included neither in DFT-AnharmoniCaOs nor inDFTB-MD modelling. We therefore estimated therotational contribution to bandwidth, computingrotational profiles for pyrene at the temperaturesof the laboratory data points. For our purposes,we used a simple rigid rotor model, with rotationalconstants taken from Malloci et al. [50]. Indeed,from high resolution vibrational spectra, with re-solved rotational structure, it appears that varia-tions of rotational constants of pyrene with vibra-tional quantum numbers, as well as due to cen-trifugal distortion, are negligible [51]. To estimatethe rotational width, we applied Equations 2 and 3given in Section 2.3 to the simulated rotational en-velopes of a-, b-, and c-type pyrene bands (pyrene isan asymmetric rotor). The rotational contributionto bandwidth is estimated to range from ∼ − at 573 K to ∼ − at 873 K for band 11 (anout of plane vibrational mode, fully c-type), andfrom ∼ − at 573 K to ∼ − at 873 Kfor band 1 (in plane C-H stretches, with similarcontribution from a- and b-type bands). To esti-mate the contribution of vibrational anharmonicityto the bandwidths measured in gas phase, we fol-lowed the approach of Pech et al. [7], subtractingthe estimated rotational contribution. The result-ing “corrected” laboratory data points are reportedin Figure 8. The absolute agreement of band-widths computed by DFT-AnharmoniCaOs with corrected laboratory data points is surprising. Inthe case of band 11, the slope of the experimentalpoints is markedly larger than the slope of DFT-AnharmoniCaOs points at lower temperature, hint-ing that it may change just around ∼
550 K. Acloser examination of Figure 3, indeed, shows sev-eral strong hot bands beginning to appear aroundthe main band at around ∼
500 K. Their increasedrelative contribution to the band is thus likely tocause a jump in the bandwidth. No shift was ap-plied to DFTB-MD data points in Figure 8. Thecomparison in this figure shows that DFTB-MDoverestimates bandwidths at the highest frequen-cies. Table 4 shows the slopes χ of bandwidths,similarly to Table 3 for band positions. The valuesof χ deduced from DFTB-MD results, even if qual-itatively consistent with experiment, appear less ac-curate than DFT-AnharmoniCaOs in the very fewcases in which comparison is possible. More exper-imental data would be needed to draw firm conclu-sions on this.The overall picture is consistent with DFTB-MDbehaving worse with unresolved band clusters pos-sibly involving resonances and better with isolatedbands, which is the case for pyrene at low frequen-cies. A tentative explanation is that when severalunresolved transitions contribute to an unresolvedband, the large intensity errors due to DFTB re-flects in a large error in the average position result-ing from their blend. Also, if the individual transi-tions involve states with the same symmetry, theyundergo resonance, and this depends critically onthe relative positions of the resonating states. Thisresults in large errors on the predicted effect of theresonance.In summary, both laboratory data and DFT-AnharmoniCaOs simulations find that band posi-tions are fairly stable below ∼
300 K, and have anapproximately linear asymptotic behaviour at hightemperatures. As to bandwidths, they also have anapproximately linear asymptotic behaviour at hightemperatures, and some of them are relatively sta-ble, like positions, below ∼
300 K, while others keepdecreasing down to very low temperatures, but witha different slope. There are no “one size fits all” val-ues that can be used to predict anharmonic bandpositions and widths for all bands, the behaviour ofindividual bands must be investigated, by theoreti-cal and/or laboratory experiments. Current modelsusing a fixed shift with respect to 0 K calculations,and/or fixed “standard” band widths to computethe emission of “astronomical” PAHs do not seem15herefore appropriate. Their inaccuracy is masked,when comparing with astronomical observations, bythe unresolved contribution of many different (un-known) individual molecules to astronomical AIBs,that does not allow one to separate the contributionof anharmonicity from those of chemical diversityand rotational broadening. On the other hand, theinaccuracy becomes apparent if one applies mod-els to cases in which laboratory data of individualmolecules are available, as we did here.Mackie et al. [52], proposed a simplified scheme forthe calculation of simulated PAH emission spectrafor comparison with AIBs, based on detailed the-oretical modelling of cascade spectra that hintedthat most bands did not show apparent shifts atdifferent excitation energies, and suggested that inmost cases only variations in bandwidths neededto be included for astronomical purposes. How-ever, this conclusion appears difficult to reconcilewith the linear shift of all band positions at hightemperatures observed in laboratory spectra, andreproduced by theoretical modelling in the presentwork.Based on our modelling, and on available labo-ratory data, we suggest the following simplifiedrecipe for fast, efficient, but still sufficiently accu-rate astronomical modelling of PAH emission, inthe framework of the thermal approximation: • below 300 K assume band positions tobe approximately constant; interpolate lin-early bandwidths between reference values at ∼
300 K and very low temperature (like our14 K data points); use as a reference the values(given in Tables 3 and 4 for pyrene) obtainedfrom either laboratory data or detailed simu-lation at these temperatures; • above 300 K assume band positions and widthsevolve following a linear trend, whose slopes χ and χ can be obtained by the best avail-able data: experimental, detailed modelling(such as e. g. DFT-AnharmoniCaOs here orSPECTRO [32]), molecular dynamics basedon full DFT or better, molecular dynamicsbased on faster approximate methods (such asDFTB used here); this choice should be an ed-ucated guess on a case-by-case basis, since e. g.DFTB-MD appears to provide fairly good re-sults for strong, isolated bands, and converselybad ones where band clusters merge; • estimate the statistical thermal distribution of the given PAH as a result of the radiationfield it is subjected to, and obtain the resultingspectrum as an appropriate weighted averageof spectra at each temperature estimated fol-lowing the two above points.We are in the process of developing such a model, inview of the forthcoming data from the James WebbSpace Telescope.
6. Conclusions
The difficulty to directly measure in laboratoryexperiments the vibrational emission of isolatedPAHs excited by the absorption of UV photons,in order to study astronomical AIBs, has led todevelop models to simulate the IR cooling cascadesof PAHs including anharmonic effects [5, 6, 7].Such modelling must be as simple as possible, to beapplicable in a systematic way, yet include enoughdetails to achieve sufficient accuracy. These details,the inputs of the models, must be as thoroughlyvalidated against experiments as possible, testingthem against all available laboratory data.Now, with the James Webb Space Telescopeexpected to be coming online in the near future,the concept of “sufficient accuracy” related to PAHemission models needs to be revisited, especially ifwe want to extract the maximum information fromthe unprecedented level of spectral and spatialdetails it will make available. This is the right timeto extend models to really include anharmonicityin a systematic way. Mackie et al. [52] showedthat they can perform detailed calculations, usinga method similar to our AnharmoniCaOs. How-ever, as with AnharmoniCaOs, fully modellinganharmonic PAH emission is a computationallyexpensive undertaking for every single molecule,which does not make it easily applicable in asystematic way to existing PAH databases [50, 3],in view of using them to interpret the detailed,subtle spectral structures that JWST data willunveil. On the other hand, cruder approximations,such as disregarding the dependence of bandpositions on excitation, even at the very highvibrational temperatures reached by an isolatedPAH after the absorption of an UV photon, appearincompatible with already available laboratorydata, as we have shown here. We thus propose herea recipe for an intermediate level of detail, whichwould be computationally less expensive than anall-out full-quantum anharmonic calculation, and16
500 1000 1500 2000
Temperature (K) B a n d p o s i t i o n ( c m ) Band1
DFT/ACDFTB/MDGas-phase
Temperature (K) B a n d p o s i t i o n ( c m ) Band2
DFT/ACDFTB/MD
Temperature (K) B a n d p o s i t i o n ( c m ) Band11
DFT/ACDFTB/MDGas-phase
Temperature (K) B a n d p o s i t i o n ( c m ) Band13
DFT/ACDFTB/MDGas-phase
Figure 7: Evolution of band positions with temperature obtained from DFT-AnharmoniCaOs (filled circles), DFTB-MD (emptycircles) calculations and gas-phase experimental from Joblin et al. [9] (black stars). DFTB-MD data points have been rigidlyshifted to smoothly join DFT-AnharmoniCaOs ones near ∼
550 K. These shifts are respectively +139, -230, +34 and +17.5 cm − for band 1, band 2, band 11, and band 13. Bands are labelled following Table 2. For band 13 the gas-phase band positionsoverlaps with DFTB-MD band positions therefore we have shifted all gas-phase data by +5 cm − for clarity.Table 2: Calculated band positions of four representative IR bands of pyrene. Experimental values are listed from the work ofChakraborty et al. [18] for the low temperature part and Joblin et al. [9] for the high temperature measurementsTheory ExperimentDFT DFTBHarmonic AnharmoniCaOs Harmonic MD Solid[18] Gas[49]Temp 0 K 0 K 300 K 523 K 0 K 600 K 300 K 570 KBand1 3164.6 3040.0 3040.9 3040.5 2960.3 2897.3 - 3052.0Band2 1635.4 1600.0 1585.3 1585.8 1856.3 1842.2 - 1597.0Band11 860.5 843.0 844.4 843.9 818.4 814.1 840.0 840.0Band13 755.6 739.0 741.1 739.7 725.1 721.7 749.4 740.0
500 1000 1500 2000
Temperature (K) B a n d w i d t h ( c m ) Band1
DFT/ACDFTB/MDGas-phase without rotaional widthGas-phase with rotational width
Temperature (K) B a n d w i d t h ( c m ) Band2
DFT/ACDFTB/MD
Temperature (K) B a n d w i d t h ( c m ) Band11
DFT/ACDFTB/MDGas-phase without rotaional widthGas-phase with rotational width
Temperature (K) B a n d w i d t h ( c m ) Band13
DFT/ACDFTB/MD
Figure 8: Evolution of bandwidths with temperature obtained from DFT-AnharmoniCaOs (filled circles), DFTB-MD (emptycircles) calculations. Bandwidths from laboratory data [9] are shown as they are (filled squares), and after subtraction of theestimated contribution due to molecular rotation in thermal equilibrium (black stars, see text for details). Bands are labelledfollowing Table 2
DFT-AnhCaOsDFTB-MD
Wavenumber (cm ) R e l a t i v e i n t e n s i t y Figure 9: Theoretical infrared spectrum of pyrene from 3200 to 650 cm − at 523 and 600 K, simulated using DFT-AnharmoniCaOs and DFTB-MD Wavenumber (cm ) R e l a t i v e i n t e n s i t y DFT-AnharmoniCaOs
Wavenumber (cm ) R e l a t i v e i n t e n s i t y DFTB-MD
Figure 10: Temperature evolution of the highest intensity band of pyrene at 844.4 cm − (cf. Table 2) simulated using DFT-AnharmoniCaOs and DFTB-MDTable 3: Empirical anharmonicity factors ( χ in 10 − cm − K − ) of the major IR bands of pyrene derived from the linear fittingof the band positions in different temperature ranges (Fig.7 and Figs 1 and 2 in supplementary material) calculated usingDFT-AnharmoniCaOs and DFTB-MD methods. Experimental values are listed from the condensed-phase measurements byChakraborty et al. [18] and the gas-phase measurements by Joblin et al. [9] No Theory Expt.Position AC MD Solid Gas µ m cm − Fit range χ AC Fit range χ MD Fit range χ rec [18] Fit range χ gas (K) (K) (K) (K)1 3.3 3040.9 300 - 523 − − − − − − − − − − − − − − − ± − − − ± − − ± − − − ± − − − ± − − − − − − − − − − ± − − − − − able 4: Empirical anharmonicity factors ( χ in 10 − cm − K − ) of the major IR bands of pyrene derived from the linear fittingof the bandwidths in different temperature ranges (Fig.8 and Figs 3 and 4 in supplementary material) calculated using DFT-AnharmoniCaOs and DFTB-MD methods. For DFT-AnharmoniCaOs we also report two reference values of the bandwidthat 14 and 300 K. Experimental values of χ are listed from the condensed-phase measurements by Chakraborty et al. [18] andthe gas-phase measurements by Joblin et al. [9]. “WR” and “WOR” labels the empirical anharmonicity factors respectively asderived from the gas-phase data or after subtraction of the estimated thermal rotational broadening (see text for details). No Theory Expt.Position AC MD Solid Gas µ m cm − Bandwidth Fit range χ AC Fit range χ MD Fit range χ rec [18] Fit range χ gas
14 K 300 K (K) (K) (K) (K) WR WOR1 3.3 3040.9 9.9 18.1 300 - 523 2.2 600 - 1600 1.2 - - 573 - 873 3.4 3.42 6.3 1585.3 5.6 6.3 373 - 523 0.2 600 - 1600 0.9 - -3 6.8 1472.4 6.5 10.7 300 - 523 1.5 600 - 1600 0.3 250 - 723 1.54 7.0 1428.9 4.0 6.3 300 - 523 0.7 600 - 1600 0.5 300 - 723 -5 7.6 1309.9 2.0 6.3 300 - 523 1.3 600 - 1600 0.5 - -6 8.1 1239.2 1.3 6.9 250 - 523 1.6 600 - 1600 0.8 200 - 723 0.9 ± ± ± thus more suitable for systematic use. This recipeconsists on assuming that we can describe bandpositions and widths as broken lines, with a firstlinear section defined by a data point at verylow temperature and one at ∼
300 K, and anotherlinear trend, defined by asymptotic slopes χ and χ , at higher temperatures, in agreement withall available experimental data. We showed thatwe can have several ways to access χ and χ forall bands of a given species, and we can combinethem to incrementally improve modelling, goingfrom the cheapest and less accurate DFTB-MD,to the more expensive fully quantum treatmentof DFT-AnharmoniCaOs or similar codes [e.g.Spectro 32], to the most accurate, but limitedin number of studied species, direct gas-phaselaboratory measurement.The natural development of this work will beto fill in anharmonic calculations (or suitablelaboratory measurements), at least at ∼
14 K and ∼
300 K, of absolute band positions and widths,and of χ and χ values, in a systematic way,for species in available PAH spectral databases[50, 3]. We are already in the process of developingan anharmonic PAH emission model that makesuse of this information to improve the accuracyof its predictions. In parallel, with the naturalincrease of available computational resources,we will extend AnharmoniCaOs calculations to higher temperatures, to check whether we alreadyachieved the asymptotic regime at ∼
523 K and/orat what temperature the accuracy of the quarticexpansion of the PES it employs begins to breakdown.The ultimate test will be to compare the simu-lated spectra from anharmonic AIB models withdirect measurements of UV-excited PAH infraredemission. These are very complex experiments andno progress could be achieved since the pioneeringexperiments of the 1990ies [13, 14, 15, 16, 17].While awaiting for further experimental results,the coming JWST data is a strong motivation onits own to proceed with this modelling effort.
7. Acknowledgements
The research leading to these results was sup-ported by the funding received from European Re-search Council under the European Union’s sev-enth framework program (FP/2007-2013) ERC-2013-SyG, Grant agreement n. 610256 NANOCOS-MOS. We acknowledge high performance com-puting resources made available by the “AccordoQuadro INAF-CINECA (2018)” and by CALMIP,Toulouse (project number: P16045). We thankCyril Falvo for support in obtaining the DFT data,as well as for providing the Wang-Landau sampling,that we used for our AnharmoniCaOs runs.20 eferences [1] A. Leger, L. D’Hendecourt, D. Defourneau, Astronomy& Astrophysics 216 (1989) 148–164.[2] G. Mulas, G. Malloci, C. Joblin, D. Toublanc, Astron-omy & Astrophysics 460 (2006) 93–104. doi: . arXiv:astro-ph/0606264 .[3] J. Bauschlicher, Charles W., A. Ricca, C. Boersma,L. J. Allamandola, The Astrophysical Journal Sup-plement Series 234 (2018) 32. doi: .[4] A. L. Mattioda, D. M. Hudgins, C. Boersma,J. Bauschlicher, C. W., A. Ricca, J. Cami, E. Peeters,F. S´anchez de Armas, G. Puerta Saborido, L. J. Alla-mandola, The Astrophysical Journal Supplement Series251 (2020) 22. doi: .[5] D. J. Cook, R. J. Saykally, The Astrophysical Journal493 (1998) 793–802. doi: .[6] L. Verstraete, C. Pech, C. Moutou, K. Sellgren,C. M. Wright, M. Giard, A. L´eger, R. Tim-mermann, S. Drapatz, Astronomy & Astrophysics372 (2001) 981–997. doi: . arXiv:astro-ph/0104144 .[7] C. Pech, C. Joblin, P. Boissel, Astronomy & Astro-physics 388 (2002) 639–651. doi: .[8] G. Mulas, G. Malloci, C. Joblin, D. Toublanc, Astron-omy & Astrophysics 456 (2006) 161–169. doi: . arXiv:astro-ph/0605411 .[9] C. Joblin, P. Boissel, A. Leger, L. D’Hendecourt, D. De-fourneau, Astronomy & Astrophysics 299 (1995) 835.[10] J. Kurtz, Astronomy & Astrophysics 255 (1992) L1–L4.[11] J. Oomens, A. G. G. M. Tielens, B. G. Sartakov,G. von Helden, G. Meijer, The Astrophysical Journal591 (2003) 968–985. doi: .[12] J. Zhen, A. Candian, P. Castellanos, J. Bouwman,H. Linnartz, A. G. G. M. Tielens, The AstrophysicalJournal 854 (2018) 27. doi: . arXiv:1809.08800 .[13] I. Cherchneff, J. R. Barker, The Astrophysical JournalLetters 341 (1989) L21. doi: .[14] J. Shan, M. Suton, L. C. Lee, The Astrophysical Journal383 (1991) 459. doi: .[15] J. Brenner, J. R. Barker, The Astrophysical JournalLetters 388 (1992) L39. doi: .[16] S. Schlemmer, D. Cook, J. Harrison,B. Wurfel, W. Chapman, R. Saykally, Sci-ence 265 (1994) 1686–1689. URL: https://science.sciencemag.org/content/265/5179/1686 . doi: . arXiv:https://science.sciencemag.org/content/265/5179/1686.full.pdf .[17] D. J. Cook, S. Schlemmer, N. Balucani, D. R.Wagner, J. A. Harrison, B. Steiner, R. J.Saykally, The Journal of Physical ChemistryA 102 (1998) 1465–1481. URL: https://doi.org/10.1021/jp9724434 . doi: . arXiv:https://doi.org/10.1021/jp9724434 , pMID:11542815.[18] S. Chakraborty, G. Mulas, K. Demyk, C. Joblin, J.Phys. Chem. A 123 (2019) 4139–4148.[19] D. Rouan, A. Leger, A. Omont, M. Giard, Astronomy& Astrophysics 253 (1992) 498–514.[20] G. Mulas, Astronomy & Astrophysics 338 (1998) 243–261.[21] G. Malloci, G. Mulas, C. Joblin, D. Toublanc, I. Porceddu, in: A. Wilson (Ed.), ESA Special Pub-lication, volume 577 of ESA Special Publication , pp.385–386.[22] N. T. Van Oanh, P. Parneix, P. Br´echignac, J. Phys.Chem.A. 106 (2002) 10144–10151.[23] M. Rapacioli, A. Simon, C. C. M. Marshall, J. Cuny,D. Kokkin, F. Spiegelman, C. Joblin, J. Phys. Chem. A119 (2015) 12845–12854.[24] B. Joalland, M. Rapacioli, A. Simon, C. Joblin, C. J.Marsden, F. Spiegelman, The Journal of PhysicalChemistry A 114 (2010) 5846–5854. URL: https://doi.org/10.1021/jp911526n . doi: . arXiv:https://doi.org/10.1021/jp911526n , pMID:20394399.[25] A. Simon, M. Rapacioli, M. Lanza, B. Joal-land, F. Spiegelman, Phys. Chem. Chem. Phys. 13(2011) 3359–3374. URL: http://dx.doi.org/10.1039/C0CP00990C . doi: .[26] A. Simon, M. Rapacioli, J. Mascetti, F. Spiegel-man, Phys. Chem. Chem. Phys. 14 (2012) 6771–6786.doi: {10.1039/c2cp40321h} .[27] A. Simon, C. Iftner, J. Mascetti, F. Spiegelman, J.Phys. Chem. A 119 (2015) 2449–2467. URL: https://doi.org/10.1021/jp508533k . doi: .[28] M. Piccardo, J. Bloino, V. Barone, In-ternational Journal of Quantum Chem-istry 115 (2015) 948–982. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/qua.24931 . doi: https://doi.org/10.1002/qua.24931 . arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.24931 .[29] O. Pirali, M. Vervloet, G. Mulas, G. Malloci, C. Joblin,Physical Chemistry Chemical Physics (IncorporatingFaraday Transactions) 11 (2009) 3443. doi: .[30] M. Basire, P. Parneix, F. Calvo, T. Pino, P. Br´echignac,Journal of Physical Chemistry A 113 (2009) 6947–6954.doi: .[31] A. K. Lemmens, D. B. Rap, J. M. M. Thunnissen,C. J. Mackie, A. Candian, A. G. G. M. Tielens,A. M. Rijs, W. J. Buma, Astronomy & Astrophysics628 (2019) A130. doi: . arXiv:1907.09351 .[32] C. J. Mackie, A. Candian, X. Huang, E. Maltseva,A. Petrignani, J. Oomens, A. L. Mattioda, W. J. Buma,T. J. Lee, A. G. G. M. Tielens, J. Chem. Phys. 145(2016) 084313. doi: .[33] F. Calvo, M. Basire, P. Parneix, The Journal of Physi-cal Chemistry A 115 (2011) 8845–8854. URL: https://doi.org/10.1021/jp202935p . doi: .[34] G. Mulas, C. Falvo, P. Cassam-Chena¨ı, C. Joblin,J. Chem. Phys. 149 (2018) 144102. doi: . arXiv:1809.05669 .[35] E. Maltseva, A. Petrignani, A. Candian, C. J. Mackie,X. Huang, T. J. Lee, A. G. G. M. Tielens, J. Oomens,W. J. Buma, The Astrophysical Journal 814 (2015) 23.doi: . arXiv:1510.04948 .[36] E. Maltseva, A. Petrignani, A. Candian, C. J. Mackie,X. Huang, T. J. Lee, A. G. G. M. Tielens, J. Oomens,W. J. Buma, The Astrophysical Journal 831 (2016) 58.doi: .[37] F. Wang, D. P. Landau, Phys. Rev. E 64(2001) 056101. doi: . arXiv:cond-mat/0107006 .[38] F. Wang, D. P. Landau, Phys. Rev. Lett. 86(2001) 2050–2053. doi: . rXiv:cond-mat/0011174 .[39] G. Seifert, D. Porezag, T. Frauenheim, Int. J. QuantumChem. 58 (1996) 185–192.[40] D. Porezag, T. Frauenheim, T. K¨ohler, G. Seifert,R. Kaschner, Phys. Rev.,
B 51 (1995) 12947–12957.[41] M. Elstner, D. Porezag, G. Jungnickel, J. Elsner,M. Haugk, T. Frauenheim, S. Suhai, G. Seifert, Phys.Rev. B 58 (1998) 7260–7268.[42] T. Heine, M. Rapacioli, S. Patchkovskii, J. Frenzel,A. M. K¨oster, P. Calaminici, H. A. Duarte, S. Es-calante, R. Flores-Moreno, A. Goursot, J. U. Reve-les, D. R. Salahub, A. Vela, demonnano, 2020. 2020,http://demon-nano.ups-tlse.fr/.[43] J. Frenzel, A. F. Oliveira, N. Jardillier, T. Heine,G. Seifert, Semi-relativistic, self-consistent chargeSlater-Koster tables for density-functional based tight-binding (DFTB) for materials science simulations.,2004-2009.[44] M. Rapacioli, F. Spiegelman, D. Talbi, T. Mineva,A. Goursot, T. Heine, G. Seifert, J. Chem. Phys. 130(2009) 244304–10. URL: http://link.aip.org/link/?JCP/130/244304/1 .[45] M. P. Gaigeot, M. Sprik, J. Phys. Chem.B 107 (2003)10344–10358. doi: {10.1021/jp034788u} .[46] W. G. Hoover, Phys. Rev. A 31 (1985) 1695–1697.[47] S. Nos´e, J. Chem. Phys. 81 (1984) 511. URL: http://scitation.aip.org/content/aip/journal/jcp/81/1/10.1063/1.447334 . doi: .[48] M. Tuckermana, R. April, G. J. Martyna, M. L. Klein,J. Chem. Phys. 97 (1992) 2635–2644.[49] C. Joblin, L. D’Hendecourt, A. Leger, D. Defourneau,Astron. Astrophys 281 (1994) 923–936.[50] G. Malloci, C. Joblin, G. Mulas, Chemical Physics 332(2007) 353–359. doi: . arXiv:astro-ph/0701254 .[51] B. E. Brumfield, J. T. Stewart, B. J. Mc-Call, The Journal of Physical Chemistry Let-ters 3 (2012) 1985–1988. URL: https://doi.org/10.1021/jz300769k . doi: . arXiv:https://doi.org/10.1021/jz300769k .[52] C. J. Mackie, T. Chen, A. Candian, T. J. Lee,A. G. G. M. Tielens, J. Chem. Phys. 149 (2018) 134302.doi: . arXiv:1810.01975 . upporting information for:Modeling Anharmonic Infrared Spectra of ThermallyExcited Pyrene (C H ): the combined view of DFTAnharmonicCaOs and approximate DFT moleculardynamics Shubhadip Chakraborty a , Giacomo Mulas a,b, ∗ , Mathias Rapacioli c , ChristineJoblin a a Institut de Recherche en Astrophysique et Plan´etologie, Universit´e de Toulouse (UPS),CNRS, CNES, 9 Av. du Colonel Roche, 31028 Toulouse Cedex 4, France. b Istituto Nazionale di Astrofisica (INAF), Osservatorio Astronomico di Cagliari, 09047Selargius (CA), Italy c Laboratoire de Chimie et Physique Quantiques (LCPQ/IRSAMC), Universit´e de Toulouse(UPS),CNRS, 118 Route de Narbonne, 31062 Toulouse, France
Table 1:
List of harmonic frequencies of pyrene using DFT, DFTB and anharmonicfrequencies calculated using AnharmoniCaOs package at 0 K.
Sym No DFT DFTBFreq Int Freq Int(cm − ) (km mol − ) (cm − ) (km mol − )A g ν u ν g ν u ν g ν u ν u ν g ν ∗ Corresponding author
Preprint submitted to Journal of Molecular Spectroscopy February 15, 2021 a r X i v : . [ a s t r o - ph . GA ] F e b able 1 – continued from previous page Sym No DFT DFTBFreq Int Freq Int(cm − ) (km mol − ) (cm − ) (km mol − )B g ν u ν g ν u ν u ν g ν g ν g ν u ν u ν u ν u ν g ν g ν g ν g ν u ν u ν g ν g ν u ν u ν g ν g ν u ν able 1 – continued from previous page Sym No DFT DFTBFreq Int Freq Int(cm − ) (km mol − ) (cm − ) (km mol − )B g ν u ν g ν u ν g ν u ν u ν g ν u ν g ν u ν u ν g ν u ν g ν g ν g ν u ν g ν u ν u ν u ν g ν g ν u ν able 1 – continued from previous page Sym No DFT DFTBFreq Int Freq Int(cm − ) (km mol − ) (cm − ) (km mol − )B g ν g ν g ν u ν u ν g ν g ν u ν u ν g ν g ν u ν u ν u ν able 2: Integration range for deriving band positions in spectra generated byAnharmoniCaOs package
No AnharmoniCaOsPosition at 300 K Range1 3040 3060-30202 1585 1605-16563 1473 1493-14534 1424 1444-14045 1311 1331-12916 1241 1261-12217 1180 1200-11608 1089 1109-10699 996 1016-97610 964 984-94411 843 863-82312 816 836-79613 739 760-73014 710 730-690Table 3:
Integration range for deriving band positions in spectra generated bydeMonNano package
Band position Temperature (K)at 300K 200 300 400 600 800 1000 1200 1400 1600723.7 740 - 715 740 - 715 740 - 712 740 - 709 740 - 706 740 - 703 740 - 700 740 - 697 740 - 696814.1 830 - 770 830 -770 830 - 770 830 - 770 830 - 770 830 -770 830 - 770 830 -770 830 -7701089.2 1098 - 1084 1095 - 1081 1092 - 1078 1090 - 1073 1088 - 1068 1084 - 1062 1082 - 1058 1080 - 1053 1078 - 10511333.1 1354 - 1323 1354 - 1321 1354 - 1319 1354 - 1312 1354 - 1308 1354 - 1304 1354 - 1298 1354 - 1291 1354 - 12911465.8 1481 - 1462 1476 - 1454 1468 - 1453 1458 - 1436 1453 - 1430 1446 - 1425 1440 - 1414 1430 - 1400 1430 - 14001628.5 1643 - 1624 1640 - 1612 1637 - 1604 1631 - 1597 1623 - 1583 1617 - 1567 1603 - 1552 1591 - 1552 1580 - 15401585.5 1597 - 1582 1592 - 1576 1590 - 1572 1583 - 1566 1568 - 1533 1566 - 1522 1561 - 1517 1547 -1506 1539 - 15001842.4 1862 - 1840 1858 - 1829 1855 - 1826 1847 - 1779 1835 - 1779 1829 - 1769 1813 - 1740 1804 - 1736 1804 - 17142923.9 3000 -2700 3000 - 2700 3000 - 2700 3000 - 2700 3000 - 2700 3000 - 2700 3000 - 2700 3000 - 2700 3000 - 2700 . Illustration of the variation of band positions with temperature ob-tained from AnharmoniCaOs (DFT) and deMonNano (MD) sim-ulationAnharmoniCaOs Temperature (K) B a n d p o s i t i o n ( c m ) Figure 1: Evolution of band positions with temperature obtained from AnharmoniCaOs. Thelinear trend of the band position with temperature is applicable only for the high temperaturerange. Therefore only the points from and beyond 300 K have been considered for the fitting.Band position at 300 K are marked in each panel of the figure.
100 200 300 400 50010861088109010921094 1090.1 0 100 200 300 400 500117811801182118411861188 1183.40 100 200 300 400 50012361238124012421244 1239.2 0 100 200 300 400 50013061308131013121314 1309.90 100 200 300 400 50014261428143014321434 1428.9 0 100 200 300 400 500146814701472147414761478 1472.40 100 200 300 400 50015751580158515901595 1585.3 0 100 200 300 400 50030303035304030453050 3040.9
Temperature (K) B a n d p o s i t i o n ( c m ) Figure 1: Evolution of band positions with temperature obtained from AnharmoniCaOs. Thelinear trend of the band position with temperature is applicable only for the high temperaturerange. Therefore only the points from and beyond 300 K have been considered for the fitting.Band position at 300 K are marked in each panel of the figure. (Continued from the previouspage) eMonNano Temperature (K) B a n d p o s i t i o n ( c m ) Figure 2: Evolution of band positions with temperature obtained from MD (demonNano)simulation. This method is dedicated to simulate the band positions at high temperature.Therefore only the points from and beyond 600 K have been considered for the fitting. Bandposition at 300 K are marked in each panel of the figure.
250 500 750 1000 1250 1500 1750 20001400142514501475150015251550 1465.8 0 250 500 750 1000 1250 1500 1750 20001550157516001625165016751700 1628.50 250 500 750 1000 1250 1500 1750 20001500152515501575160016251650 1585.5 0 250 500 750 1000 1250 1500 1750 20001750177518001825185018751900 1842.40 250 500 750 1000 1250 1500 1750 20002800282528502875290029252950 2923.9
Temperature (K) B a n d p o s i t i o n ( c m ) Figure 2: Evolution of band positions with temperature obtained from MD (demonNano)simulation. This method is dedicated to simulate the band positions at high temperature.Therefore only the points from and beyond 600 K have been considered for the fitting. Bandposition at 300 K are marked in each panel of the figure. (Continued from the previous page)
2. Illustration of the variation of bandwidths with temperature ob-tained from AnharmoniCaOs (DFT) and deMonNano (MD) sim-ulationAnharmoniCaOs
100 200 300 400 500024681012 711.6 0 100 200 300 400 500024681012 741.10 100 200 300 400 500024681012 817.8 0 100 200 300 400 500024681012 844.40 100 200 300 400 500024681012 961.5 0 100 200 300 400 500024681012 994.9
Temperature (K) B a n d w i d t h ( c m ) Figure 3: Evolution of bandwidths with temperature obtained from AnharmoniCaOs. Thelinear trend of the band position with temperature is applicable only for the high temperaturerange. Therefore only the points from and beyond 300 K have been considered for the fitting.Band position at 300 K are marked in each panel of the figure.
100 200 300 400 500024681012 1090.1 0 100 200 300 400 500024681012 1183.40 100 200 300 400 500024681012 1239.2 0 100 200 300 400 500024681012 1309.90 100 200 300 400 500024681012 1428.9 0 100 200 300 400 500024681012 1472.40 100 200 300 400 5000246810 1585.3 0 100 200 300 400 5001012141618 3040.9
Temperature (K) B a n d w i d t h ( c m ) Figure 3: Evolution of bandwidths with temperature obtained from AnharmoniCaOs. Thelinear trend of the band position with temperature is applicable only for the high temperaturerange. Therefore only the points from and beyond 300 K have been considered for the fitting.Band position at 300 K are marked in each panel of the figure. (Continued from the previouspage) eMonNano
600 800 1000 1200 1400 1600 1800051015202530 723.7 600 800 1000 1200 1400 1600 1800051015202530 814.1600 800 1000 1200 1400 1600 1800051015202530 1089.2 600 800 1000 1200 1400 1600 1800051015202530 1331.1
Temperature (K) B a n d w i d t h ( c m ) Figure 4: Evolution of bandwidths with temperature obtained from MD (demonNano) simula-tion. This method is dedicated to simulate the band positions at high temperature. Thereforeonly the points from and beyond 600 K have been considered for the fitting. Band position at300 K are marked in each panel of the figure.
00 800 1000 1200 1400 1600 1800051015202530 1465.8 600 800 1000 1200 1400 1600 1800051015202530 1628.5600 800 1000 1200 1400 1600 1800051015202530 1585.5 600 800 1000 1200 1400 1600 1800051015202530 1842.4600 800 1000 1200 1400 1600 180040455055606570 2923.9