Toward Simulation of Fe(II) Low-Spin → High-Spin Photoswitching by Synergistic Spin-Vibronic Dynamics

A new theoretical approach is presented and applied for the simulation of Fe(II) low-spin (LS, singlet, t2g6eg0) → high-spin (HS, quintet, t2g4eg2) photoswitching dynamics of the octahedral model complex [Fe(NCH)6]2+. The utilized synergistic methodology heavily exploits the strengths of complementary electronic structure and spin-vibronic dynamics methods. Specifically, we perform 3D quantum dynamics (QD) and full-dimensional trajectory surface hopping (TSH, in conjunction with a linear vibronic coupling model), with the modes for QD selected by TSH. We follow a hybrid approach which is based on the application of time-dependent density functional theory (TD-DFT) excited-state potential energy surfaces (PESs) and multiconfigurational second-order perturbation theory (CASPT2) spin–orbit couplings (SOCs). Our method delivers accurate singlet–triplet–quintet intersystem crossing (ISC) dynamics, as assessed by comparison to our recent high-level ab initio simulations and related time-resolved experimental data. Furthermore, we investigate the capability of our simulations to identify the location of ISCs. Finally, we assess the approximation of constant SOCs (calculated at the Franck–Condon geometry), whose validity has central importance for the combination of TD-DFT PESs and CASPT2 SOCs. This efficient methodology will have a key role in simulating LS → HS dynamics for more complicated cases, involving higher density of states and varying electronic character, as well as the analysis of ultrafast experiments.


S2 Vibronic Coupling Parameters and Spin-Orbit Couplings
All below parameters were obtained by using DFT/TD-DFT (B3LYP*/TZVP) quantum chemistry for the following diabatic Hamiltonians: i) VCHAM-B3LYP* and ii) LVC-B3LYP*. For i), we introduced the upper S 4 −S 6 singlet states for VCHAM fitting along the two antisymmetric modes ν 13 and ν 14 ; this is required because the 1 T 1g and 1 T 2g manifolds interact along ν 13 and ν 14 . The dynamics simulations were performed with the 1 T 1g singlet states only, as the 1 T 2g states are energetically inaccessible. The utilized spin-orbit coupling (SOC) matrix 1 is given in Tables S9−S14 and the ab initio data for the calculation of SOCs along the ν 13 and ν 15 Fe-N stretching modes are given in ref. 1. Table S2: Zeroth-order coefficients ε (α) , given in eV. The tabulated ε (α) values are the B3LYP*calculated (DFT/TD-DFT) vertical excitation energies at the Franck-Condon (FC) geometry. Mode κ for LVC-B3LYP* (≥ 0.001 eV) triplet states: for LVC-B3LYP* (≥ 0.001 eV) quintet states:

