Reactive Force Field (ReaxFF) interatomic potentials can be used in MD simulations to describe bond breaking (van Duin et al., 2001). As illustrated by Eq. 2, ReaxFF includes several contributions to the potential energy, such as two-body energy terms which depend on interatomic distance and electrostatic, torsional and bond order dependent contributions. E system = E bond + E lone pair + E overcoordination + E undercoordination + E valence angle + E penalty + E three body conjugation + E C 2 correction + E triple bond energy + E torsional + E four body conjugation + E hydrogen bond + E vd Waals + E Coulomb (2) Bond order dependent potentials allow description of different sp ⁿ (n = 1, 2, 3) orbital hybridizations of carbon atoms. The bond order is found by analysing the chemical environment and interatomic distances around an atom. For a thorough description of the different energy terms in Eq. 2 we refer to the main paper written by van Duin et al. (2001). Tapered ReaxFF (Furman and Wales, 2020) eliminates discontinuities in the energy, which is inherent to ReaxFF potentials that use bond order terms and torsion angles in the potential energy expressions. The tapered ReaxFF approach is implemented in the Amsterdam Modeling Suite program (Franchini et al., 2021), which is used for the present study. Several ReaxFF potentials are available in this environment for hydrocarbon MD simulations. In our study we use the 2008-C/H/O parametrisation, which is capable of describing C-H bonds and their dissociation (Chenoweth et al., 2008).

Before the MD calculations are started, the indene geometry is optimized. We then wish to simulate the dynamics after absorption of a photon. The relaxation from high lying electronic states to low lying electronic states is typically fast (Allain et al., 1996a; Allain et al., 1996b). In this work we investigate the case where the vibrational relaxation occurs to the ground electronic state before interesting dynamics happens. It is thus assumed that after the absorption of a photon, the photon energy is quickly transferred to the nuclear motions and evolution proceeds on the lowest adiabatic PES, which is described by the ReaxFF potential that we use.

The assumption of fast energy transfer to the nuclei is expected to work better as the size of a system increases. A study of Ghanta et al. (2001) shows rates for non-radiative decay of low lying states of naphthalene and anthracene cations to occur within a few 10th of femtoseconds. In our calculations the energy is initially distributed among the vibrational modes of the molecule in two different ways, either according to a thermal Boltzmann distribution or alternatively the energy is put into a single C-C bond. The simulations are always started from the equilibrium geometry for a non-rotating spatially immobile indene molecule. Thermal initialization of velocities leads to slight overall rotation and translation of the molecule. These two types of motion are subtracted, so initially the total energy is completely in the form of vibrational kinetic energy, and measured from the bottom of the potential well. We report results from the molecular dynamics simulations for various total energies. For the single C-C bond excitation, equal oppositely directed momenta are assigned to the two atoms in the bond while all other atoms are left at rest, thus leaving the center momentum and rotational momentum equal to zero. Since the photo-excitation energy that we wish to simulate overwhelmingly exceeds the experimental thermal energy, neglecting the initial rotational motion is justified (and center of mass motion is as usual not relevant). The MD simulations continue in time until the molecule dissociates. There are however also cases where the dissociation is so slow that we could not follow the time evolution until dissociation occurs.

The MD calculations are performed for up to 10 ns, with a 0.01 fs time-step. Calculations longer than 1 ns may lead to non-physical heat up, thus increasing the total energy of the system by several percent or more. Such trajectories are discarded in the final analysis of the decay rate. In our calculations this problem typically appears if the molecule has considerable amount of kinetic energy thereby allowing angular parameters in the ReaxFF to often correspond to near linear configurations. Decreasing the energy cutoffs for the torsional angles seem to result in smoother changes of the potential energy landscape. We therefore suggest that reparameterising the ReaxFF potentials with lower energy cutoff values might increase the numerical stability and improve the conservation of energy in the MD calculations. Such reparametrization is however outside the scope of this work.

In order to employ RRKM theory to calculate microcanonical reaction rate constants, information about the energetics and rovibrational properties of the reactants and the transition states are required. Second-order Møller–Plesset (MP2) Perturbation Theory (Møller and Plesset, 1934) along with the 6-311+G(d,p) basis set is used to optimize the stationary points on the potential energy surface, i.e., reactants, intermediates, transition states and products. Rotational constants and vibrational frequencies are also obtained at the same level of theory. To compute more reliable relative energies with respect to the reactants, the composite quantum method CBS-QB3 is utilized (Montgomery Jr et al., 2000). Note that the reported energies include harmonic zero-point energies calculated at MP2/6-311+G (d,p) level of theory. Gaussian 16 (Frisch et al., 2016) is employed for all electronic structure calculations.

The obtained optimized geometries (Z-Matrices), vibrational frequencies and rotational constants are presented in the Supplementary Material.

