Mechanism of Ultrafast Intersystem Crossing in 2‐Nitronaphthalene

Abstract Nitronaphthalene derivatives efficiently populate their electronically excited triplet states upon photoexcitation through ultrafast intersystem crossing (ISC). Despite having been studied extensively by time‐resolved spectroscopy, the reasons behind their ultrafast ISC remain unknown. Herein, we present the first ab initio nonadiabatic molecular dynamics study of a nitronaphthalene derivative, 2‐nitronaphthalene, including singlet and triplet states. We find that there are two distinct ISC reaction pathways involving different electronic states at distinct nuclear configurations. The high ISC efficiency is explained by the very small electronic and nuclear alterations that the chromophore needs to undergo during the singlet–triplet transition in the dominating ISC pathway after initial dynamics in the singlet manifold. The insights gained in this work are expected to shed new light on the photochemistry of other nitro polycyclic aromatic hydrocarbons that exhibit ultrafast intersystem crossing.


S1.1 Computational Details
The excited-state dynamics of 2-nitronaphthalene (2NN) in gas phase were simulated using a local version of the surface hopping including arbitrary couplings (SHARC) approach. [1][2][3] Density functional theory (DFT) and time-dependent DFT (TDDFT) calculations were performed at the PBE0 4,5 /DZP 6 level of theory using the ADF2016 program package 7 and the recently implemented SHARC/ADF interface. 8 Initial Conditions. Before simulating the excited-state dynamics, first a groundstate geometry optimization was performed and normal modes and vibrational frequencies were calculated at the minimum-energy geometry ( Figure S1). Then, 1000 initial conditions were sampled around the optimized geometry using a harmonic Wigner distribution 9 for a temperature of T = 300 K. For each structure, the 10 lowest excited singlet and triplet states were calculated in the Tamm-Dancoff approximation (TDA) using a Becke integration grid and a ZlmFit grid of normal quality. Scalar relativistic effects were included in the zeroth-order regular approximation (ZORA) and spin-orbit couplings (SOCs) were calculated perturbatively.  Figure S1: Ground-state minimum-energy geometry of 2NN in gas phase optimized at the PBE0/DZP level of theory.
Absorption Spectrum. Using the geometries of the initial conditions (see above), the absorption spectrum and the density of states ( Figure S2) were simulated by Gaussian convoluting the oscillator-weighted and unweighted stick spectra, respectively, using a FWHM of 0.1 eV. Good agreement to experimental reference spectra can be obtained for 2NN in methanol (MeOH) using PBE0. 10 Up to our knowledge, there is no experimental absorption spectrum of 2NN in gas phase; in lieu, we show in Figure S2 the experimental absorption spectrum of 2NN recorded in n-heptane, 11 where the good agreement between the positions and intensity ratios of the calculated and experimental absorption bands over a large energy range is apparent. We note that the calculated spectrum does not show the fine structure of the experimental absorption bands. We encountered a similar situation in our study of 2NN in MeOH; there, however, we could explain the appearance of a shoulder in the third experimental absorption band by the presence of different electronic states, even though the shoulder was not visible in our calculated spectrum. 10 Figure S2). Accordingly, 51 trajectories were started in the S 1 state and 54 trajectories were started in the S 2 state. Based on the density-of-states at the initial conditions, three singlet states (S 0 , S 1 , and S 2 ) and six triplet states (T 1 -T 6 ) were included in the SHARC simulations and two additional singlet and triplet states were calculated as inactive states. To illustrate the characters of the excited-states included in the dynamics simulation, we show their corresponding states calculated at the FC geometry in Table S1. As mentioned in the main paper,

