Coronene: a model for ultrafast dynamics in graphene nanoflakes and PAHs

Assuming a delta pulse excitation, quantum wavepackets are propagated on the excited state manifold in the energy range from 3.4-5.0 eV for coronene and 2.4-3.5 eV for circumcoronene to study the time evolution of the states as well as their lifetimes. The full-dimensional (102 and 210 degrees of freedom for coronene and circumcoronene respectively) non-adiabatic dynamics simulated with the ML-MCTDH method on twelve coupled singlet electronic states show that the different absorption spectra are only due to electronic delocalisation effects that change the excited state energies, but the structural dynamics in both compounds are identical. Breathing and tilting motions drive the decay dynamics of the electronic states away from the Frank-Condon region independently of the size of the aromatic system. This promising result allows the use of coronene as a model system for the dynamics of larger polycyclic aromatic hydrocarbons (PAHs) and graphene one dimensional sheets or nanoflakes.


I. INTRODUCTION
Photoinduced processes play a crucial role in many fields of chemistry, physics, biology, and material sciences, ranging from their involvement in vision to anticancer therapies and optoelectronics.[1][2][3][4][5] Understanding such processes is essential to be able to design new compounds and tune the optical and electronic properties of materials, enabling a desired functionalisation.Particularly in the field of optoelectronic devices such as solar cells, lightemitting diodes (LEDs), and lasers, photoinduced processes are fundamental, influencing their performance and efficiency.To conveniently tune their performance, one must understand the behavior of the excited states of the target materials and study their time evolution.Studying these processes enhances our understanding of fundamental physical phenomena, facilitates the development of new theoretical models, and guides the design of novel materials with tailored properties.[6][7][8] Polycyclic aromatic hydrocarbons (PAHs) are π-conjugated organic molecules that consist of multiple fused aromatic rings.They are naturally found as products from the incomplete combustion of organic materials, forest fires, and volcanic eruptions.[9,10] Leaving aside their threat to humans health, [11] in the last years they have become highly relevant in industrial applications such as the production of dyes and pigments and pharmaceuticals [12,13] due to their photostability and photophysical properties.[14,15] These compounds have also been found to be a key piece in the understanding the chemistry of the interstellar medium.[16] Moving towards the solid state, graphene quantum dots (GQDs) are nanostructures composed of small fragments of graphene, which is a two-dimensional carbon allotrope, typically showing lateral dimensions ranging from a few to several hundred nanometers.[17] GQDs possess unique optical, electronic, and chemical characteristics, distinct from bulk graphene, presumably arising from their size, edge structures, and quantum confinement effects.[18][19][20] In recent years, it has been proposed that their tunable photoluminescence may be related to the method of synthesis, attributed to doping effects, surface functional groups, and the presence of molecular fluorophores, [21][22][23][24][25][26] which has been widely investigated also by theoretical computations, recently reviewed by M. Langer et al. in Ref. 27.As a consequence of their ashtonishing properties, GQDs have applications in optoelectronics, [22,[28][29][30] energy storage, [31,32] sensing, [20,33] and biomedicine.[34,35] Graphene nanoflakes, also known as graphene nanoribbons or nanostrips, are narrow and elongated structures of graphene, and similarly to GQDs, they can be synthesized with different widths, lengths, and edge structures, resulting in diverse electronic properties.[36,37] Their edge-state-dependent properties are due to quantum confinement along the width of the nanoflakes.In general, they are larger than GQDs, and therefore the quantum confinement effects may be less pronounced.
Coronene and circumcoronene are both examples of polycyclic aromatic hydrocarbons that are often used as models for studying the properties of larger PAHs, GQDs, and graphene nanoflakes.[38] Coronene is a PAH consisting of a planar molecule with a hexagonal shape, composed of seven fused benzene rings.It has a well-defined structure and is often used as a model system to study the properties of extended aromatic compounds due to its relative simplicity.Coronene is known for its fluorescence properties and has been extensively studied for its electronic, optical, and chemical properties.[39] It is also used as a benchmark molecule in theoretical calculations and experimental studies.[40][41][42][43] UV-VIS results show a band at 4.1 eV, with a shoulder displaced to the blue (4.3 eV).A forbidden band, allowed only by symmetry breaking, appears at 3.6 eV.
Circumcoronene is a larger PAH composed of a circular arrangement of benzene rings.
It is formed by the fusion of multiple coronene units, resulting in a larger and more complex structure.It contains exactly double the amount of carbon and hydrogen atoms than coronene (48 carbons and 24 hydrogens).Due to its size and extended π conjugation, circumcoronene exhibits unique electronic and optical properties.It is often employed as a model system to investigate the properties of graphene nanoflakes because of its structural similarity (and its larger size compared to coronene) and the possibility of tuning its electronic structure through chemical modifications.The synthesis of circumcoronene has always remained a challenge and it has been recently synthesized in solution, [44] although a UV-VIS spectrum is not yet available.
Both coronene and circumcoronene have been widely studied in various fields, including materials science, organic electronics, and nanotechnology, as they provide valuable insights into the behavior of PAHs and graphene nanoflakes.Due to their reduced size compared to the material structures of interest, one can use more accurate theoretical methods to deeply understand their properties.Therefore, they serve as important model systems for understanding the fundamental properties and potential applications of larger carbon-based materials.[38] Although there are excited state dynamical studies on smaller PAHs, such as benzene, [45,46] naphthalene dimer, [47] anthrazene [48] and pyrene, [49] the field is lacking an exhaustive work on the excited state behavior of coronene and circumcoronene.The purpose of this article is to compare the excited state dynamics on coronene and circumcoronene molecules, to assess the feasibility of using the coronene molecule as a model system for larger polycyclic aromatic hydrocarbons.

A. Electronic structure details
To compute vertical energies, gradients and couplings needed throughout this work, we benchmarked several TD-DFT functionals and basis sets against the wave function method EOM-CCSD/cc-pVDZ for the coronene molecule.For the CCSD ground state reference, the symmetry was chosen to be D 2h , the highest symmetry available.The Q-Chem package [50] was used to compute the reference EOM-CCSD vertical energies, whereas and Gaussian 16 [51] was chosen to compute the vibrational frequencies of both systems.For every other calculation, i.e., Wigner-displaced geometries, parametrisation of the vibronic coupling Hamiltonian and on-the-fly dynamics, we used the ORCA 5.0 electronic structure program.[52,53] The details about the benchmark we performed are included in the Supplementary Information, and they are supported by a previous work by B. Shi et al. [38] The TD-B3LYP-D3/def2-SVP method was proven to provide an excellent agreement with respect to experimental results while keeping the computational time feasible.For circumcoronene, we used the same combination to compute the parametrised potentials: TD-B3LYP-D3/def2-SVP.

B. The vibronic coupling model
To describe the potential energy surfaces of the systems, we constructed an approximate Hamiltonian as a first order Taylor expansion of the vertical electronic state energies calculated at the Frank-Condon (FC) point.
A diabatic basis ensures a smooth description of the potential energy surfaces and massand frequency-weighed normal mode coordinates Q α (α = 1, . . ., f ), where f is the number of vibrational degrees of freedom (DOFs) of the system, are chosen as the coordinate basis to expand the potential energy surfaces.Tuning modes carry the gradient, driving the molecules away from the FC region, whereas coupling modes mix pairs of electronic states and are responsible for the population transfer.When this expansion is truncated in the first order, the method is known as the linear vibronic coupling (LVC) model, and the Hamiltonian looks as follows: (1) . ( The zeroth-order Hamiltonian H H H (0) is defined by the diagonal kinetic energy operator of the N electronic states to be included, A simple way of defining the zeroth-order potential matrix W W W (0) is to use the harmonic approximation shifted by the vertical excitation energies E i at the Frank-Condon (FC) geometry while W W W (1) , the first-order corrections of the potential matrix are defined using the following expansion: The κ α linear intra-state vibronic couplings carry the gradients (along the α tuning modes) of the i excited potential energy surface at the FC point.The λ (ij) α linear inter-state vibronic constants measure the coupling between each pair of electronic states i and j along the normal mode α.
In this work, we used the LVC model as implemented in the SHARC molecular dynamics package [54][55][56], which generates all the LVC model parameters [57] from (a) the calculation of the ground-state frequencies at the gas-phase optimised geometries of coronene and circumcoronene and (b) 2(3N − 6) + 1 excited state energy calculations of displaced geometries in positive and negative directions along every vibrational mode from the equilibrium configuration.From these computations, diabatic state gradients and couplings were computed from numerical derivatives, as detailed in Ref. 57.As we are dealing with systems of nonabelian symmetry groups, a rotation of the Hamiltonian was applied to the Hamiltonian to correct spurious couplings, as described in Ref. 58.

C. The ML-MCTDH wavepacket propagation method
The nuclear dynamics of the system can be studied by solving the time-dependent Schrödinger equation (TDSE), where the Hamiltonian operator H can be generalized to describe a system characterized by a set of coupled potential energy surfaces in the diabatic picture, like in Eq. 1.The total wavefunction of the system is a function of the normal modes and its definition depends on the method we use to solve this first-order differential equation.
A common approach to solve the TDSE is the multiconfigurational time dependent Hartree (MCTDH) method, [59,60] which expands the total wave function in a set of time- which are known as single-particle functions (SPFs) and are expanded in a time-independent primitive basis set χ (α) , constructed using a DVR [61] grid, Due to the exponential scaling of the computational effort with respect to the dimensionality of the system, MCTDH is normally used for systems with no more than 30 vibrational degrees of freedom.For higher-dimensional systems, the multi-layer MCTDH (ML-MCTDH) [62] approach can be used, where a recursive layering structure is constructed in a way that the SPFs are expanded in a set of time-dependent SPFs and only the last layer is expanded using the time-independent DVR primitive basis, still ensuring perfect convergence to the exact quantum dynamics result at the nuclear basis limit.
All these methodologies are the core of the Quantics package [63,64], used throughout this work to study the quantum dynamical properties of coronene and circumcoronene.

D. Trajectory Surface Hopping
The Tully surface hopping (TSH) [65,66] method propagates the nuclei following the Newton equations of motion where the atom A is moving in a single potential energy surface (V i (R(t))) pushed by forces corresponding to its gradient.Compared to MCTDH, instead of a full nuclear wavepacket, TSH consists of individual trajectories that move completely uncoupled from the bundle.
The propagation from one-time step to another is performed using the velocity Verlet algorithm, [67] where velocities v A (t) and the atomic coordinates R A (t) change according to The nuclei are assumed to move adiabatically most of the time, and in some particular regions where the coupling with other electronic states is high, undergo non-adiabatic transitions known as hops, thus the name surface hopping.[65,66] At every time step, the nuclei have the choice to change the electronic state based on a probability to hop that depends on the coefficients of the electronic states involved.The electronic wavefunction is therefore expressed as a sum of different electronic states with time-dependent coefficients using the Born-Huang expansion.Inserting this wavefunction ansatz in the TDSE yields to the following equations for the propagation of the coefficients: where the Hamiltonian term is directly calculated from the precomputed PES or via electronic structure methods, and the second term is calculated from the NAC between electronic states and velocities: Throughout this work, we have used the SHARC implementation of TSH [68] both calculating electronic quantities on-the-fly and on precomputed potentials, as well as the more accurate wavepacket method ML-MCTDH on the same potentials.The reason for using the less-accurate method TSH was twofold: first, to assess the role of triplets since the spin orbit coupling was negligible at the Frank-Condon geometry and, second, to compare whether dynamics on precomputed potentials were a good approximation to the dynamics computed on-the-fly.The TSH calculations used the "fewest switches" algorithm, [66] along with the decoherence correction of Granucci and Persico.[69] The time-step was chosen to be 0.5 fs.The hopping probabilities were obtained from state overlaps.[70] 400 initial structures and velocities were generated from a Wigner distribution.
More system-specific computational details are provided at the begining of the corresponding sections.

III. RESULTS: QUANTUM DYNAMICS OF CORONENE A. Vertical absorption and static results
Using the level of theory introduced in Sec.II A, we computed the corresponding vertical excitations from the singlet ground state (X, S 0 ) of coronene, collected in Table I in energetic order following the TD-B3LYP-D3 results.Besides the experimental reference, [41] Table I includes also the electronic character of every excited state, providing the molecular orbitals (MOs) with the largest contributions.The graphical representations of these MOs are included in Fig. 3, where the last 4 occupied orbitals are in the bottom figure, and the first 4 virtual orbitals are in the top row.The MOs represented in Fig. 3 are degenerate in the D 2h symmetry group in the following pairs of MOs: (LUMO+3, LUMO+2), (LUMO+1, LUMO), (HOMO, HOMO−1), (HOMO−2, HOMO−3).Among all the vertical transitions involving these orbitals, only two of them show oscillator strength at 4.4 eV, which corresponds to the excitation to the 2B 2u (S 6 ) and 2B 3u (S 7 ) singlet excited states, matching the maximum absorption peak of the experimental spectrum around 4 eV.These two bright states are degenerate in D 2h , and they correspond to a E 1u state in D 6h .[38] There is another absorption peak in the experimental spectrum at 3.6 eV which appears dark in the FC region and corresponds to the excitation to the 1B 3u (S 2 ) state.
The transition to this state is forbidden by symmetry, but it becomes accessible through the vibronic couplings between states.This effect has been observed by generating a set of initial conditions displaced from the equilibrium configuration using a Wigner distribution based on the S 0 state frequencies computed at the TD-B3LYP-D3/def2-SVP level of theory.To study the dynamical evolution of the diabatic state manifold, excited state dynamics have been carried out with different methods from both bright excited states, S 2 and S 6 /S 7 , as well as from the dark states to show how well it compares to non-adiabatic dynamics in the circumcoronene system.
For the case of the TSH, the 54 on-the-fly trajectories were initialised based on oscillator strengths using an energetic window of 3.0 to 3.75 eV to excite to the first "forbidden" band and performed on four coupled singlet electronic states (S 0 -S 4 ) and six triplets (T 1 -T 6 ) to minimise the computational effort needed.The TSH dynamics on LVC potentials started from the same geometries as the on-the-fly dynamics and included 54 trajectories on twelve singlet and eleven triplet electronic states.The ML-MCTDH dynamics were initialised as the hartree product of the vibrational eigenstate of every degree of freedom.The dynamics on the twelve coupled electronic states was propagated using Runge-Kutta 5 for the nuclear degrees of freedom and short iterative Lanczos (SIL) for the electronic one.4.6×10 7 primitive basis functions were used for the grid and the calculations included 1.2×10 35 configurations.
For the 8D reduced MCTDH dynamics the Bulirsch-Stoer extrapolation scheme was used for the propagation of the SPFs and the SIL integrator for the A-coefficients and the electronic propagation.The total size of the grid for the converged results in this system was 3.25×10 9   and used 789 nuclear configurations.
The experimental gas phase absorption spectrum of coronene was recorded by S. Hirayama et al., [41] where the first absorption band of coronene was reported around ∼3.6 eV, which corresponds to excite the system to the second singlet excited state, 1B 3u (S 2 ).
The population transfer from this excited state is shown in Fig. 2   In Fig. 2, the main difference between panels (a) and (b) is the PESs on which the classical nuclear trajectories were propagated.The more similar those panels, the better the LVC approximation in the PES construction.It is important to note that diabatic populations in panel (a) are not truly diabatic, but that the states have been labelled according to the diabatic state distribution that was present at the initial geometries, in a diabatic by ansatz manner.To clarify this fact, the combination of the S 1 and S 2 adiabatic populations is included in both (a) and (b) panels of Fig. 2, were a good agreement can be found.TSH in this case does not give the same diabatic results on-the-fly or on precomputed potentials and both panels are substantially different.This may be surprising since the LVC model has been proven to be a good approximation in rigid systems with impeded torsions, as it is the case In Fig. 2, the difference between panels (b) and (c) lies in the quantum treatment of the atomic nuclei in MCTDH in contrast to classical propagations during the TSH/LVC approach.From the comparison of these two panels, we can conclude that the TSH/LVC approach is a good approximation that covers qualitatively the population changes.However, the fine print of the dynamics of coronene is only provided by the most accurate method ML-MCTDH.
For the ML-MCTDH method, all 102 vibrational modes of the system have been included.
The ML-tree structure, available on the Supplementary Information, has been optimized by looking at the largest frequency-weighted κ α gradients and λ (ij) α couplings among states which involved the initial excited state 1B 3u (S 2 ).[71] By using the same arguments, we can analyze those vibrations (normal modes) that contribute the most to the structure relaxation and drive the nuclear dynamics.
The identification of the leading vibrations in the dynamics is crucial, not only for unraveling the nuclear dynamics of coronene but for validating this methodology for future work with larger systems like graphene nanoflakes, and to understand the absorption and emission properties of these systems in general.For the case of coronene, the decay from S 2 is mostly due to the totally symmetric breathing mode ν 24 , which shows a negligible diabatic coupling between S 2 and S 1 but large gradient on the S 2 state.Another leading mode in this dynamics is ν 19 , the scissoring mode exhibiting the largest diabatic coupling.Two cuts of the diabatic PES along these two modes, along with their corresponding diabatic couplings are represented in the top panels of Fig. 5.The full-dimensional (full-D) dynamics have been compared with several reduced-dimensionality models, where only the most important modes have been included in a MCTDH calculation.We found the most relevant 8 modes which are the main responsible for the spectral features of the main band and during the first 10 fs of the dynamics.We systematically added modes to this 8D model trying to find a larger system to recover the full-D dynamics up to 50 fs.This was not possible, since there are many modes which contribute with small couplings to the Hamiltonian and cannot be left out.The spectrum resulting from the 8D dynamics and state populations compared to the full D results are depicted in the SI.The eight modes included are the three main breathing modes (full breathing v24, inner benzene breathings v60 and v77), four tilting modes (v19, v64, v79 and v88) and a CH stretching mode (v108).Note that the numeration of the modes started at v7 after removing the six translations and rotations.
After the vertical excitation from the ground state (GS) to the 1B 3u (S 2 ) state, the diabatic population transfer from this excited state is represented in Fig. 6 (top left panel) to up to 50 fs, where it can be seen that the population is equally distributed between S 2 and S 1 states already after the first 15 fs.When classical trajectories are used for the nuclei motion, a population transfer is also observed from the S 2 to S 1 state.However, the time scales observed in panels b) and c) of Fig. 2 are slightly different.Slower time scales are often observed when running classical trajectories, since they need to move towards the intersection, whereas a wavepacket has grid functions i.e., the possibility of population, in every part of the conformational space.From this comparison, it is suggested that the quantum effects play an important role in the nuclear motion over the PESs.We could also extract that electronic populations may not be the best description when dealing with high symmetry molecules, since arbitrary rotations impede a correct assignment of state characters.