RRKM theory assumes that the molecular system behaves ergodically and that vibrational energy redistribution is faster than reaction (Nordholm and Bäck, 2001; Leitner et al., 2003). These conditions are not always fulfilled, but are typically good approximations. RRKM theory requires sums and densities of states. To find these, we make the separable harmonic oscillator—rigid rotor assumption. It should be noted that at high energies the separable harmonic oscillator approximation might not be good (Nguyen and Barker, 2010).

In RRKM theory, the energy-dependent rate constant, k(E), for a unimolecular reaction can be written (Holbrook et al., 1996; Steinfeld et al., 1999) k E = N ‡ E h ρ E (3) in which h is the Planck constant. N ^‡(E) and ρ(E) are the sum of states for the transition state and the density of states for the reactant, respectively. The precision of the rate constants predicted by RRKM theory relates to how accurately the sums and densities of states are evaluated. In the present study N ^‡(E) and ρ(E) are calculated using a direct count of vibrational quantum states but also using a classical expression to justify comparison to the classical MD simulations. Considering a molecular system which includes s vibrational degrees of freedom with non-degenerate vibrational frequencies (ν _i ), which applies to indene, we may use the following expression to find the sum of states classically (Steinfeld et al., 1999): N E Classical = E s s ! ∏ i = 1 s h ν i (4)

We point out that in this classical expression we choose the bottom of the potential energy well as the reference level, which we also do in the MD simulations. For the quantal RRKM results we choose the zero point level as the reference level for the energy.

In this study, the code Densum in the MultiWell software suite (Barker et al., 2022) was used to perform the direct counts in order to obtain the sums and densities of states for separable harmonic vibrational degrees of freedom. Energy step sizes of 100 cm⁻¹ were used and the densities of states were determined from the number of states in these energy bins. The classical density of states is simply evaluated as the derivative of Eq. 4.

Besides the vibrational motions, external rotations are also considered active in the present work. The rotational states for an asymmetric top like indene cannot be calculated using a simple mathematical expression. Here, we therefore approximate the indene molecule as a prolate symmetric top, for which the rotational constants A, B and C are related as A ≠ B = C. Here we set a rotational constant D = (BC)^1/2 and use it in the rotational energy expression for a prolate symmetric top which becomes (Holbrook et al., 1996; Steinfeld et al., 1999): E rot J , K = D J J + 1 + A − D K 2 (5) where J is the total angular momentum quantum number (the ”J-rotor”) and K is the quantum number for the projection of J on the molecular symmetry axis. The vibrational sums and densities of states are combined with the contribution from the rotational states in the rate constant calculations.

There are three chemically different C-H bonds which can be broken in the five-membered ring of indene. These H elimination reactions have loose transition states. The energies relative to indene for these reaction channels are illustrated in Figure 1. The most favorable hydrogen loss channel, i.e., the one with lowest dissociation energy (D₀), breaks one of the C-H bonds in the CH₂ group. The dissociation energy for this is 3.60 eV at CBS-QB3 level. This value is at least 1.6 eV lower than the energy required for loss of any other H atom in indene. The dissociation energy for the elimination of hydrogen atoms in benzene, where all carbon atoms are sp ² hybridized, has been reported to be close to 5 eV (Kislov et al., 2004).

If there is a concerted H₂ loss channel, it would be expected to occur via a transition state centered on the CH₂ group, similar to the H₂ loss channel in the dissociation of benzene (Kislov et al., 2004). Despite our searches, such an H₂ loss channel was not detected at MP2 level. Several possible pathways for dissociation of indene are shown in Figure 2, with relative energies at CBS-QB3//MP2/6-311+G(d,p) level of theory.

Breaking the C₁-C₂ bond leads to INT1 via TS1 (at 4.92 eV), see Figure 2. INT1 is 4.85 eV less stable than indene and can react to eliminate ethene (C₂H₂) and cyclopropene (c-C₃H₄). INT1 must surmount TS2 (at 6.47 eV) to form C₂H₂ and bicyclo[4.1.0]hepta-1,3,5-triene. It can alternatively pass through TS3 (at 7.11 eV), breaking the C₄-C₅ bond, to form C₃H₄ and benzyne.

The fission of the C₄-C₅ bond in indene produces INT2, with an energy 4.36 eV above indene. At MP2 level, a saddle-point is found in between indene and INT2, but at the CBS-QB3 level the energy estimated for this saddle-point is nearly 0.04 eV lower than INT2. The nature of this transition state is loose and its motion along the reaction coordinate corresponds to the very low vibrational frequency 98 cm⁻¹. INT2 can react to form c-C₃H₄ and benzyne by passing TS5 (at 7.39 eV).

