Elucidating the reaction kernel and probing the effect of anharmonicity in the ring-closing reaction of fulgide single crystals

Abstract

Chemistry involves dynamics that transform chemical structures from one form to another. However, among the vast milieu of quantum vibrations in a molecule, it boils down to a few key motions that drive the system across the transition state. It is the anharmonicity at the transition state or barrier-crossing region that couples normal modes, leading to localized motions and reduced dimensionality. The interplay of strongly anharmonic local modes collectively drives the system across the barrier-crossing region, forming a photoproduct. Ultrafast broadband transient absorption spectroscopy has revealed the effect of reduced dimensionality in a prototypical ring-closing reaction in fulgide single crystals. The relatively large anharmonicity at the reactive crossing and the strong reaction forces experienced during the chemical transformation provide a significant driving force for the vibrational modes, revealing a new mechanism of coherent vibrational energy transfer between molecular modes. This effect is observed as a non-impulsive growth of modulation in the amplitude of an 80 cm⁻¹ mode coupled to the reaction coordinate. Our study sheds light on the lattice-coupled reaction dynamics owing to specific system–bath interactions and provides new insight into utilizing lattice alignment for chemical transformation in a solid-state crystalline environment.

---

Introduction

Determining the molecular modes in a barrier-crossing process has always been challenging. Out of the 3N-6 vibrational degrees of freedom for a molecule with N atoms, there happen to be only a few key reactive modes that define the reaction kernel at the barrier-crossing region.^1–3 The “magic” of chemical transformation is the reduction of the defining moment of chemistry, involving nuclear passage over a barrier separating one stable structure (reactant) from another (product), to only a few key reactive modes. Observation of these key reactive modes gives a detailed understanding of the reaction coordinate and the competition among the myriads of inconsequential quantum vibrations or orthogonal coordinates to the reactive passage. These reaction modes effectively define the reaction mechanism, which reproducibly executes the same chemistry on related molecular systems, which in turn makes chemistry a transferable concept. This collapse in reactive nuclear degrees of freedom arises from the highly anharmonic motions in the barrier-crossing region that strongly couple various normal modes, leading to the localized motions defining the reaction coordinate.¹

Several time-resolved ultrafast studies have been done over the past two decades, but very few have been successful in resolving the full complement of reaction modes for a given class of reactions.^4–14 Most studies of chemical reaction dynamics have been conducted in the solution-state, where there is a broad distribution of orientations for the reactant molecules, along with a wide array of possible solvent configurations in relation to the system–bath coupling of the reaction coordinate. The inherent broad distribution of initial conditions causes uncertainties in the phase correlations of the molecular vibrations involved in the reactive crossing, masking the full extent of the reaction kernel in directing motions through the reactive crossing point.^6,15,16 On the other hand, the crystalline state gives us the opportunity to carry out reactions in a solvent-free environment with well-defined initial conditions. The crystal lattice, along with its constrained environment, can promote specific system–bath interactions owing to the spatial arrangement of the reactants within the lattice, which offers the opportunity to control, to a certain degree, the anharmonic coupling and the length scales involved in nuclear reorganization along a reaction coordinate.^{6–8,10–12,16–22} However, solid-state reactions are subject to issues pertaining to the reversibility of the initial excitation condition to enable sufficient sampling to spectroscopically resolve the reaction dynamics. This limitation can be bypassed using sufficiently low repetition rates for the excitation laser and using photochromic systems where the final product state has a well-separated absorption band that can be optically pumped back to regenerate the reactant state.²³ Here, we use the model photochromic fulgide system,^24,25 Aberchrome 540, to provide a single crystal environment to resolve the reaction kernel, lattice interactions, directed bond-making, and subsequent structural reorganization involved in this ultrafast electrocyclic ring-closing reaction.

Photochromism involves reversible changes in the chemical, structural, and physical properties of a molecule upon absorption of light.^26,27 This reversible feature makes these systems suitable candidates for optical switches, memory devices, photoactive molecular machines, tools in biological imaging, and so on.^28–33 Ground-breaking work by Irie et al. has demonstrated the role of photochromism in chemical systems.^34,35 Even though most technological applications require a solid-state environment, there have been very few ultrafast studies carried out in the solid-state, predominantly due to reversibility issues.^36,37 Aberchrome 540 (herein referred to as fulgide) undergoes a ring-closing reaction upon absorption of ultraviolet (UV) light around 340 nm to form a closed-ring photoproduct with an absorption maximum of 520 nm (see Fig. 1a and c). The reverse reaction of ring opening can be initiated with visible light for a considerable number of cycles with minimum fatigue resistance. Absorption of UV light also causes isomerization, which acts as a competing pathway in the ring-closing reaction.^38,39 Previous ultrafast studies in the solution-state have shown that the fulgide ring-closing reaction occurs on the timescale of 100 fs, and a similar study on the fulgimide derivative showed a slightly slower sub-picosecond timescale with a competing isomerization process in the former case.^40–42 However, in a crystalline environment, because of confinement, the isomerization pathway is inhibited.^43–45 Moreover, the electrostatic environment in the crystal is vastly different from the solution-state, strongly influencing the system–bath interactions. The influence of crystal lattice on the reaction dynamics has been previously shown, and it is known to be influenced by the spatial arrangement of the reactant.⁶ Only a few ultrafast studies have been performed on electrocyclic photochromic molecules in their crystalline state, and specifically, no previous work on fulgide single crystals.^11,12,23,46 This work addresses this shortfall to specifically uncover the mechanism of the photoinduced ring-closing reaction in the fulgide single crystal, i.e., resolve the specific reaction modes directing motions through the reactive crossing point from the excited state reactant surface to the product surface. This passage is defined by the conical intersection (CI) between the open and the closed-ring potential energy surfaces (PES), which is a central feature of a large body of photochemical processes involving structural transitions.⁴⁷

Fig. 1 (a) Experimental absorption spectrum of the fulgide single crystal. The spectrum of the open-ring form is multiplied by a factor of three. (b) Visualization of the fulgide single crystal unit cell along the crystallographic b axis, representing the optical beam (340 nm) axis for transient absorption spectroscopy. (c) Schematic representation of the photochromic reaction in fulgide transitioning from the open-ring (upper structure) to the closed-ring form (lower structure), which involves a 1.8 Å change in the Cx–Cy distance from the crystal structure⁴⁴ (the Cx–Cy atoms are highlighted in yellow; carbon atoms are depicted as black balls, oxygen atoms as red balls, and hydrogen atoms as pink balls). (d) Calculated HOTO (highest occupied transition orbital) and LUTO (lowest unoccupied transition orbital) for the open-ring isomer. There is a change in electron density in the excited state between the labile carbon atoms that provides the force leading to a biased sampling towards bond formation in the excited state, clearly observable in the HOTO to LUTO transition.