S3 Diabatic and Adiabatic/Spin-Diabatic Potentials
We here analyze the diabatic potentials and coupling strengths for the three different models LVC-B3LYP*, VCHAM-B3LYP*, and VCHAM-CASPT2, which the dynamics simulations are based on. We first consider the antisymmetric Fe-N stretching modes. Figure S1 shows the diabatic PESs along the component of the degenerate antisymmetric Fe-N stretching vibrations (ν 13 or ν 14 ), which accounts for stronger vibronic effects in the diabatic potentials, quantified by the linear on-diagonal κ (α) i constants. As seen in Figure  S1, the singlet and triplet PESs along ν 13 /ν 14 (light blue and green) split, and importantly, the 1 T 1g and 3 T 2g potentials cross in regions that are easily accessible by the wavepacket in the initially excited 1 T 1g states. This is the reason for the fast (∼100−300 fs) singlettriplet ISC. The overall agreement for the PESs shown in Figure S1 is rather good, which supports the adequacy of B3LYP*, as well as the LVC approach for constructing diabatic states. In order to analyze the singlet-triplet ISC, it is also important to consider the magnitude of nonadiabatic couplings between the singlet excited states; this is also shown in Figure S1 for the two antisymmetric Fe-N stretching modes, ν 13 and ν 14 . As seen in the figure, the magnitude of singlet couplings is comparable for the three methods, with the weakest coupling strengths for LVC-B3LYP*. However, even in the case of these weaker couplings (0.01 − 0.05 eV), the LVC-B3LYP* 1 T 1g / 3 T 2g crossings are easily accessible, as reflected by the ∼200−300 fs singlet-triplet ISC timescales obtained by the corresponding TSH and QD simulations (see Figures 4a,c in the article).
In Figure S2, we present the PESs along the totally symmetric Fe-N stretching mode ν 15 , which, as discussed in Section 4.1 in the article, is the key mode for the quintet states (double occupation of e * g orbitals). Similarly to the case of the antisymmetric Fe-N stretching modes, the overall agreement is rather good. There are two aspects, though, for which the DFT/TD-DFT (B3LYP*) and CASPT2 results differ: i) the separation of the two triplet manifolds 3 T 1g and 3 T 2g is roughly twice as large for CASPT2 ( Figure S2, panel c) than for B3LYP* ( Figure S2, panels a and b), and ii) in contrast to CASPT2, for B3LYP*, the perfect triple degeneracy of the quintet 5 T 2g states along ν 15 is not maintained, which is a consequence of the description of the lowest quintet by as single Slater determinant (unrestricted DFT) and the other two by quintet TD-DFT. In addition, for LVC-B3LYP* ( Figure S2, panel a), the HS-LS energy gap along ν 15 is lowered by ∼0.2−0.3 eV. However, as is clear from Figure 4 in the article, these variations in the PESs along ν 15 do not lead to any qualitative differences in the simulated dynamics. We note that nonadiabatic couplings are negligible for mode ν 15 . 2 Figures S3-S5 and S6-S8 show the B3LYP* (DFT/TD-DFT) adiabatic/spin-diabatic PESs of [Fe(NCH) 6 ] 2+ , obtained by the VCHAM and LVC methods, respectively. In the case of VCHAM, both quantum chemistry data points and energy curves obtained by the diagonalization of the diabatic Hamiltonian (fits) are shown. We note again that for VCHAM, modes ν 13 , ν 14 (antisymmetric modes), we had to include three additional singlet excited states (i.e., the three components of the 1 T 2g manifold), as the 1 T 1g -1 T 2g states interact along ν 13 and ν 14 ; this was not necessary for LVC, as the quantum chemistry calculations were only performed in the vicinity of the FC geometry, where no  Figure S9 shows the 3D QD populations dynamics obtained by the VCHAM-B3LYP* and VCHAM-CASPT2 1 models. As is clear from the figure, the triplet IC is affected by the utilized quantum chemistry method. This is explained by the nearly twice as large 3 T 1g -3 T 2g energy gap (see Table S2), which is consistent with the differences observed in Figure S9.