C. Excited state quantum dynamics from every electronic state
To complete the dynamical study while considering the results obtained from the analysis of the quantum dynamics from the 1B 3u (S 2 ), ML-MCTDH has been carried out to unravel the quantum dynamics from every electronic excited state of coronene in the relevant energy range (below 5 eV).By exciting the system to each of the first 8 singlet states, we followed the corresponding diabatic populations.The corresponding figures of these population transfers are collected in the Supplementary Information (Figs.S9-S16).
Since the excited states 2B 3u (S 6 ) and 2B 2u (S 7 ) are responsible of the most intensive peak in the experimental absorption spectrum, the quantum dynamics from these states are especially relevant.Due to the symmetry reduction in this calculations, 2B 3u and 2B 2u , which are degenerate in D 2h , correspond to the same state in the real symmetry group D 6h .This is confirmed by the fact that the quantum dynamics from both excited states are almost identical.After the excitation to these states, ∼ 70% of the population is transferred during the first 50 fs.From 2B 3u (S 6 ), the decay seems slightly faster as compared to starting the dynamics on the 2B 2u (S 7 ).Similar behavior is found for the diabatic population transfer from 1B 1g (S 3 ) and 1A g (S 4 ) states, which are also degenerated in D 2h .The likeness between them is stronger during the first 10 fs, going from S 3 and S 4 to the degenerated pair S 10 and S 11 respectively, caused by the energy proximity and coupling among these states after vibrational relaxation.The first two excited states are also populated but the decay is faster from S 3 .The diabatic transfer from 2A g (S 5 ) and 2B 1g (S 8 ) are the fastest, where more than 70% of the population is lost before 10 fs.