Herein, we present transient absorption (TA) measurements of the ring-closing reaction in a fulgide single crystal (see Fig. 1b) performed with a temporal resolution of 35 fs (σ of the IRF, Instrument Response Function, SI Section 1.1, eqn (5)). Through a combination of global analysis methodology, time-frequency analysis, and Fourier filtering methods, we have resolved the reaction kernel of the ring-closing reaction in the fulgide single crystal along the CI in a barrierless transition. The present work has also revealed new details in the reaction mechanism, which involve a non-impulsive growth of vibration amplitude on an ultrafast timescale that light up as coherent photoinduced motions directing the chemical transformation to the overall relaxed photoproduct state. This new insight is observable in a solid crystalline environment where the initial conditions are well-defined with specific system–bath interactions, highlighting the highly anharmonic nature of the product surface near the reactive crossing or CI. Complementary ab initio calculations support the multidimensional energy landscape of the ring-closing reaction in the fulgide single crystal.

Results and discussions

Ultrafast spectroscopy on the fulgide single crystal

The open-ring fulgide molecule has a strong absorption band in the UV region, as shown in Fig. 1a. Time-dependent Density Functional Theory (TDDFT) calculations were performed to gain insight into the nature of the transition at 350 nm, as shown in the calculated absorption spectrum (see Fig. S1, SI). The absorption band at 350 nm is a HOTO (Highest Occupied Transition Orbital) to LUTO (Lowest Unoccupied Transition Orbital) transition as shown in Fig. 1d. Upon UV photoexcitation of fulgides in the single-crystal TA experiments (see Material and method section for experimental details), we see a broad excited state absorption (ESA) feature spanning the probe spectrum from 430 nm to 680 nm in Fig. 2a and b. We performed a global analysis procedure on the TA data to understand the underlying dynamics. We fitted the TA data into a sequential model (A→ B→ C→ D→ E→), where each compartment corresponds to a certain spectral species (SI, Section 1.1; see Fig. S2 for SVD filtering, and Fig. S3 and S4a for the temporal evolution of the different spectral species). Fig. 2c shows the Evolution-Associated Difference Spectra (EADS) for the time windows depicted in Fig. 2a and b. We found spectral components associated with 68 fs, 135 fs, 1.13 ps, and 3 ps time constants, along with a non-decaying component. On closer inspection, we find the EADS feature associated with the 135 fs (Fig. 2c) time constant blue-shifted relative to the 68 fs EADS feature. This blue shift highlights the relaxation of the initially prepared photoexcited state (species A) to an intermediate state (species B) in 68 fs, followed by the transformation of the intermediate state to species C in 135 fs. Species C spectrally corresponds to the photoproduct absorption band that appears in 135 fs (Fig. 2c; see Fig. S4b for the spectral feature of species C), indicating the ultrafast formation of the closed-ring photoproduct.

Fig. 2 (a) and (b) TA data of the fulgide single crystal for two different time windows with step sizes of 20 fs and 50 fs, respectively, for the time delays (note: the saturated line at 530 nm represents the background signal from the reversion pulse; see Fig. S5, SI). (c) Evolution Associated Difference Spectra (EADS) were obtained after sequential analysis for the time windows covering (a) and (b). (d) The calculated potential energy for the ground state S₀ and the excited state S₁ along the reaction path of the fulgide ring-closing reaction. The top panel illustrates the evolution of the Cx–Cy distance (highlighted in Fig. 1c) along the reaction pathway. The bottom x-axis indicates the root-mean-square distance (RMSD) between the molecular geometries in the reaction pathway and the reactant, which is the ground state equilibrium geometry. The highlighted line illustrates the presence of a local minimum in the excited state surface that would create an intermediate state.

We also calculated the PES along the Cx–Cy reaction path for the ground S₀ and excited S₁ states, as shown in Fig. 2d (see Materials and methods for computation details). From the calculations and the global analysis procedure, we associate the ∼70 fs process with the relaxation of the photoexcited molecules in the S₁ PES as they evolve out of the Franck–Condon (FC) region to form the intermediate state. During this evolution, the Cx–Cy distance changes from 3.4 Å (FC region) to 2.19 Å (intermediate state). This is followed by the crossing of the S₁–S₀ CI to form the closed-ring photoproduct (species C; Fig. S4b, SI) in 135 fs. The theoretical calculations (Fig. 2d) predict an overshoot in contraction followed by a gradual increase in the Cx–Cy distance from 2.19 Å to 2.24 Å as the system evolves from the intermediate state to the crossing point (S₁–S₀ CI). The transformation to the closed-ring photoproduct leads to a Cx–Cy distance of 1.76 Å, highlighting the formation of a carbon–carbon bond (see Fig. 1c and 2d). We found an additional 1.13 ps EADS component, and we assign this to the structural reorganization in the photoproduct PES after the Cx–Cy bond formation. Since the ring-closing reaction occurs in 135 fs, which involves a substantial change in the Cx–Cy bond distance (3.4 Å to 2.24 Å at the CI) in a constrained crystalline environment, the closed-ring structure is not fully relaxed to the global minimum of the photoproduct. Hence, it must undergo structural reorganization to reach the potential energy minimum of the photoproduct geometry. This relaxation process is followed by vibrational cooling of the photoproduct molecules, as seen in the narrowing of the photoproduct band in 3 ps (Fig. 2c and S4c, SI). From our EADS analysis, we also obtained a long component, which indicates non-decaying photoexcited molecules within the time window of our measurement. However, they all relax to their initial ground state before the arrival of the next actinic pump pulse in 8 ms, determined experimentally by a repetition rate dependence and incorporated into the experimental protocol (Fig. S5, SI).

Vibrational dynamics of fulgide molecules in the single crystal environment

The example of oscillatory beating patterns is shown in Fig. S6a (SI), which is a 2D residual map from the sequential global analysis procedure (the FFT spectrum over the wavelength range is shown in Fig. S6b, SI). We have identified three wavelengths for the vibrational analysis; they are 465 nm, 520 nm, and 560 nm, respectively (as depicted in Fig. 3a–c). The coherent oscillations in the TA signal represent the vibrational dynamics coupled to the optical absorption band in which the excited state PES provides a displacive force along these vibrational coordinates.^48–51 These oscillations include reactive nuclear dynamics directly participating in the ring-closing reaction as well as the FC active ground/excited state modes populated through the resonant impulsive Raman process.^52,53 To identify the vibrational modes participating in the ring-closing reaction, we applied time-frequency analysis to isolate frequency components of interest^4,6 (see SI, Section 1.2).