S1.2 Problems Encountered in the Simulations
As mentioned above, a small number of the calculated trajectories had to be excluded for the analysis. Five trajectories from the T = 300 K ensemble were excluded as they crashed prematurely (t < 332 fs); all other trajectories could be propagated until t = 500 fs. An additional trajectory was excluded because it hopped to a (previously) inactive state * . The remaining 99 trajectories were considered for the analysis.
Additionally, we comment on two critical cases encountered for a small number of trajectories. The first is given by trajectories that exhibit a "unusually" large fluctuation of the total energy which occured in some cases where the S 0 and S 1 were nearly degenerate. Then the total energy could increase momentarily by up to 0.2-0.5 eV before it fell back to the previous range after the states separate in energy again. An example of this situation is shown in Figure S3(a), and one may speculate that it is due to the single-reference nature of TD-DFT that can only give an incomplete description of the multi-configurational nature that the ground state S 0 ingests when it is nearly degenerate with the S 1 . † We note, however, that this is not the general behavior when the S 0 and the S 1 are in close in energy, as can be seen in the exemplary situation depicted in Figure S3(b) where the total energies remains almost constant in the near-degeneracy region. In total, we found the large energy fluctuation due to a S 1 /S 0 near-degeneracy for 8 trajectories. We did not exclude these trajectories from the analysis for the following reasons. First, the DFT/TDDFT calculations do converge in this region of the potential energy surface (PES), and as the total energy goes back to the range outside of this region, we can conclude that the the trajectory can safely leave this region. And second, one may speculate that a trajectory entering the region of near-degeneracy should hop to the ground state; however, absence of this event does not indicate an error as it is in general possible for a trajectory to hop to * the SHARC program is able to follow the character of the MCH states by monitoring wave function overlaps 17 of the calculated states † the trajectory in Figure S3 the ground state without showing the fluctuation in the total energy, see Figure S3(b).
Potential problems that TDDFT surface hopping can exhibit in the region of S 1 /S 0 near degeneracies are discussed by Tavernelli and co-workers in Ref. 18.
The second critical case is provided by situations where hops occur between states with a large energy gap (> 1 eV). An example is shown in Figure S4. There, a trajectory is in the S 1 state until t = 204.5 fs, when it hops to the S 2 state which lies 1.63 eV higher in energy. In order to conserve the total energy of the trajectory, this sudden jump in potential energy is compensated by a rescaling of the momenta of the system decreasing the kinetic energy, as can also be seen in Figure S4. The rescaling is performed by changing the momenta of all atoms in the same manner, i.e., via multiplication by the same factor so that the total energy is conserved. The probability for the large energy hop is very small amounting only to P ∼ 0.001, however, this hop still occured. ‡ For such an event to be statistically significant, one should propagate at least N = 1/P trajectories -so that the event occurs stochastically once. Such low- ‡ to decide whether or not a hop occurs, a random number ∈ [0, 1] is sampled and compared to the hopping probability; when the probability is larger than the random number, the trajectory hops to a different state probability hops occurred 12 times. Despite this is a considerably large number for the comparatively small numbers of trajectories (N = 99), we believe that the occurrence of the low-probability hops does not have a significant impact on our results. All lowprobabilities hops occurred either from the S 1 or the T 1 to higher-lying excited singlet and triplet states, respectively. Before the low-probability hop, most of the trajectories had been for a considerable simulation time in the S 1 or T 1 state, and this (initial) state was quickly retrieved within the remaining simulation time (average retrieval time τ ret = 18 fs). Thus, we can safely include these trajectories with low-probability hops in our analysis.

S11
S2 Mechanism of Excited-State Dynamics S2.1 Major Reaction Channels Figure S5 shows the time evolution of the populations of the MCH states of 2NN in gas phase. Figure 2 in the main paper is similar, but shows only the sum of the excited triplet state populations T n (n = 2-6) and the lowest-energy triplet state T 1 , allowing for a general reaction mechanism, where the population is transferred from the S 1 to some higher-lying triplet state T n before subsequent deactivation to T 1 . This can be   Table S3: Rate constants k i (in fs −1 ) as well as corresponding mean time constants τ i = 1/k i and errors ∆τ i (in fs) calculated using the bootstrap method for the time evolution of the MCH populations according to the mechanism shown in eq. (1).
Following eq. (1), we have calculated functions for the time evolution of the S 2 /S 1 /T n /T 1 populations using a nonlinear least-squares fit, which are shown in Figure 2 in the main paper. The rate constants k S /k ISC /k T and corresponding time constants τ i = 1/k i obtained in these fits are shown in Table S3.