IV. RESULTS: EXCITED STATE QUANTUM DYNAMICS OF CIRCUMCORONENE AND COMPARISON TO CORONENE
If we think about coronene as the smallest portion of a graphene sheet including all the symmetry operations (D 6h ), then circumcoronene is the next step towards the generalization of graphene retaining the highest symmetry.Thus, we extended our study to this system, with twice as many atoms (72) as coronene and 210 vibrational modes.

A. Vertical absorption and static results
For consistency, the same methodology and level of theory introduced in Sec.II A, and used for coronene in Sec.III A, has been used to characterize the diabatic PES of circumcoronene.The vertical excitations have been computed also for this bigger graphene nanoflake using TD-DFT and included in Table II, similarly to the analysis of coronene in Table I, except for the absence of the EOM-CCSD and experimental references.Although this system has been synthesized recently by R. Gershoni-Poranne et al., [72] there is not an experimental absorption and emission spectrum, yet.One very interesting outcome of this work is the fact that the molecular orbitals involved in the excitations are equivalent to those in coronene, however, the excited state energies lay one electronvolt beneath due to delocalization effects.  in the coronene molecule.This analysis was done by looking at the MOs involved in every transition and comparing those orbitals with the ones involved in the coronene excitations.
The graphical representation of these MOs is included in Fig. 3, where the four highest occupied and four lowest virtual MOs of coronene and circumcoronene have been plotted.
These MOs in circumcoronene are pairwise degenerated, similar to the coronene orbitals.
The symmetry of each molecular orbital has been deluded by comparison with respect to the equivalent MO in coronene.The correlation between the MOs in the two systems has been represented by crossing lines in Fig.The circumcoronene absorption spectrum, represented in Fig. 4, has been simulated in the same manner we did for coronene in Fig. 1, where the upper panel shows the static spectrum with the vertical absorption peak allowed by symmetry (S 3 and S 4 ) and including then vibronic effects by computing vertical transitions from 200 initial configurations generated from a Wigner distribution, which reveals a smaller absorption peak at lower energies, corresponding to excite the system to S 2 .The same procedure has been employed to study the out-of-equilibrium dynamics of circumcoronene to observe the effect of the enlargement of the system in the molecular plane.
By exciting the system to the second singlet excited state (1B 3u , S 2 ) the dynamic absorption spectrum has been recovered from the Fourier Transformation of the autocorrelation function and plotted in the bottom panel of Fig. 4 to compare better with respect to the static absorption.The quantum dynamics for this system have been studied using ML-MCTDH on the LVC precomputed potentials calculated at the TD-B3LYP-D3/def2-SVP including all the vibrational degrees of freedom.During these dynamics, the grid sizes were 5.91×10 219 and 2.43×10 50 nuclear configurations were included.Due to its large size, on-the-fly TSH dynamics is out of reach for this system.
We also performed a vibrational mode analysis by exploring the analogous modes that were relevant in coronene.The scissoring mode in coronene (ν 19 ) is equivalent to ν 34 in circumcoronene, which is equally important for the de-excitation dynamics from S 2 in this molecule.In the same way, the totally symmetric breathing ν 24 in coronene is equivalent to ν 28 in circumcoronene, which is also responsible for leading the molecules away from the Frank-Condon region on the S 2 state.In Fig. 5, these four cuts of the adiabatic PES are represented in order to show the clear correlation between these important vibrations in both systems.The right panels correspond to the most important symmetric breathing modes with negligible coupling and large gradients, while the left panels correspond to the most relevant scissoring mode in both systems, exhibiting large diabatic couplings.The diabatic coupling between 1B 2u (S 1 ) and 1B 3u (S 2 ) is plotted in grey.
Although the behavior of both PES is quite similar, it is important to emphasize here that they are displaced to lower energies and that there is a higher density of states in the case of circumcoronene, which can be explained by the size difference and the larger electronic delocalization in this system.
In Figure 6  FIG. 6. Population transfer from ML-MCTDH full-D dynamics starting in every excited state for coronene (left column) and circumcoronene (right column).The order of states has been altered to keep the same order as in Table II so that the dynamics are comparable.