Fig. 3 The panels (a)–(c) represent the residuals obtained from the global analysis (sequential) procedure along with their corresponding Fast Fourier Transform (FFT) in the bottom plots for 465 nm, 520 nm, and 560 nm, respectively. Panel (d) shows the decay of the reactive modes (234 cm⁻¹ and 296 cm⁻¹) and tuning mode (370 cm⁻¹) obtained from fitting the Fourier filtered TA residual data (explained in SI Section 1.2). Decay time τ is shown in the inset.

Our analysis shows that the key reactive modes correspond to 296 cm⁻¹ and 234 cm⁻¹, decaying in 163 fs and 130 fs, respectively (Fig. 3d). The vibrational modes that decay on a similar timescale as the ring-closing reaction (∼135 fs) are assigned as the key reactive modes driving the ring-closing reaction (see SI, Section 2). Among all the vibrational degrees of freedom in fulgide, these key reactive modes define the reaction kernel or the epicenter for the reaction forces at the CI between the photoproduct and the excited state surfaces, which results in the enormous reduction of dimensionality of the reactive space to these few key modes. Thus, the CI is not a singularity but serves as a high-dimensional seam that acts as a gateway or a portal through which molecules can either reach the product or the ground electronic surface.⁵⁴ In contrast to the solution-state TA measurement, two vibrational modes (59 cm⁻¹ and 120 cm⁻¹) were observed and interpreted as excited state vibrational coherences.⁴⁰ However, the analysis fails to characterize the nature of the vibrations and their projection onto the reaction coordinate.

In addition to the aforementioned vibrations, we also observe the presence of a 370 cm⁻¹ mode decaying on a much faster timescale of 69 fs (see Fig. 3d). From our transient absorption global analysis (Fig. 2c and d), we have already noted a 70 fs component, which defines the evolution of the system from the FC region to the intermediate state. On this basis, the 370 cm⁻¹ mode is assigned as the tuning mode, which reduces the energy gap between the electronic states participating in crossing and brings the system to the intermediate state.^55,56 In this context, we have also calculated the excited state Raman spectrum at the FC region and found strong Raman activity in the 1600 cm⁻¹ frequency region, which involves the stretching of the ethylenic double bonds participating in the ring-closing reaction (Fig. S16a and b, SI). Given that the C=C stretches in the electrocyclic moiety are strongly perturbed by the change of electron density during the ring-closing event (see Fig. 1d), these ethylenic stretches with bond elongation and softening can modulate the excited state PES. In analogy with the photoisomerization of vision,^4,55 the bond elongation of the C=C stretch modes would lead to larger torsional motion in bringing the labile Cx–Cy carbon atoms into wavefunction overlap. From the excited state Raman calculation, we assigned the nature of the 370 cm⁻¹ mode (see Fig. S16c, SI). The vectorial representation of the calculated vibration (421 cm⁻¹) in Fig. S16c shows the movement of the highlighted Cx–Cy atoms in the same direction, which is important for the ring-closing reaction. This motion directs the reduction of the Cx–Cy bond distance (3.4 Å to 2.19 Å, Fig. 2d) during the initial evolution from the FC region to the intermediate state. Thus, the low-frequency 370 cm⁻¹ mode brings the system to the intermediate state, and the reactive modes (296 cm⁻¹ and 234 cm⁻¹) allow the system to cross the CI to form the closed-ring photoproduct.

In addition to the reactive modes, we also observe vibrational modes decaying on a sub-picosecond timescale (Fig. 4a). These vibrations highlight the structural reorganization occurring in the photoproduct PES after the ring-closing reaction. The intrinsic constraint of the crystalline environment ensures a well-defined starting point in the reaction trajectory. Moreover, the solvent-free environment in the crystal lattice ensures that the molecular vibrations are not statistically washed out, allowing us to track vibrations without the loss of phase correlation. This enabled us to explore the reaction mechanism beyond the initial reactive crossing point. We observe a series of frequencies ranging from 14 cm⁻¹ to 204 cm⁻¹ (listed in Tables S1–S5 and Fig. 4a) that participate in the relaxation of the closed-ring photoproduct to its lowest potential energy geometry in a sub-picosecond timescale. Since the decay time is <1 ps and is greater than the timescale of the ring-closing reaction (135 fs), these modes act as assistive modes or dissipative channels for the system to relax to the photoproduct minimum. These modes are not part of the reaction coordinate per se but collectively define the reorganization coordinate and provide dissipation to make the crossing irreversible, highlighting the multidimensional reaction phase space associated with the fulgide ring-closing reaction in a single-crystal environment. A similar multidimensional reactive space was hypothesized for the fulgimide derivative but remained unexplored with respect to the vibrational coordinate.⁴¹ We also calculated the total vibrational frequency spectrum, which is the autocorrelation velocity spectrum of molecular dynamics (MD) trajectories from ref. 57, as shown in Fig. S17 (SI).⁵⁷ The calculated 125 cm⁻¹ mode matches closely with the assistive mode observed experimentally (136 cm⁻¹ in Fig. 4a) and the 318 cm⁻¹ vibration with the experimentally observed reactive mode (296 cm⁻¹ in Fig. 3d; see Fig. S18–S21, SI, and Materials and methods). The vectorial representation of the calculated reactive and assistive modes (Fig. S18–S21, SI) strongly corroborates the change in electron density upon the electronic transition (see Fig. 1d), reflecting the idea of nuclei responding to the resulting gradient or forces in the excited state potential.

Fig. 4 (a) Oscillation amplitude of the reorganization/assistive modes decaying on a sub-picosecond timescale (decay time τ in the inset). The two panels in (b) show the proposed non-linear mixing of the assistive modes to coherently generate the observed 80 cm⁻¹ non-impulsive vibration (black line). The modes implicated are unique combinations within the observed spectral modulation that could coherently drive a polarization at 80 cm⁻¹ as a part of overall relaxation to the closed-ring photoproduct (see SI Section 3). (c) Wavelet analysis of the Fourier filtered mode around 80 cm⁻¹ from experimental TA residual signal at 520 nm (see SI Section 1.3 and Fig. S7).

Upon closer inspection, we find an intense peak of 80 cm⁻¹ in the FFT power spectrum at 520 nm TA residual (see Fig. 3b). On Fourier filtering the TA residual (520 nm; see Fig. 3b) in the frequency domain around 80 cm⁻¹, we find unexpected temporal dynamics (Fig. S7, SI). The 80 cm⁻¹ mode appears to increase its amplitude until ∼1.7 ps, followed by a decrease in its amplitude with time, as depicted in Fig. 4b and c (this behavior is also seen in the wavelet analysis of the unfiltered TA residual of 520 nm; see Fig. S22, SI). This behavior is unlike vibrational wavepacket dynamics that are excited by the actinic pump through the linear superposition of states whose temporal evolution starts immediately after photoexcitation. This delayed response of the 80 cm⁻¹ vibration clearly shows that this mode is not excited by the actinic pump laser. Moreover, when we looked at the temporal evolution of the 80 cm⁻¹ mode across the probe wavelength range, we found that this effect is only prominent in the product absorption band, and the phase of the oscillation is π shifted from the ones that are not associated with the photoproduct absorption band (Fig. S23, SI). This unusual, non-impulsively driven, 80 cm⁻¹ mode highlights that the dynamics of this mode is strongly correlated to the dynamics of the photoproduct state. The calculated Raman spectrum of the isolated closed-ring fulgide molecule does not show Raman activity in the <150 cm⁻¹ frequency region (Fig. S24, SI), highlighting that the effect predominantly arises because of the crystalline environment. The previous work on the solution-state ultrafast spectroscopy of fulgide has not alluded to this effect.^40,42