S2.2 Minor Reaction Channels
The mechanistic model of eq. (1) only included the major reaction channels in the first 500 fs of the excited-state dynamics of 2NN. The model does not include direct relaxation from the S 2 to the triplet states, the relaxation to the ground state S 0 , and the individual relaxation processes within the manifold of the higher-lying triplet states T n (n = 2-6). The direct S 2 → T n relaxation is a minor channel populating the triplet states compared to the major route S 2 → S 1 → T n (see hops in Table S2).
Based on our simulations, for this minor channel the time constant § would amount to τ = 1099 ± 8540 fs, i.e., possessing a huge error due to the small number of trajectories.
Since this time constant is not reliable, we can only conclude that the S 2 → T n process is much slower that the competing S 2 → S 1 process. Likewise, the relaxation to the ground state S 0 from the S 1 state is also not well represented by a sufficient number of trajectories and a time constant is meaningless, as it occurs on a slower time scale than our simulation time (τ S 1 →S 0 = 10205 ± 10629 fs. § ). Also note that due to the small number of trajectories between the individual T n states, it is not possible to calculate individual reaction rates without problems of overfitting. Instead, we resort to describe intersystem crossing from (mainly) the S 1 state to the higher-lying triplet states T n (n = 2-6).

S2.3 ISC Pathways
In the main paper, we stated that ISC occurs via two different pathways, a 1 nπ * → 3 ππ * route and a 1 ππ * → 3 nπ * route. This was established after manual inspection of the natural transition orbitals of the singlet and triplet excited states mediating the ISC process for a small fraction of all ISC hops occurring in the dynamics. ¶ To subsequently analyze all ISC hopping events, we calculated the mean atomic electron/hole § arithmetic mean of 100 bootstrap copies with error margins ¶ the small fraction consisted of ca. 30 ISC hops out of a total of 418 ISC hops occurring in all trajectories difference population of the oxygen atoms of the nitro group for all singlet states using the TheoDORE program package. [20][21][22][23] If the difference was negative, i.e., electron density at the oxygen is lower in the excited state than in the ground state, the state was characterized as the 1 nπ * state, whereas if the difference was positive, it was characterized as the 1 ππ * state (see Figure 3(b) in the main paper). The character of the receiver triplet state for the ISC follows the El-Sayed's rules, 24 which we verified for the fraction of ISC hops that we inspected individually.
Having characterized the ISC hops as either 1 nπ * → 3 ππ * or 1 ππ * → 3 nπ * , we have calulated the mean atomic hole/electron difference populations for both sets of states which are shown in Figure 4 in the main paper. Alongside we also show the main natural transition orbitals (NTOs) involved in the excited states at the FC geometry, which resembled the main NTOs describing the excited states that we inspected individually.
Note that the actual NTOs in the transition differ from those depicted in Figure 4 to a certain degree, e.g., because of the distortion of the structures due to admixture with other orbitals. As such, the NTOs in Figure 4 in the main paper are shown only for illustration purposes. As the structures are taken from the dynamics simulations, naturally they differ from each other and from the optimized FC geometry. Thus it is not straightforward to define, e.g., average NTOs that describe different (although similar) excitations at different points of the trajectories. In contrast, the hole/electron difference populations are atomic properties that can easily be averaged over different structures and are therefore more convenient to analyze. Note that, unfortunately, as atomic properties, the hole/electron populations do not contain any "directional" information, i.e., they do not allow a direct discrimination between nπ * and ππ * states if the n and π orbitals are located at the same atoms -see, e.g., the similar populations for the T LE (π π * ) and T LE (nπ * ) states in Figure 4 in the main paper. However, having verified El-Sayed's rules at a sizeable number of excited states in the ISC hops, we feel assured to apply the rule generally, allowing us to safely differentiate all 1 nπ * → 3 ππ * from 1 ππ * → 3 nπ * hops.
Based on the discrimination of the two ISC pathways, we determined their importance in the excited-state dynamics. Their relative contribution was shown in Figure 4 in the main paper, which was calculated as the fraction of trajectories that entered the triplet manifold via the respective pathway and stayed in the a triplet state after a simulation time of 500 fs. These numbers are collected in Table S4. The actual number of trajectories ending up in a triplet state via the minor pathway is small, i.e., 4 trajectories, raising the question on the statistical significance (see above). Since the ratio for the contributions of the minor and major pathway is similar within the number of trajectories (ca. 9 %) and the number of hops (ca. 12 %), we believe that the relative contributions of the minor and major pathways determined by the number of trajectories ending in a triplet represent an acceptable estimate.
The total numbers of hops from singlet to triplet states are much larger than the number of trajectories ending up in a triplet state as it occurs quite frequently that a trajectory hops from a singlet to triplet state in one time step only to hop back to a singlet state at the next time step. Thus, a trajectory needs in average ca. 10 "attempts" to hop from a singlet to a triplet state in order to remain in a triplet state.