V. CONCLUSIONS
We have benchmarked the level of theory to describe the electronic potential energy surfaces in coronene, finding that TD-B3LYP-D3/def2-SVP provides a good agreement with experimental results.A parametrised LVC model explains the symmetry-forbidden vibronic band at around 3.4 eV.The full dimensional excited state dynamics is driven by breathing and tilting modes.No reduced mode combination, however, leads to the full dimensional result, indicating that most of modes contribute to the dynamics with small couplings.On-the-fly excited state dynamics where electronic quantities are computed along the trajectories show that the contribution of triplet states is negligible and indicates that the conformational molecular changes are correctly described displaced harmonic potentials.As the D 2h assumed symmetry mixes electronic states that would be degenerated at the FC region in the D 6h group, differences in population transfer arise in on-the-fly dynamics starting in the S 2 state when compared to precomputed potentials.The same LVC/TD-B3LYP-D3/def2-SVP methodology was applied for the circumcoronene molecule with 210 nuclear degrees of freedom.Although the absorption spectrum lays 1 eV below the spectrum of coronene, the spectral features are very similar.In fact, the excited state dynamics started from S 2 to S 8 involving twelve coupled electronic states and 210 DOFs agree very well with the coronene dynamics.This result shows that only one electronic structure calculation at the FC in the larger system and a state mapping to the smallest unit -coronene-are needed to describe the early ultrafast dynamics of polyaromatic hydrocarbons and their derivatives.