Interestingly, we found that the decay of the assistive modes corresponds with the rise of the 80 cm⁻¹ mode (Fig. 4b and c). Furthermore, we observe the existence of pairs of assistive modes whose difference frequency is approximately 80 cm⁻¹, as illustrated in Fig. 4b. This suggests a possible mechanism of energy transfer involving the non-linear difference frequency mixing of assistive modes to coherently generate the 80 cm⁻¹ mode. The observation that the 80 cm⁻¹ mode does not relax on the sub-picosecond time scale, typically associated with intramolecular vibrational relaxation (IVR), indicates stronger coupling between the assistive modes compared to other competing channels for energy redistribution. This surprising observation also sheds light on the selectivity with which the assistive modes can mix non-linearly, highlighting the coordinate dependence of the anharmonicity in the photoproduct PES (see SI Section 3) and the far-from-equilibrium geometry after the ring-closing event. Previous reports on the non-impulsive generation of molecular wavepackets, which was seen as a growth in the amplitude of an oscillation, have been explained in the context of coherent formation of the bimolecular intermediate species in triiodide single crystal,⁶ sequential activation of the molecular breathing during spin-crossover reaction in Fe(II) complex⁹ and spin-vibronic mechanism in Pt–Pt dimer complexes in solution.^58,59 Chergui and co-workers have interpreted the growth of an oscillation amplitude in the Pt(pop) complex as the transfer of vibrational coherence to the triplet state, where the ISC (intersystem crossing) time is faster than the vibrational dephasing time and the transfer is mediated by an intermediate triplet state.⁶⁰ Chen and co-workers have shown that the non-impulsive growth of the Pt–Pt vibration is correlated to the proximity of the crossing point to the FC region.⁵⁸ In the case of fulgide, the crossing point is far away from the FC region, as shown in the theoretical PES calculation of the ring-closing reaction (Fig. 2d). This makes the idea of vibrational coherence transfer to the product state elusive in the fulgide single crystal. However, the possibility of non-linear mixing between assistive modes highlights the anharmonicity of the PES at the crossing point and opens new opportunities for a coherent vibrational energy transfer mechanism. The anharmonic coupling in this region is much stronger than in typical nonradiative relaxation processes, which is strongly enhanced by the large changes (1.8 Å) in the labile C–C positions (Cx–Cy distance highlighted in Fig. 1d) during the rapid passage through the CI. The resulting motions and changes in PES curvature would provide much larger anharmonicity and driving force for the non-linear mixing between the assistive modes relative to other points in the PES. It is the reaction forces involved in the ring-closing reaction that have revealed the anharmonic coupling between the dissipative assistive modes coupled to the reorganization coordinate.

Polarization dependence on observation of the 80 cm⁻¹ mode

The polarization dependence of the probe beam on the TA spectrum at 2 ps is plotted in Fig. S25a (SI). The FFT spectrum of the residual at 575 nm for the 20 degree polarization data only shows the presence of an intense 15 cm⁻¹ mode (Fig. S25d and e, SI). This indicates that the 80 cm⁻¹ mode has a strong orientation sensitivity in the crystal and is only observable along a particular probing direction of the crystal. This is further substantiated by the absence of the sub-picosecond component in the 20 degree polarization TA data (Fig. S25c, SI). The concomitant absence of the 80 cm⁻¹ mode, along with the sub-picosecond reorganization component in the global analysis of the 20 degree polarization TA data (Fig S25b–d, SI), strongly suggests that the reorganization of the molecular system is highly anisotropic, and the ring-closing reaction generates the growth of the 80 cm⁻¹ mode. The observation of the 15 cm⁻¹ mode (Fig. S25d, SI) highlights how the lattice or the molecular framework responds to the ring-closing reaction because of the induced strain in the crystal due to volume expansion (1310 to 1316 ų from open-ring to closed-ring structure).⁴⁴

Discussion

From our ultrafast spectroscopic measurements, we observe a unique ring-closing pathway, which distinguishes our results from the previous solution-state work.^38–40,42 Solid-state chemical processes provide unique insight into the degree of vibrational coupling and the influence of the surrounding lattice or bath on the reaction pathway. The molecular crystal, with its discrete nature and constrained motions, gives a well-defined geometry in comparison to the solution-state studies. In the solution-state reaction dynamics, the solvent defines the surrounding bath, which acts as a dissipative sink to accommodate relaxation to a global minimum. There are many different open conformations in the solution-state (existence of multiple isomers – α, β), which leads to parallel and competing E–Z isomerization pathways.^40,42 The ability of the molecule to undergo conformational changes with solvent fluctuations, leading to cavity fluctuations, allows the system to undergo different isomerization or conformations as observed in the solution-state. In the case of a single crystal, the lattice defines the system–bath interactions and the specific spatial relationships that are conserved along the reaction pathway. Specifically, the crystal packing and the geometry of the molecules in the crystal lattice inhibit the E–Z isomerization pathway.^24,43–45 It gives a well-defined starting geometry for the ring-closing reaction, which ultimately leads to the conservation of the phase correlation of the motions during the ultrafast ring-closing reaction. In the solution-state, the random orientation of the solvent and its dynamics would enable sampling other ring-closing/conformation pathways, destroying the coherence or phase correlations of motion. This point is evident as the same effects are not observed in the solution-state studies.^38–40,42 The coherent vibrational motions report directly on the anharmonicity of their respective potential, which is evolving along the reaction coordinate. The large driving force for the ring-closing reaction (1.5 eV, see Fig. 2d) associated with a 3.5 to 2.19 Å Cx–Cy bond distance change in 70 fs puts the system in a highly non-equilibrium and anharmonic geometry after the ultrafast curve-crossing event through the CI. These intricate details of the anharmonicity of the product PES surface are lost due to the random solvent and conformational motions of the fulgide system in solution. The random solvent cavity fluctuations within the solvation shell occur on the same timescale (sub-picosecond) as that of the structural reorganization (IVR) timescale in the ring-closing reaction and thereby wash out the details on the anharmonicity in the reaction surface.^61–64 We observe a well-defined single reaction coordinate due to the rigidity of the surrounding solid-state lattice as opposed to the evolution along multiple conformational isomerization pathways with different degrees of anharmonicity in solution-state processes that interfere with respect to the vibrational coherences.