It is possible to form a quite stable isomer of indene, viz. 2-methylphenylacetylene which is 1.29 eV less stable than indene itself. Starting from indene this can happen through cleavage of the C₃-C₄ bond by passing the loose transition state TS6 4.4 eV above indene, to first form the bi-radical INT3. INT3 (at 4.18 eV) is only weakly bound and can via TS7 (at 4.28 eV) easily proceed to form 2-methylphenylacetylene (at 1.29 eV). The newly formed molecule can react to eliminate C₂H (at 7.49 eV) or CH₃ (at 5.95 eV).

The harmonic vibrational frequencies, rotational constants and energy barrier heights found above are used to obtain sums and densities of states and then finally energy-dependent (microcanonical) rate constants, k(E)′s, from Eq. 3.

Figure 3 shows k(E)′s that we have obtained by RRKM using direct count for indene fragmentation. The hydrogen elimination reaction channels are barrierless, i.e., there is no intrinsic saddle point along the reaction coordinate. For such reactions, the transition state can be found by a variational approach that finds the minimum flux (sum of states) along the reaction coordinate. The details of these calculations are discussed in the Supplementary Material.

A hydrogen atom is most easily lost from the CH₂ group due to the low dissociation energy for this process. Thus, in indene such hydrogen atoms require less energy for elimination than hydrogen atoms bonded to sp ² hybridized carbon atoms. Similarly loss of a hydrogen atom from the CH₂ group requires less energy than loss of C₂H₂ or C₃H₄. The observation in Figure 3 that hydrogen atom loss has the largest rate constant is thus expected. An experimental study on the photo dissociation of hydrogen atoms from the CH₂ group in indene obtained a rate constant of 7.4 × 10⁶ s⁻¹ at 193 nm (Yi et al., 1991) corresponding to a photon energy of 6.4 eV. This rate value is close to 1.7 × 10⁷ s⁻¹ predicted by our RRKM calculations.

As seen in Figure 3, C₂H₂ loss is faster than C₃H₄ loss. This is in accordance with the 6.47 eV barrier for C₂H₂ loss, being lower than the barriers for C₃H₄ loss. C₃H₄ loss occurs through two pathways with the barrier heights 7.11 and 7.39 eV respectively. The rate for formation of the molecule 2-methylphenylacetylene (INT4) due to isomerization of indene is also plotted in Figure 3. This molecule can further dissociate to for instance eliminate C₂H and CH₃, but with we have not further investigated these channels.

Figure 4 illustrates energy resolved rate constants obtained by our MD simulations using the reactive force field ReaxFF (2008-C/H/O parametrization). Two different sets of simulations were run in which either the initial energy was randomly distributed among all the vibrational modes or deposited into one specific bond of the indene molecule, viz. the C₁-C₂ bond in Figure 1.

Thermally excited indene fragments through several competing channels involving C-H bond breaking, C-C bond opening and rearrangement of chains of carbon atoms. At the lower energies studied, the dissociation is dominated by the loss of one of the two hydrogen atoms in the CH₂ group. At higher energies, product molecules/radicals like C₂H₂ or C₂H₃ are also observed. Overall, around 500 simulations at random energies between 8 eV and 25 eV were performed. The obtained total rate constants are fitted to Eq. 1. This gives E ₀ = 2.54 eV and k ₀ = 3.73 × 10¹⁴ s⁻¹. The error bars in Figure 4 include 67% of the data points and are calculated using a quantile regression (Koenker et al., 2005).

About 500 simulations were performed in which the energy was deposited into the C₁-C₂ bond. At high energies this leads to ring opening followed by C₂H₂ or C₂H₃ elimination. This rate of dissociation is much larger than the corresponding one for thermally excited indene. The fast dissociation of the molecule occurs before energy is redistributed among the vibrational modes. This is a manifestation of apparent non-RRKM behaviour not described by the statistical approach. At lower energies, the opened ring closes again right away and energy is redistributed within the molecule leading to subsequent dissociation of a hydrogen atom, primarily from the CH₂ group.

Our electronic structure calculations show that it requires 10.2 eV energy to pass the barrier for C₂H₂ elimination. This agrees with the sudden drop in C₂H₂ loss seen in Figure 4 as the energy goes below about 10 eV. Further, the single bond excitation rate constants at energies below 10 eV are in accord with the average dissociation rate constants obtained with thermal excitation. This shows that both direct bond excitation and thermal excitation produce similar results for the lower energies where hydrogen elimination dominates.

Our study does not consider ionization of the molecule. It is expected that for larger molecules the importance of inclusion of ionization diminishes (Zettergren et al., 2021). PAH ions have been studied with on the fly dynamics, see Simon et al. (2017). To include charged species in our MD calculations, reparameterization of ReaxFF for potential terms involving C and H might be necessary in order to achieve quantitative predictions.