S2.4 Spin-Orbit Couplings in the ISC Pathways
In this section we include a brief discussion of the SOCs during the dynamics. For this, we first show the time evolution of two trajectories undergoing ISC via the major and minor pathway ( Figure S6) along with the size of the SOCs between the S 1 and the lowest-lying triplet states T n . The SOCs were calculated as the vector sum of the SOC matrix element between the S 1 singlet state and all three components of the respective triplet state, i.e., The trajectories in Figure S6 were selected as they show a S 1 -T n hop very early in the simulation time and afterwards stay for the remaining simulation in the triplet manifold; we only show an interval of 20 fs around the hopping event.
Hop to T 4 Figure S6: Time evolution of spin-orbit couplings between the S 1 and the lowestlying triplet states T n for a trajectory undergoing ISC via the (a) major pathway or (b) minor pathway. The grey areas denote the time steps between which ISC hops occur.  Figure 4).
The trajectory undergoing ISC via the minor pathway [ Figure S6(b)] initially is in the S 1 which is the S CT (ππ * ) state. Although there are several triplet states close in energy to this singlet state, the SOCs are considerably smaller (< 10 cm −1 ) than for the trajectory undergoing ISC via the major pathway. Still, eventually at 30.5 fs, the trajectory hops from the S 1 state to the T 4 state, which is the T LE (nπ * ) state. The SOC during this hop is 9 cm −1 which is again of similar size as the average value of SOCs reported for the minor pathway in the main paper, i.e., SOC ∼ 8 cm −1 (see Figure 4). Shortly, after the trajectory has hopped into the triplet manifold, the MCH singlet states change in character and the S LE (nπ * ) becomes the S 1 at 33 fs. This state possesses larger SOCs with the low-lying triplet states and explains the sudden increase in S 1 -T n -SOCs shown in Figure S7.
The size of the SOCs is related to the probability of singlet-to-triplet ISC occuring via the minor or major pathway, but it is also related to the probability for the trajectory to hop back from the triplet to the singlet state. This usually occurs within a few time steps after the singlet-triplet hop, where the trajectory is still in the same region of the potential energy surface, meaning that the energy gaps and SOC between the states are approximately the same and, thus, so is the hopping probability. To demonstrate this, we show in Figure S7 the time evolution of a trajectory that hops back and forth six times between singlet and triplet states during the last 200 fs of the 500 fs simulation time. the end our simulation. This example demonstrates that a large SOC is not sufficient for ISC, but the trajectory must also be able to quickly escape the singlet-triplet hopping region while being in the triplet state.