Another issue that pertains to solid-state chemistry is the extent of excitonic delocalization upon photoexcitation. The absorption spectrum of the single crystal fulgide is broadened relative to the isolated gaseous state (Fig. 1a and S1, ESI) and the solution-state spectrum.^38–40,42 The decoherence dominates the absorption lineshape. Regardless, the initial absorption would be delocalized over several molecular sites. The excited state undergoes extremely fast decoherence due to both random fluctuations of the surrounding sites and sampling the reaction pathway, which localizes the excited state.^65,66 The subsequent sampling of the reactive pathway would completely destroy coherent coupling between sites, as the energy dissipation (>1000 cm⁻¹, 1.5 eV) would well exceed the excitonic coupling (order of 100 cm⁻¹ for typical lattice site distances between chromophores).²² Further, the fluctuations and energy dissipation processes are much larger than the thermal rms motions (300 cm⁻¹ at room temperature) of the intermolecular bath, which already leads to sub-100 fs decoherence timescales. The initial fluctuations at any one site along the reaction pathway would exceed the excitonic coupling and collapse the wavefunction to this site due to both the associated electronic decoherence and the initial energy dissipation involved in the atomic motions along the reaction coordinate. These processes would break the required condition of resonance for spatially extended states. Due to the localization of the electronic states, the phase correlation of a spatially extended wavefunction cannot be maintained over multiple sites. These bath fluctuations and more significantly, motion along the reaction coordinate would kill spatially extended coherences or phase relationships between multiple sites. Thus, the moment a spatially extended electronic wavefunction has a nuclear displacement along a reaction coordinate due to a change in electronic distribution, the local energy would completely decohere the spatially extended wavefunction, resulting in the collapse of the wavefunction to one molecular site. The same effect would break up vibronic couplings at the intramolecular level, which require resonant conditions.⁶⁷ In contrast, the spatially extended picture of vibrational modes at the intermolecular level, described as lattice phonons, is not broken by decoherence of the electronic excitation. These are spatially extended collective modes by definition and are not affected by the highly localized forces at the intramolecular level. These intermolecular modes are nevertheless involved in energy transfer via intermolecular coupling during vibrational relaxation and dissipation of excess energy along the reaction or reorganization coordinate. This energy dissipation process is the same physics that is involved in the vibrational cooling of excited molecules observed in molecular crystals.^68–70 It is possible to directly observe this coupling of reaction dynamics or dissipation to lattice phonon modes by looking at low q (reciprocal space vector) scattering between Bragg peaks in ultrafast electron diffraction measurements,^71–73 which is planned for future experiments.

Conclusion

In summary, we investigated the reaction dynamics for the fulgide photochromic system and identified the reaction kernel for the ring-closing reaction in a single crystal. The reaction kernel comprises two key reactive modes, 296 cm⁻¹ and 234 cm⁻¹. The physical picture is summarized in Fig. 5.

Fig. 5 Summary of the ultrafast ring-closing reaction in the fulgide single crystal.

The enormous reduction of dimensionality at the reaction saddle point, as depicted in Fig. 5, directs atomic motions along specific pathways, leading to a reproducible, effectively nonstochastic reaction mechanism that makes chemistry a transferable concept. In this regard, we were able to identify the tuning mode comprising a low-frequency 370 cm⁻¹ vibration responsible for bringing the system from the FC region to the intermediate state, and the other key reaction modes (296 cm⁻¹ and 234 cm⁻¹) that propagate the system from the intermediate state to the product surface through the CI. A similar intermediate state was observed in the ultrafast electron diffraction study of diarylethene with a shortened Cx–Cy distance of the reactive carbon atoms in a sub-picosecond timescale (3.71 Å at the FC region to 2.19 Å in the intermediate state).¹² In fulgide, the transformation to the intermediate state occurs in 70 fs, highlighting the relatively large forces acting on the reactive carbon atoms and the adjacent chemical moieties. This rapid relaxation involves significant atomic motions, which lead to a far-from-equilibrium position during the reactive crossing, unveiling the highly anharmonic part of the photoproduct PES. Thus, the surprising observation of a non-impulsive, coherent, vibrational energy transfer to an 80 cm⁻¹ mode profoundly highlights the anharmonicity at the reactive seam in the PES. This focuses the reaction along a single allowed pathway, constrained by larger barriers to lattice rearrangement in the solid-state relative to the solution-state, which enables the observation of changes in vibrational coherences that report on the anharmonicity with the unusual observation of the coherent mode mixing to generate the difference frequency (80 cm⁻¹) along the reorganization coordinate. The observation of this frequency during the ring-closing reaction in the crystal lattice effectively represents a spatially extended reorganization of the molecular geometry within the intermolecular lattice potential. This observation is strongly correlated to the reactive geometry in the crystalline environment and benefited from the relatively large atomic displacement of the reaction modes involved in the ring-closing reaction, unlike the solution-state measurements in fulgides, where the broad distribution of solvent configurations would wash out this feature to the reaction mechanism.^40,42 This discovery could open new ways of perceiving chemical reactions, with yet another mechanism collapsing the system onto a few highly non-linearly coupled coordinates over different length scales. The unique features of discrete reactive geometries within lattice confinement may inspire new synthetic strategies to allow scaling to mesoscale structures beyond what is possible in solution-state processes.

Materials and methods

Sample preparation

The fulgide was prepared following a published procedure (Cai et al.).²⁴ Characterization of fulgide by NMR (Nuclear Magnetic Resonance) spectroscopy matches the previously reported data (Cai et al.).²⁴

Optical spectroscopy