DATA AVAILABILITY STATEMENT
The benchmark for the vertical energies in coronene, the comparison from spectra computed from 8D dynamics and full dimensionality for coronene, the ML-MCTDH trees which show the nuclear base function distribution among the 102 and 210 degrees of freedom and the detail from the dynamics from every electronic state for both molecules can be found in the Supporting Information file.The frequency calculations at the optimised geometries of both compounds, quantics inputs, operator files and initial conditions for the on-the-fly TSH dynamics are in the comp_data.zipfile.

3g 4 .
FIG. 1. Computed absorption spectrum of coronene.Top panel (a) shows the static spectra computed by vertical excitations from the FC region using TD-DFT (black line) and including vibronic coupling by computing vertical excitations from 200 initial conditions selected from a Wigner distribution.The experimental absorption spectrum is taken from Ref. 41. Lower panel (b) represents the dynamic absorption spectra obtained from the Fourier Transformation (FT) ⟨Ψ(t)|Ψ(t = 0)⟩ of the autocorrelation function of the quantum dynamics starting from 1B 3u (S 2 ) excited state computed with ML-MCTDH and including the 102 normal modes of coronene.
using the three different methods for propagating the nuclear dynamics: (a) on-the-fly TSH classical trajectories, (b) TSH classical trajectories on the LVC potentials and (c) ML-MCTDH quantum dynamics on the same LVC potentials.