S5 Description of Vibrational Dynamics by TSH and QD
We now assess how the TSH and QD methods describe the nuclear dynamics of photoexcited [Fe(NCH) 6 ] 2+ . For this analysis, we utilize our LVC-B3LYP* Hamiltonian, which allows direct comparison of the TSH and QD results. Figure S10 shows the TSH nuclear dynamics along ν 13 and ν 15 . As seen in Figure S10, panel a, and also the single trajectory in Figure 5a in the article, 1 A 1g → 1 T 1g excitation immediately activates ν 13 . Propagation occurs equally towards positive and negative q 13 values, which is due to the antisymmetric character of ν 13 , which breaks the threefold symmetry of 1 T 1g . This has the consequence that the average displacement cancels out (two components with the same magnitude but opposite sign), and thus ⟨q 13 ⟩ is zero (thick black line in Figure  S10a). This is in line with a bifurcating wavepacket created by nonadiabatic coupling between the 1 T 1g components. The ensemble of trajectories dephases and broadens in ∼100 − 200 fs, which is in agreement with our QD simulations on CASPT2 PESs. 1 From ca. 400 fs, weak, slowly growing oscillations appear in the average ⟨q 13 ⟩, which we attribute to the growing quintet population; this is supported by the asymmetry in the distribution of triplet-quintet hopping geometries projected to ν 13 shown in Figure 7a in the article.
The trajectories along the breathing mode ν 15 are displayed in Figure S10, panel b, which are very different from those along ν 13 , discussed above. The dynamics along ν 15 are characterized by coherent oscillations (breathing vibrations) with a period of ca. 100 fs with rather slow coherence decay. Initially, the oscillations turn at ca. q 15 = −2 and q 15 = 12, which, according to the PESs along ν 15 ( Figure S8) are consistent with vibrations along ν 15 in the singlet ( 1 T 1g ) and triplet ( 3 T 1g , 3 T 2g ) excited states, respectively. At the second oscillation period, a new component appears with an outer turning point at ca. q 15 = 20, which is assigned to the quintet ( 5 T 2g ) states. At the next inner turning point, a third component arises between ca. q 15 = −7 and q 15 = −8, which we assign to the ground state 1 A 1g . These ν 15 structural features are clearly characteristic of the electronic character (shown in Figure S10, panel b), reflecting the e * g occupation ( 1 T 1g , 3 T 1g , 3 T 2g − single, 5 T 2g -double, 1 A 1g -zero occupation). The coherently oscillating average ⟨q 15 ⟩ slowly shifts towards higher q 15 values, which, as seen from Figure S10, panel b, is due to the 5 T 2g component, whose relative weight (i.e., population) gradually increases in line with rising quintet population. We note that oscillations with the same period appear in the corresponding quintet population curve ( Figure 4a in the article, red), evidencing that the triplet-quintet transition is controlled by breathing mode ν 15 . Figure S11 shows the QD vibrational dynamics, characterized by the wavepacket centroids ⟨q i ⟩ and widths dq i , along ν 13 and ν 15 . Here, ⟨q i ⟩ and dq i were calculated by a weighted sum over the electronic states, with the weights being the relative populations. As in the case of TSH ( Figure S10), the QD vibrational dynamics are dominated by the coherently oscillating breathing motion ν 15 , with the same period. The shift towards higher ⟨q 15 ⟩ values due to the population of quintet states is again observed, and loss of coherence is even slower than for the TSH results shown in Figure S10, panel b, in fact, it is negligible. This is reasonable because our TSH model is full dimensional, which thus allows interaction of all modes, while our QD Hamiltonian includes only the three Fe-N stretching modes. Interestingly, our VCHAM-CASPT2 QD simulation with these three modes leads to a significantly faster coherence decay. 3 This means that the in the present case, the influence of vibronic effects, which are stronger for CASPT2, on the coherence decay is stronger than the presence of the bath of vibrational modes. Finally, the wavepacket along the antisymmetric mode ν 13 broadens similarly and on the same timescale as in the case of TSH dynamics. On the other hand, the centroid ⟨q 13 ⟩ does not cancel initially but exhibits small oscillations that are slowly damped, in contrast to the slowly growing oscillations for TSH in Figure S10, panel a. The reason for this difference is that the preparation of the initial state is different: while a single diabatic 1 T 1g component for QD, an adiabatic, i.e., a mixture of 1 T 1g states for TSH. Note, S33 however, that due to the degeneracy of 1 T 1g components within octahedral symmetry, this difference in the excitation scheme should not show up in any physical observable; indeed, the corresponding TSH and QD population dynamics show very good agreement (Figures 4a and 4c in the article). Figure S11: Nuclear dynamics along ν 13 and ν 15 , as obtained from the LVC-B3LYP* QD simulations. ⟨q i ⟩ and dq i denote the centroid and width of the wavepacket, respectively.