The steady-state absorption spectrum of fulgide was obtained using a homemade absorption set-up which included a fibre-coupled deuterium halogen lamp (DH-2000-BAL) and an Ocean Optics device as the detector. The transient absorption spectroscopy was performed in a home-built transient absorption (TA) spectrometer,²³ explained in ref. 23. The spectroscopy set-up used a Ti:sapphire laser (Coherent Legend Elite USP) with 40 fs, 800 nm output pulses running at 1 kHz. The 340 nm excitation pulses used in this work were generated using the sum frequency mixing of a Noncollinear Optical Parametric Amplifier (NOPA) signal around 570–680 nm and 800 nm fundamental pulses. The NOPA signal was generated using the doubled fundamental at 400 nm as the pump and the broadband white light as the seed. The 1 mJ, 800 nm output from the Ti:sapphire laser was initially split in a 70 : 30 ratio using a beam splitter. 30% of the split 800 nm beam was later used for sum frequency generation (SFG), and the remaining 700 μJ was used for NOPA signal generation. The 700 μJ arm was again split into a 90 : 10 ratio using a beam splitter where the 90% was used for second harmonic generation using a BBO crystal, barium beta borate (θ = 29.3°), and the remaining 10% 800 nm beam was used to generate white light seed by focussing the 800 nm to a 3 mm thick sapphire disk (half-inch diameter). The beam focusing was achieved with a spherical curved mirror. Later, by varying the distance between the focusing and collimating mirrors, the seed beam was refocused on the NOPA crystal (θ = 31.4°). The pump and the seed overlapped at the NOPA crystal at an external angle of ∼6°. The pump power at the NOPA crystal was of the order of 80–90 μJ. This gave a NOPA signal around the wavelength range of 570–680 nm, depending on the delay between the pump and the seed beam, which was later collimated using a spherical mirror. The NOPA pulse energies of 12–13 μJ were obtained, which were further compressed using an SF10 prism pair with a Brewster angle cut. After prism compression, the NOPA beam was directed to the 0.3 mm-thick SFG BBO crystal (θ = 31.5°) along with the 800 nm beam, which was initially split in the 70 : 30 ratio before the NOPA setup (see Fig. S26, SI). The 800 nm beam was path compensated with respect to the NOPA output beam to ensure temporal overlap (along with the spatial overlap) at the SFG crystal. The sum frequency signal is strongly dependent on the prism compression settings of the NOPA beam. The SFG signal was collimated using a spherical mirror and sent to a UV prism compression setup (UV fused silica prism) before transient absorption measurements. The SFG signal was around 340 nm with a 12–15 nm bandwidth (see Fig. S27, SI) and a pulse energy of 1 μJ (before prism compressor). This UV beam was used as the actinic pump for the fulgide single-crystal TA experiment. The pump beam was focused to a 95 μm beam size (FWHM, full-width half maximum) at the sample plane. The probe pulses were generated by focusing the 800 nm fundamental laser beam into a 2 mm-path-length quartz cuvette filled with ultrapure deionized water. The supercontinuum spanned from 400–750 nm spectral range. The probe beam at the sample plane had a focal diameter (FWHM) of 80 μm. The transmitted probe was routed to a home-built Czerny–Turner spectrograph coupled to a linear array CCD (Hamamatsu, S11155-2048-01 with C11165-01 driver).

For the transient absorption experiment, the fulgide single crystal was taped (Kapton) to a one inch quartz window, and the thickness of the crystal was sufficiently thin to allow the transmission of the broadband white-light probe. The fulgide crystals were cut along the ac plane perpendicular to the b axis and placed on the quartz window (Fig. 1b). Photoexcitation of the fulgide molecules with the UV pulses triggers an electrocyclic ring-closing reaction, forming a new closed-ring photoproduct (Fig. 1d). To ensure that the same initial conditions are maintained for each pump excitation event, we used a continuous wave (CW) laser (λ = 528 nm) to convert all the molecules back to their original form and reduced the repetition rate of the actinic pump laser to 125 Hz so that the interaction with the CW laser was long enough, in this case 8 ms, to convert all the closed-ring molecules at the focal spot to the starting open-ring form (Fig. S5, SI). It is important to mention that the closed-ring molecules absorb the visible part of the spectrum so that they can be reversibly transformed through visible light excitation back to the original open-ring form (Fig. S1a, SI). The excitation fluence of 2.73 mJ cm⁻² was used during the measurements, and it fell within the linear regime of photoexcitation (Fig. S28, SI) to avoid any multiphoton effect.⁷⁴ The photoexcitation was done with a circularly polarized 340 nm pulse and probed with a linearly polarized broadband white light ranging from 400 nm to 750 nm. Polarization of the white light was set to maximize the photoproduct signal, which was at 80 degrees with respect to the horizontal polarization of the fundamental 800 nm beam.

Computation

The crystal structure of the open and closed-ring fulgide was obtained from ref. 44, Harada et al. The optimization of atomic positions of open and closed-ring fulgide coordinates was performed with Density Functional Theory⁷⁵ (DFT) calculation using B3LYP functional⁷⁶ and 6-31G(d,p)⁷⁶ as the basis set. The absorption and the excited-state Raman spectrum (Time-Dependent DFT, TDDFT)^77,78 calculations were performed using the B3LYP functional and Def2SVP as the basis set. The calculated absorption spectrum, along with the oscillator strength, is shown in Fig. S1 (SI), and the excited-state Raman spectrum is shown in Fig. S16 (SI). All the calculations were done in Gaussian 016, Revision C.01.⁷⁹

When calculating the reaction pathway of the fulgide ring-closing reaction, we first optimize the minimum energy conical intersection (MECI) geometry of the fulgide between S₀ and S₁ states. Secondly, the NEB method⁸⁰ with frozen endpoints is used to search the reaction pathway from reactant to MECI, and from MECI to product. All the calculations are based on the complete active space self-consistent field (CASSCF) theory using the Terachem package,⁸¹ with 2 active electrons in 2 active orbitals. The 6-31G(d,p) basis set is used.

The total vibrational spectrum of molecular dynamics (MD) trajectories (trajectories obtained from ref. 57) is shown in Fig. S17 (SI). We can find the peaks around 125 cm⁻¹ and 318 cm⁻¹. To analyze the principal vibrational components contributing to each peak of the vibrational frequency spectrum, we performed Principal Component Analyses (PCA) of the MD trajectories, and the four principal component vectors are shown in Fig. S18–S21 (SI). To study the range of vibrational frequency distributions for each PCA component, we calculated the autocorrelation velocity spectrum of each PCA component from the projections of the velocity of MD trajectories to each PCA component. The result is also shown in Fig. S18–S21 (SI). The peak of the vibrational frequency spectrum around 130 cm⁻¹ is contributed by the first PCA component, and the third PCA component exhibits a peak of around 150 cm⁻¹. The 296 cm⁻¹ peak is contributed from the first, second, and fourth PCA components around 300 cm⁻¹. The real vibrational vector is a superposition of these PCA components.

Author contributions

RJDM formulated and supervised the project. The planning and discussion of the experiments were done by SM, SFB, and RJDM. The transient absorption experiment, TDDFT calculations, and data analysis were performed by SM. MZ, XD, and ZL performed the NEB calculation of the PES and the MD simulation trajectory analysis. JC and SM (Scott Murphy) provided the single crystal fulgide molecule. The manuscript was written with contributions from all authors.

Conflicts of interest

There are no conflicts of interest.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5sc03764f.

Acknowledgements

This work was supported by the Natural Sciences and Engineering Research Council of Canada (RJDM, RSM). Support for the early phase of the work was also provided by the Max Planck Society (RJDM). Support for ZL was provided by the National Natural Science Foundation of China (Grants No. 12174009, 12234002, 92250303), the National Key R&D Program of China (Grant No. 2023YFA1406801), and the Beijing Natural Science Foundation (Grant No. Z220008). The authors thank Compute Canada (Niagara and Cedar) for the calculations presented in this manuscript.