FIG. 2 .
FIG. 2. Diabatic population transfer from the 1B 3u (S 2 ) excited state of coronene using three different propagations: (a) On-the-fly Tully Surface Hopping classical trajectories (TSH).The diabatisation has been made just tracing back to the state character of the initial populated state and populations are collected as the square of the quantum amplitudes of the states, (b) TSH classical trajectories using the LVC model for the PES and, (c) ML-MCTDH using also the LVC model potentials.Insets of molecular structures of coronene visited during both the on-the-fly and on LVC potentials using the TSH method are included in the figures (a) and (b) with 8892 and 2213 overlapped structures respectively.
of the coronene molecule.The performance of the LVC approximation has been also tested by inspecting the molecular configurations visited by the nuclear dynamics, represented as overimposed xyz structures in Fig.2right hand side.This figure shows a fairly harmonic behavior of the coronene geometry during the classical trajectories of the nuclei both over the LVC potential and the 'real' ab-initio potential.Another proof-of-principle of the LVC approximation in this system is the fact that the absorption spectrum obtained with the PES is in good agreement with respect to the experimental, represented in the lower panel (b) of Fig.1.As we can see in Table.I, the MOs involved in S 1 and S 2 are the same (with different relative weights), and they are degenerated in D 2h .To unravel why the diabatic population transfer differs, we checked the orbital character at the geometries visited along the on-the-fly dynamics and concluded that the differences arise due to a trivial rotation of the degenerated MOs.