S2.5 Initial Dynamics in the Singlet States
In the main paper, we analyze the S CT (ππ * ) → S LE (nπ * ) dynamics. To do so, we discriminate the states as either S CT or S LE by their dipole moment µ using µ > 8 D

S20
for S CT and µ < 8 D for S LE , respectively. The arbitrary value of 8 D is chosen because it lies in between the initial range of dipole moments (ca. 10-12 D) and the value towards which the ensembles converges (µ = 4 D, see Figure 5(a) in the main paper). Using µ = 8 D as a criterion to distinguish the states S CT and S LE , we obtained the time constant of τ = 81 fs for the population transfer of S CT → S LE . Despite the arbitrariness of µ = 8 D, we note that slightly different reasonable choices of µ give similar values for τ , e.g., ranging from 90 fs (µ = 6 D) to 78 fs (µ = 9 D). We consider these values reasonable as they also lie between the dipole moments the S CT (ππ * ) and S LE (nπ * ) states at their optimized geometries (Table S5). Table S5: Coordinates inÅ of the optimized geometries of the S CT (ππ * ) and S LE (nπ * ) states as well as their energies and dipole moments. In the above analysis, we assume that the initial dynamics takes place within the singlet states and thus all the population is in the singlet states. Figure S8 shows the average dipole moment for the singlet, triplet and all the states juxtaposed. As it can be seen, the average dipole moment of all trajectories µ all (t) (green line) resembles very closely to that of the singlet states µ S (t) (blue line), even at longer times when the triplet states already have gained substantial amount of population.   Figure S9 shows the normal mode analysis (NMA) 25,26 for different sets of geometries: all geometries in the initial conditions (labelled as "Init"), the hopping geometries corresponding to the major/minor ISC pathway (labelled as "Major"/"Minor"), and the geometries where the active state is a singlet state with a dipole moment µ < 8 D/µ > 8 D (labelled as "S LE (nπ * )"/"S CT (ππ * )"). The corresponding vibrational frequencies are listed in Table S6 Figure S9: Normal mode analysis of the initial conditions, all hopping geometries in the major and minor pathway, and all geometries in the 1 nπ * and 1 ππ * states. Normal modes 1-6 refer to translation and rotation and are not depicted here. See Table S6 for the vibrational frequencies of all normal modes.

S2.7 Nuclear Motion
In the main paper we discuss the nuclear motion in terms of three internal coordinates that are responsible to drive the system towards intersystem crossing (recall Figure 6).
Trajectories in the singlet or triplet state are discriminated with the spin expectation value of S 2 ; those in either the S LE (nπ * ) or S CT (ππ * ) state are discriminated with the dipole moment criterion as explained above. time; this is e.g., the case at early times for the triplet states and at later times for the S CT (ππ * ) state. Accordingly, caution is required to interpret the behavior of the average internal coordinates at these points. For completeness, the evolution of the torsional angle τ CCNO and its standard deviations are also shown in Figure S11, as the torsion of the nitro group has been discussed in Refs. 12, 27, 28. As it can be seen, the average of τ CCNO is never zero, but it is small, in agreement with the fact that no dissociation is observed for 2NN.
Note that the curves shown in both Figures 6(b-d) in the main paper and S11 are not pure averages but Bezier fits calculated with GNUPLOT. We have excluded 10 data points in the beginning of the simulation from the Bezier fit of the average internal coordinates of the triplet trajectories as there were no trajectories in a triplet state yet. In both Figures 6 and S11 we use the average of both bond distances r NO and the four torsional angles τ CCNO that can be derived from the molecular structure. Figure S11: Time evolution of the averages and standard deviations of internal coordinates of the trajectories in S LE (nπ * ) (blue), S CT (ππ * ) (green) or any triplet state (red), averages and standard deviation of internal coordinates at the hopping geometries of the major (violet) and minor (orange) ISC pathways, and reference internal coordinates at the FC geometry (gray).