References

1. A. A. Ischenko, P. M. Weber and R. J. D. Miller, Chem. Rev., 2017, 117, 11066–11124.
2. Z. Li, S. Gyawali, A. A. Ischenko, S. Hayes and R. J. D. Miller, ACS Photonics, 2020, 7, 296–320.
3. R. J. D. Miller, Annu. Rev. Phys. Chem., 2014, 65, 583.
4. P. J. M. Johnson, A. Halpin, T. Morizumi, V. I. Prokhorenko, O. P. Ernst and R. J. D. Miller, Nat. Chem., 2015, 7, 980–986.
5. J. C. Owrutsky, D. Raftery and R. M. Hochstrasser, Annu. Rev. Phys. Chem., 1994, 45, 519–555.
6. R. Xian, G. Corthey, D. M. Rogers, C. A. Morrison, V. I. Prokhorenko, S. A. Hayes and R. J. D. Miller, Nat. Chem., 2017, 9, 516–522.
7. T. Ishikawa, S. A. Hayes, S. Keskin, G. Corthey, M. Hada, K. Pichugin, A. Marx, J. Hirscht, K. Shionuma, K. Onda, Y. Okimoto, S. Koshihara, T. Yamamoto, H. Cui, M. Nomura, Y. Oshima, M. Abdel-Jawad, R. Kato and R. J. D. Miller, Science, 2015, 350, 1501–1505.
8. Y. Jiang, L. C. Liu, A. Sarracini, K. M. Krawczyk, J. S. Wentzell, C. Lu, R. L. Field, S. F. Matar, W. Gawelda, H. M. Müller-Werkmeister and R. J. D. Miller, Nat. Commun., 2020, 11, 1530.
9. M. Cammarata, R. Bertoni, M. Lorenc, H. Cailleau, S. Di Matteo, C. Mauriac, S. F. Matar, H. Lemke, M. Chollet, S. Ravy, C. Laulhé, J.-F. Létard and E. Collet, Phys. Rev. Lett., 2014, 113, 227402.
10. M. Gao, C. Lu, H. Jean-Ruel, L. C. Liu, A. Marx, K. Onda, S. Y. Koshihara, Y. Nakano, X. Shao, T. Hiramatsu, G. Saito, H. Yamochi, R. R. Cooney, G. Moriena, G. Sciaini and R. J. D. Miller, Nature, 2013, 496, 343–346.
11. H. Jean-Ruel, R. R. Cooney, M. Gao, C. Lu, M. A. Kochman, C. A. Morrison and R. J. D. Miller, J. Phys. Chem. A, 2011, 115, 13158–13168.
12. H. Jean-Ruel, M. Gao, M. A. Kochman, C. Lu, L. C. Liu, R. R. Cooney, C. A. Morrison and R. J. D. Miller, J. Phys. Chem. B, 2013, 117, 15894–15902.
13. C. G. Elles and F. F. Crim, Annu. Rev. Phys. Chem., 2006, 57, 273–302.
14. B. K. Carpenter, J. N. Harvey and A. J. Orr-Ewing, J. Am. Chem. Soc., 2016, 138, 4695–4705.
15. V. A. Apkarian and N. Schwentner, Chem. Rev., 1999, 99, 1481–1514.
16. R. Field, L. C. Liu, W. Gawelda, C. Lu and R. J. D. Miller, Chem.–Eur. J., 2016, 22, 5118–5122.
17. A. J. Orr-Ewing, Chem. Soc. Rev., 2017, 46, 7597–7614.
18. J. Harada and K. Ogawa, Chem. Soc. Rev., 2009, 38, 2244.
19. R. Crespo-Otero, Q. Li and L. Blancafort, Chem.–Asian J., 2019, 14, 700–714.
20. R. P. Xian, S. A. Hayes, G. Corthey, C. A. Morrison, A. Marx, H. Daoud, C. Lu and R. J. D. Miller, J. Am. Chem. Soc., 2025, 147, 28973–28980.
21. V. Tiwari, M. Gallagher-Jones, H. Hwang, H.-G. Duan, A. I. Kirkland, R. J. D. Miller and A. Jha, ACS Phys. Chem. Au, 2024, 4, 660–668.
22. H. Hwang, A. Mackenzie, M. A. Kochman, M. Gallagher-Jones, V. Tiwari, S. F. Bittmann, F. Tellkamp, S. Botchway, A. I. Kirkland, H.-G. Duan, R. J. D. Miller and A. Jha, J. Phys. Chem. B, 2025, 129, 7982–7994.
23. K. M. Siddiqui, G. Corthey, S. A. Hayes, A. Rossos, D. S. Badali, R. Xian, R. S. Murphy, B. J. Whitaker and R. J. D. Miller, CrystEngComm, 2016, 18, 7212.
24. J. Cai, A. Farhat, P. B. Tsitovitch, V. Bodani, R. D. Toogood and R. S. Murphy, J. Photochem. Photobiol., A, 2010, 212, 176–182.
25. E. Uhlmann and G. Gauglitz, J. Photochem. Photobiol., A, 1996, 98, 45–49.
26. M. Lrie, in Molecular Switches, 2001, pp. 37–62.
27. Y. Yokoyama, in Molecular Switches, 2001, pp. 107–121.
28. V. Balzani, A. Credi, B. Ferrer, S. Silvi and M. Venturi, Top. Curr. Chem., 2005, 262, 1–27.
29. T. lkeda and A. Kanazawa, in Molecular Switches, 2001, pp. 363–397.
30. A. P. de Silva, N. D. McClenaghan and C. P. McCoy, in Molecular Switches, 2001, pp. 339–361.
31. S. Shinkai, in Molecular Switches, 2001, pp. 281–307.
32. I. Willner and B. Willner, in Molecular Switches, 2001, pp. 165–218.
33. T. Mrozek, A. Ajayaghosh and J. Daub, in Molecular Switches, 2001, pp. 63–106.
34. M. Irie, Chem. Rev., 2000, 100, 1685–1716.
35. M. Irie, T. Fukaminato, K. Matsuda and S. Kobatake, Chem. Rev., 2014, 114, 12174.
36. N. Tamai and H. Miyasaka, Chem. Rev., 2000, 100, 1875–1890.
37. M. Hada, Y. Nishina and T. Kato, Acc. Chem. Res., 2021, 54, 731–743.
38. F. Renth, M. Foca, A. Petter and F. Temps, Chem. Phys. Lett., 2006, 428, 62–67.
39. F. Renth, R. Siewertsen, F. Strübe, J. Mattay and F. Temps, Phys. Chem. Chem. Phys., 2014, 16, 19556–19563.
40. R. Siewertsen, F. Renth, F. Temps and F. Sönnichsen, Phys. Chem. Chem. Phys., 2009, 11, 5952.
41. B. Heinz, S. Malkmus, S. Laimgruber, S. Dietrich, C. Schulz, K. Rück-Braun, M. Braun, W. Zinth and P. Gilch, J. Am. Chem. Soc., 2007, 129, 8577–8584.
42. R. Siewertsen, F. Strübe, J. Mattay, F. Renth and F. Temps, Phys. Chem. Chem. Phys., 2011, 13, 3800.
43. T. Tayu and S. Kurita, J. Phys. Chem. Solids, 1996, 57, 475–482.
44. J. Harada, R. Nakajima and K. Ogawa, J. Am. Chem. Soc., 2008, 130, 7085–7091.
45. S. Kobatake and M. Irie, Chem. Lett., 2004, 33, 904–905.
46. H. Port, P. Gärtner, M. Hennrich, I. Ramsteiner and T. Schöck, Mol. Cryst. Liq. Cryst., 2005, 430, 15–21.
47. M. S. Schuurman and A. Stolow, Annu. Rev. Phys. Chem., 2018, 69, 427–450.
48. M. Liebel, C. Schnedermann, T. Wende and P. Kukura, J. Phys. Chem. A, 2015, 119, 9506–9517.
49. G. D. Scholes, G. R. Fleming, L. X. Chen, A. Aspuru-Guzik, A. Buchleitner, D. F. Coker, G. S. Engel, R. van Grondelle, A. Ishizaki, D. M. Jonas, J. S. Lundeen, J. K. McCusker, S. Mukamel, J. P. Ogilvie, A. Olaya-Castro, M. A. Ratner, F. C. Spano, K. B. Whaley and X. Zhu, Nature, 2017, 543, 647–656.
50. S. Takeuchi, S. Ruhman, T. Tsuneda, M. Chiba, T. Taketsugu and T. Tahara, Science, 2008, 322, 1073–1077.
51. K. Ishii, S. Takeuchi and T. Tahara, J. Phys. Chem. A, 2008, 112, 2219–2227.
52. J. Chesnoy and A. Mokhtari, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3566–3576.
53. S. Mitra, S. R. K. Ainavarapu and J. Dasgupta, J. Phys. Chem. B, 2022, 126, 5390–5399.
54. T. J. Martinez, Nature, 2010, 467, 412–413.
55. C. Schnedermann, M. Liebel and P. Kukura, J. Am. Chem. Soc., 2015, 137, 2886–2891.
56. A. Kühl and W. Domcke, J. Chem. Phys., 2002, 116, 263–274.
57. M. A. Kochman, T. Gryber, B. Durbeej and A. Kubas, Phys. Chem. Chem. Phys., 2022, 24, 18103–18118.
58. S. R. Rather, N. P. Weingartz, S. Kromer, F. N. Castellano and L. X. Chen, Nature, 2023, 620, 776–781.
59. P. Kim, S. Roy, A. J. S. Valentine, X. Liu, S. Kromer, T. W. Kim, X. Li, F. N. Castellano and L. X. Chen, Chem. Sci., 2024, 15, 14766–14777.
60. R. Monni, G. Capano, G. Auböck, H. B. Gray, A. Vlček, I. Tavernelli and M. Chergui, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 6396–6403.
61. V. F. Crum, L. M. Kiefer and K. J. Kubarych, J. Chem. Phys., 2021, 155, 134502.
62. J.-Y. Liu, W.-H. Fan, K.-L. Han, W.-Q. Deng, D.-L. Xu and N.-Q. Lou, J. Phys. Chem. A, 2003, 107, 10857–10861.
63. M. Berg and D. A. V. Bout, Acc. Chem. Res., 1997, 30, 65–71.
64. R. K. Venkatraman and A. J. Orr-Ewing, Acc. Chem. Res., 2021, 54, 4383–4394.
65. J. Cao, R. J. Cogdell, D. F. Coker, H.-G. Duan, J. Hauer, U. Kleinekathöfer, T. L. C. Jansen, T. Mančal, R. J. D. Miller, J. P. Ogilvie, V. I. Prokhorenko, T. Renger, H.-S. Tan, R. Tempelaar, M. Thorwart, E. Thyrhaug, S. Westenhoff and D. Zigmantas, Sci. Adv., 2020, 6, eaaz4888.
66. H.-G. Duan, M. Thorwart and R. J. D. Miller, J. Chem. Phys., 2019, 151, 114115.
67. V. Tiwari, W. K. Peters and D. M. Jonas, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 1203–1208.
68. J. R. Hill and D. D. Dlott, J. Chem. Phys., 1988, 89, 842–858.
69. J. R. Hill and D. D. Dlott, J. Chem. Phys., 1988, 89, 830–841.
70. D. D. Dlott, J. Lumin., 1990, 45, 397–400.
71. G. Sciaini, M. Harb, S. G. Kruglik, T. Payer, C. T. Hebeisen, F.-J. M. zu Heringdorf, M. Yamaguchi, M. H. Hoegen, R. Ernstorfer and R. J. D. Miller, Nature, 2009, 458, 56–59.
72. L. Kremeyer, T. L. Britt, B. J. Siwick and S. C. Huberman, Struct. Dyn., 2024, 11, 024101.
73. M. R. Otto, J.-H. Pöhls, L. P. René de Cotret, M. J. Stern, M. Sutton and B. J. Siwick, Sci. Adv., 2025, 7, eabf2810.
74. R. J. D. Miller, O. Paré-Labrosse, A. Sarracini and J. E. Besaw, Nat. Commun., 2020, 11, 1240.
75. K. A. Baseden and J. W. Tye, J. Chem. Educ., 2014, 91, 2116–2123.
76. J. Tirado-Rives and W. L. Jorgensen, J. Chem. Theory Comput., 2008, 4, 297–306.
77. M. A. L. Marques and E. K. U. Gross, Annu. Rev. Phys. Chem., 2004, 55, 427–455.
78. J. Schirmer, Phys. Chem. Chem. Phys., 2025, 27, 4992–5005.
79. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16, Revision C.01, Gaussian, Inc., Wallingford, 2016.
80. V. Ásgeirsson, B. O. Birgisson, R. Bjornsson, U. Becker, F. Neese, C. Riplinger and H. Jónsson, J. Chem. Theory Comput., 2021, 17, 4929–4945.
81. S. Seritan, C. Bannwarth, B. S. Fales, E. G. Hohenstein, C. M. Isborn, S. I. L. Kokkila-Schumacher, X. Li, F. Liu, N. Luehr, J. W. Snyder Jr, C. Song, A. V Titov, I. S. Ufimtsev, L.-P. Wang and T. J. Martínez, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2021, 11, e1494.

---

Footnote

† We dedicate this work to the memory of Jianxin Cai, who passed away on September 4, 2020.

---