FIG. 3 .
FIG.3.Top view of molecular orbitals of circumcoronone showed in comparison with respect to the corresponding ones in coronene.All these orbitals correspond to π * orbitals asymmetric with respect to the molecular plane.Symmetry labels are included for the coronene system, connected with lines with equivalent MO in circumcoronene, where the symmetry can be extrapolated.

FIG. 4 .
FIG. 4. Computed absorption spectrum of circumcoronene.Top panel (a) shows the static spectra computed by vertical excitations from the ground state optimised geometry using TD-DFT (black line) and including vibronic effects by computing vertical excitations from 200 initial conditions selected from a Wigner distribution, which are shown in colors after convoluting with gaussians with a width of 0.2 eV.Lower panel (b) represents the dynamic absorption spectra obtained from the Fourier Transformation (FT) of the autocorrelation function ⟨Ψ(t)|Ψ(t = 0)⟩ of the quantum dynamics starting from 1B 3u (S 2 ) excited state computed with ML-MCTDH and including the 210 normal modes.

FIG. 5 .
FIG.5.Potential energy cuts along the two leading modes in the dynamics from the excited state 1B 3u (S 2 ) of coronene (top panels) and the corresponding modes in circumcoronene (bottom panels).

TABLE II .
Vertical excitation energies, in eV, and their corresponding oscillator strengths were obtained using TD-B3LYP-D3/def2-SVP level of theory for circumcoronene in D 2h symmetry.
3. It is interesting that the MOs that contribute the most to the bright excited states in coronene, 2B 2u (S 6 ) and 2B 3u (S 7 ), correspond to the same MOs that are involved in the two excited states in cirumcoronene showing a non-zero oscillator strength, 2B 2u (S 3 ) and 2B 3u (S 4 ).
we grouped the excited state dynamics in coronene and circumcoronene from S 2 to S 8 according to their similarities.1B 2u and 1B 3u states behave in a very similar manner in coronene and circumcoronene, transforming population among them in the sub 20-fs timescale.After this time, an equilibrium is reached and both states keep half the population.1B 1g , 1A g , and 2B 1g states have an equivalent initial decay lead by asymmetric in-plane motions.The population of the states is kept constant after 20 fs.When exciting to 2A g in coronene and circumcoronene (S 5 and S 8 respectively) a complete population decay to other states is reached within 10 fs.The bright 2B 2u and 2B 3u states, however, decay at a slower pace, showing lifetimes above 250 fs.