Excited-state intramolecular proton transfer and competing pathways in 3-hydroxychromone: a non-adiabatic dynamics study

Abstract

Excited-state intramolecular proton transfer (ESIPT) is a fundamental photochemical process in which photoexcitation induces proton transfer within a molecule, leading to the formation of a tautomeric excited state. It was observed experimentally that the 3-hydroxychromone (3-HC) system exhibits two distinct proton-transfer time scales upon excitation to the lowest bright singlet excited state: an ultrafast component on the femtosecond time scale and a slower one on the picosecond time scale, largely insensitive to solvent effects. Up to now, the microscopic origin of the second time constant has only been hypothesised. Here, using mixed quantum-classical non-adiabatic dynamics simulations, we explicitly observe the two ESIPT time constants and we rationalise the origin of the second time scale by the presence of a competitive out-of-plane hydrogen torsional motion. Comprehensive analysis of the excited-state potential energy surfaces and non-adiabatic trajectories enables us to construct an explicit reaction network for 3-HC, delineating the interplay between direct ESIPT and torsion-mediated pathways. This unified mechanistic framework reconciles the coexistence of ultrafast and slower ESIPT components, offering new insights into the non-adiabatic excited-state dynamics of the system.

---

1. Introduction

Excited-state intramolecular proton transfer (ESIPT) refers to a photochemical event where photoexcitation triggers proton transfer within a molecule, producing a tautomeric excited state. The earliest description of this process was provided by Weller,^1,2 noting that the remarkable Stokes shift (of nearly 10 000 cm⁻¹) detected in the fluorescence of salicylic acid could be explained by a fast excited-state proton transfer. Indeed, the ESIPT induces a major structural reorganisation compatible with the large observed shift. Later, Ottersted³ and Kasha⁴ also focused on the interpretation of this phenomenon in other molecular systems.

Since then, significant experimental and theoretical computational efforts have been dedicated to understanding excitation-induced proton transfer in a wide range of molecular systems. This resulted in numerous applications of ESIPT systems as fluorescent probes,^5–7 photostabilizers,^8,9 organic light-emitting diodes (OLEDs),^10–15 and in recognition of its central role in several biochemical processes.^16–20 The time scale of the ESIPT process is considered rather short,²¹ on the femtosecond to picosecond time scale, depending on the molecular properties that can also affect the yield. In particular, the presence of an activation barrier along the coordinate that connects the relevant tautomers involved in the proton transfer and/or the competition with other molecular motions can lead to a notable reduction of the rate and yield in comparison to the case in which the ground-state tautomer is not stable in the excited-state and transforms, in a direct fashion, toward the photoproduct. In the case of incomplete tautomerization, the emission spectrum can exhibit two bands: one with a modest Stokes shift from the reactant state and another with a large Stokes shift from the ESIPT product. This aspect is also referred to as dual fluorescence.

Understanding the dynamics and mechanisms of such a process is important not only from a fundamental photochemical and photobiological perspective but also crucial for the rational design of molecules with tailored optical properties. The experimental investigation of the ultrafast ESIPT requires high-resolution time-resolved spectroscopic techniques. Among them, femtosecond stimulated Raman spectroscopy (FSRS) offers valuable insight into skeletal motions during proton transfer,^22,23 as well as time-resolved fluorescence, UV/vis and IR absorption spectroscopy.²⁴

Since it is difficult from experiments to obtain information about transition state structures, energy barriers, possible conical intersections between the involved excited states, or other critical features of the reaction pathway, such data, essential for unraveling the detailed mechanism of ESIPT, are often obtained from state-of-the-art ab initio excited-state calculations. In recent years, theoretical investigations have therefore been carried out on a wide variety of ESIPT systems using diverse computational strategies: among them, time dependent-density functional theory (TD-DFT),^25–27 post-Hartree–Fock approaches, such as second order approximated coupled cluster (CC2)^28–30 or algebraic diagrammatic construction approximated to the second order (ADC(2)),^31,32 and multireference methods such as complete active space self-consistent field (CASSCF), in some studies with the inclusion of the perturbative correction (CASPT2).^33–35 Along with static approaches, non-adiabatic dynamics simulations, employing multi-configuration time-dependent Hartree (MCTDH),^36–39ab initio multiple spawning (AIMS),^40,41 and trajectory surface hopping (TSH),^42–45 were performed to explicitly follow the proton transfer on the excited state surface(s).

In recent years, 3-hydroxychromone (3-HC), depicted in Fig. 1a, and its derivatives have drawn significant attention for their photochemical properties and potential applications as molecular probes that exploit their spectrally distinct dual fluorescence.²⁴ Indeed, upon excitation to the first (bright, ππ*) singlet excited-state, the 3-HC molecule can easily access the S₁ minimum and fluoresce or undergo fast ESIPT reaction. Chevalier et al.²⁴ reported a distinctive behaviour of the 3-HC molecule, namely the presence of two ESIPT rate constants: a fast component on the femtosecond time scale and a slower one on the picosecond time scale, independent of the solvent environment. In polar protic solvents, such as water, the formation of the conjugated base of the molecule (i.e., 3-HC deprotonated at the hydroxyl functional group) and the formation of solute–solvent intermolecular hydrogen bonds could be competitive with the ESIPT process, and was thus hypothesised to give rise to the slower component in the picosecond time scale. Furthermore, in aprotic apolar solvents such as methylcyclohexane, a slower ESIPT component was also observed, and attributed either to (i) intramolecular vibrational relaxation (IVR), which reduces the excitation of proton transfer-promoting modes, and/or (ii) transient access to the excited trans-enol conformation (see Fig. 1b).

Fig. 1 (a) Enol and keto tautomeric forms of the 3-hydroxychromone (3-HC). (b) (cis-)enol to trans-enol conformational equilibrium.

Among theoretical-computational studies on the 3-HC system^43,45–48 only Perveaux et al.³⁷ specifically addressed the origin of the slower (picosecond time scale) ESIPT rate constant. They localised a previously undocumented conical intersection between the first two singlet excited states (the bright, ππ* and dark, nπ*, at the Franck–Condon (FC) point) in 3-HC. This intersection exhibits a planar geometry (C_s symmetry) and is readily accessible from the FC point. Using multilayer MCTDH quantum dynamics simulations, the authors proposed that following excitation to S₁, the 3-HC nuclear wave packet can reach this S₁/S₂ conical intersection, enabling ultrafast population transfer to the nπ* state along the non-reactive direction of the ESIPT coordinate. These findings were later qualitatively supported by the quantum dynamics simulations of Anand et al.³⁸ However, the theoretical reproduction of two ESIPT time constants from non-adiabatic dynamics simulations and an explicit correlation with the hypothesised mechanism is still lacking. We present here a theoretical-computational investigation of the ESIPT reaction in the 3-HC molecule (Fig. 1a) by means of on-the-fly mixed quantum-classical nonadiabatic dynamics simulations. Our aim is to contribute to a thorough mechanistic explanation of the factors that lead to the observed double time constants associated with the ESIPT reaction. A deeper understanding of the mechanisms competing with the proton transfer reaction could aid in the design of more efficient 3-HC-based molecular probes.

The article is structured as follows: in Section 2 the computational details used in this work are presented. In Section 3, the results are exposed and discussed. Section 4 offers some concluding remarks.

2. Theoretical methods and computational details

2.1. Electronic structure methods

At the FC point, the 3-HC molecule exhibits two low-lying excited states that are nearly degenerate. Small geometric displacements from the FC point can readily alter the state character owing to this near-degeneracy of S₁ and S₂ states. This emphasizes the necessity of (i) a rigorous characterization of the energetics and state characters involved, and (ii) the explicit inclusion of non-adiabatic effects in describing nuclear dynamics, even in the vicinity of the FC point. Electronic structure calculations by Perveaux et al.³⁷ examined the reliability of TD-PBE0 in capturing the S₁ (bright, ππ*) and S₂ (dark, nπ*) states at the FC point and critical regions of the potential energy surfaces, through comparisons with the higher-level CC2, ADC(2), and EOM-CCSD methods. Anand et al.³⁸ assessed the TD-DFT approach, reporting an S₁–S₂ energy gap at the FC point between 0.01 and 0.10 eV, depending on the functional, in qualitative agreement with the 0.15 eV gap obtained from their EOM-CCSD calculations.

On this basis, TD-PBE0⁴⁹ was chosen here to compute the potential energy surfaces and their couplings for the surface-hopping simulations. This functional offers a favourable balance between accuracy and computational cost and has been shown to reliably describe the electronic structure at key points along the ESIPT coordinate.³⁷ The calculations employed the cc-pVDZ⁵⁰ basis set. To assess the influence of the basis set, the critical points along the ESIPT pathway were re-optimised using the cc-pVTZ and aug-cc-pVDZ bases and the energy profiles were re-computed. The S₁–S₂ energy gap at the FC point from our TD-PBE0/cc-pVDZ (cc-pVTZ/aug-cc-pVDZ) calculations is 0.12 (0.08/0.10) eV. The results obtained from the non-adiabatic dynamics simulations will be compared (see below) with experimental data obtained in methylcyclohexane. The influence of the solvent was also tested on the aforementioned ESIPT profile using the dielectric continuum description provided by the IEFPCM.⁵¹ All the calculations were carried out with Gaussian 16, rev. A03.⁵² The tight SCF convergence criteria and the fine integration grid were requested. Several key conical intersections were optimised with the Orca 6.1 package⁵³ using the updated branching plane (UBP) method, which circumvents the need for explicit non-adiabatic coupling matrix elements.⁵⁴ The optimisation calculations employed the same level of theory as above, namely TD-DFT with the PBE0 functional and the cc-pVDZ basis set.

2.2. Propagation of the coupled electron-nuclear dynamics

Trajectory surface-hopping simulations were carried out with the SHARC software package (version 3.0.1).⁵⁵ Initial conditions were generated from the 0 K Wigner distribution based on the ground-state vibrational frequencies of the (cis-)enol conformer. A total of 10 000 geometries were sampled, and for each, the excitation energies and oscillator strengths of the first two singlet states (S₁ and S₂) were computed at the TD-PBE0/cc-pVDZ level via the SHARC-GAUSSIAN interface. These data were employed to construct the absorption spectrum using the nuclear ensemble approach,⁵⁶ where each line was convoluted with a Gaussian function of 0.1 eV full width at half maximum. No energy shift was applied to the computed spectrum with respect to the experimental one.

From this ensemble, the initial conditions for the non-adiabatic dynamics simulations were selected based on excitation energy and oscillator strength. Specifically, an energy window of 0.2 eV centered at 3.9 eV (ca. 320 nm, the excitation wavelength used by Chevalier et al.²⁴ in their time-resolved UV-vis experiments) was applied. Within this range, the brightest geometries for each state were chosen following the stochastic selection procedure of Piteša et al.,⁵⁷ as implemented in SHARC. This yielded 311 and 86 trajectories initiated on the S₁ and S₂ states, respectively.

A nuclear time step of 0.5 fs was used. Three electronic states, S₀–S₂, were considered in the non-adiabatic dynamics simulations. The electronic wave function was propagated using the local diabatization method⁵⁸ with 25 substeps. Granucci and Persico decoherence scheme was used, with a factor of 0.1 a.u. of energy.⁵⁹ Hopping probabilities between S₁ and S₂ were computed from the time evolution of the electronic amplitudes. An energy threshold of 0.10 eV was used to enforce a hop between an active excited state and the ground state, due to the well-documented deficiency of TD-DFT in the description of conical intersections involving the electronic ground state.⁶⁰ After a hop, the nuclear velocity vector was rescaled isotropically to conserve the total energy. All trajectories exhibit total energy conservation within the selected 0.5 eV threshold in both S₁ and S₂ ensembles. No trajectories terminated prematurely due to convergence failures in the electronic structure calculations.

2.3. Analysis

To gain deeper insight into the ESIPT process and competing relaxation pathways, structural and statistical analyses on the full ensemble of trajectories were performed. A chemically intuitive and practical way to characterize the ESIPT reaction is through the distances between the transferring proton and the two oxygen atoms involved. We denote these bond distances as d(O_d–H) and d(H–O_a), corresponding to the donor oxygen–hydrogen and hydrogen–acceptor oxygen distances, respectively (see Fig. 1b). To probe the role of torsional motions in both tautomers, we additionally monitored the dihedral angles ϕ₁ and ϕ₂, defined by the atoms (H, O_d, C₁, C₂) and (H, O_a, C₂, C₁), respectively (see Fig. 1b). All structural analyses were performed with the aid of the MDTraj software package.⁶¹

Along the mixed quantum-classical trajectories, the molecule was classified as either enol or keto (see Fig. 1a) tautomer to obtain the time-dependent yield of the ESIPT reaction (see below). The time-dependent yield of the keto form was monitored by analysing the distances d(O_d–H) and d(H–O_a) along the trajectories. More precisely, at each time step, the structure was assigned to the keto form when d(H–O_a) < d(O_d–H). The (cis-)enol form is the most stable tautomer in the ground state and is the form from which the dynamics starts. In contrast, the keto tautomer is stabilized in the S₁ state and is formed along the ESIPT reaction coordinate. To further analyse the nuclear motions along the non-adiabatic trajectories, a principal component analysis (PCA)⁶² was performed on the full ensemble of structures. This approach allowed us to identify the dominant collective motions contributing to the ESIPT process and other relevant conformational changes.

3. Results and discussion

The results are now presented, starting with a description of the two lowest-lying valence excited electronic states at the ground state equilibrium geometry of 3-HC, its UV-vis absorption spectrum, and some critical points on the excited state surfaces along the ESIPT coordinate and possible competitive pathways. Then, the photodynamics of 3-HC, as described by surface hopping simulations, is discussed. In particular, we will focus on the analyses of the adiabatic electronic state populations and the key nuclear degrees of freedom involved in the ESIPT reaction, the kinetic modelling of the proton transfer, and the definition of a complete reaction network that connects all the relevant species involved in their different electronic states.

3.1. Vertical absorption spectrum

The vertical absorption spectrum of 3-HC in the gas phase was computed via the nuclear ensemble approach⁵⁶ from all (10 000) geometries sampled from the Wigner distribution. Vertical excitation energies at each geometry were evaluated at TD-PBE0/cc-pVDZ level of theory, and a Gaussian line shape with FWHM of 0.1 eV was applied. The resulting spectrum is shown in Fig. 2. The computed maximum absorption occurs at 308 nm, in good agreement with the experimental value of ca. 313 nm reported by Chevalier et al.²⁴ in methylcyclohexane. At the FC point of the ground state (cis-)enol tautomer, the first two excited states, S₁ and S₂, are close in energy, with S₁ exhibiting bright ππ* character and S₂ being dark with predominantly nπ* character. Due to the near-degeneracy and the vibronic coupling of these states, even slight distortions of the molecular geometry can lead to mixing or even re-ordering of the states, as also observed by Nag et al.⁴³ Owing to vibronic coupling between the first two excited states, the dark state can acquire intensity from the bright state when the molecular geometry is displaced along specific normal modes away from the FC region. Due to this intensity borrowing, both states are photo-excited and must be considered as initial states in the dynamics simulations. In addition, the absorption spectrum was decomposed into the contribution from individual diabatic states (i.e., ππ* and nπ* characters at the FC point) and reported in Fig. S1.

Fig. 2 TD-PBE0/cc-pVDZ vertical absorption spectrum of 3-HC obtained via the nuclear ensemble approach. The gray band represents the energy window used for the initial condition selection. No energy shift was applied to the computed spectrum. The digitalised experimental absorption spectrum (band C) from Chevalier et al.²⁴ is reported as a blue dashed line.

3.2. ESIPT coordinate and competitive pathway characterisations

The ESIPT coordinate was characterised by locating key points on the potential energy surfaces and interpolating between them to obtain the energy profile along a collective coordinate that connects them. The structures considered include the FC point, the S₁ minima corresponding to the (cis-)enol and keto tautomers (see Fig. 1a), and the S₁ transition state (TS) connecting them. Linear interpolations in internal coordinates (LIIC) were performed between these critical points to map the ESIPT coordinate. The potential energy curves of the first three electronic states along these LIIC paths, computed at the TD-PBE0/cc-pVDZ level, are shown in Fig. 3. To assess the effect of basis set size, the critical points were also optimised with the cc-pVTZ and aug-cc-pVDZ basis set, and the LIIC profiles were recomputed. Expanding the basis from cc-pVDZ to cc-pVTZ or aug-cc-pVDZ has only a minor impact on the potential energy curves along the ESIPT coordinate (see Fig. S2).

Fig. 3 Linear interpolation in internal coordinates (LIIC) along the critical points that characterise the ESIPT coordinate at TD-PBE0/cc-pVDZ level of theory. The energy barrier on the S₁ potential energy surface between the (cis-)enol minimum to the TS is indicated by the two arrows.

3-HC exhibits a small barrier of approximately 0.04 eV along the ESIPT coordinate on the S₁ potential energy surface, computed from the S₁ (cis-)enol minimum to the TS on the same potential energy surface. In contrast, a significantly larger barrier of about 0.97 eV is observed along the same coordinate on the S₂ state potential energy surface, in qualitative agreement with the 0.95 eV barrier reported by Nag et al.⁴³ at the TD-B3LYP/6-311++G** level. It is noted that optimising the ESIPT TS on the S₂ state gives a barrier of 0.46 eV, which is still approximately ten times higher than the one on the S₁ surface. In the ground-state, the keto form is not stable and all the optimisation attempts on this potential energy surface ended in the enol form. Consequently, the gradient at the FC point on the S₁ surface, together with the modest energy barrier along the proton-transfer coordinate, is expected to drive the system toward enol-to-keto tautomerization. In view of the substantially higher barrier on S₂, only a minor fraction of photoexcited molecules is anticipated to undergo ESIPT on this surface, as will be confirmed by the subsequent analysis of the non-adiabatic dynamics simulations (see below).

Previous studies^24,37 suggested that other relaxation processes may compete with proton transfer following photoexcitation of 3-HC, potentially explaining the double time constant observed for the ESIPT reaction. To investigate this, we first located a minimum-energy conical intersection (MECI) between the S₁ and S₂ states near the FC point, which was previously identified by Perveaux et al.³⁷ as a possible origin of the second ESIPT time constant. This MECI adopts a planar geometry (C_s symmetry, as the stationary points located and discussed so far) and is structurally very similar to the FC point, differing mainly by an increase of approximately 0.09 Å in the d(H–O_a) distance. Notably, we find the MECI to lie 0.10 eV below the FC point, in agreement with Anand et al.,³⁸ suggesting that it is readily accessible to the nuclear wave packet following photoexcitation. Perveaux et al.³⁷ reported that, possibly due to a different optimisation algorithm, the corresponding conical intersection lies slightly above the FC point, yet still confirms the accessibility of the S₁/S₂ seam in the vicinity of the FC point.

After crossing the MECI seam, the S₁ state, now nπ* in character, exhibits two symmetrically equivalent minima along the torsional coordinate defined by a ϕ₁ dihedral angle at approximately ±21.6 degrees. The loss of C_s symmetry can be due to the pseudo-Jahn–Teller effect, arising from the vibronic coupling of the two states in the region of (near-)degeneracy.⁶³ We refer to these stationary points as S₁ torsional minima (or in brevity, tor-enol). To characterise the energy profile along this pathway, we performed a LIIC from the MECI to one of the two symmetrically equivalent torsional minima, followed by a relaxed scan from this minimum to the planar trans-enol conformer (i.e., when ϕ₁ = ±180.0 degrees). The resulting S₁ energy curves are shown in Fig. 4. An exactly symmetric surface is expected when moving in the opposite (negative) direction along the ϕ₁ torsion due to the symmetry of the system. From the MECI, the system can relax along ϕ₁ to reach the torsional minimum without encountering any barrier. From this minimum, a modest barrier of 0.09 eV must be overcome to reach the trans-enol conformer, in line with previous calculations.³⁷ Given the energy available upon vertical photoexcitation, the system is expected to explore this torsional coordinate, as will be confirmed by the non-adiabatic dynamics simulations.

Fig. 4 Linear interpolation in internal coordinates (LIIC) between the MECI and the S₁ torsional minimum (left) and relaxed scan from the S₁ torsional minimum to the planar trans-enol conformer (right). Both energy profiles were obtained at TD-PBE0/cc-pVDZ level.

3.3. Non-adiabatic dynamics simulations

3.3.1. Electronic state populations. We begin the discussion of the coupled electron–nuclear dynamics by analyzing the electronic state population evolutions, shown in Fig. 5 for each ensemble separately and for the combined trajectories. The convergence of the electronic populations over time was verified and is reported in Fig. S3. The time evolution of the electronic populations for the full ensemble of trajectories is also reported in terms of the diabatic basis of ππ* and nπ* states at the (cis-)enol FC point in Fig. S4. For the ensemble initialised in the S₂ state, the population of the initially active state decreases rapidly within the first 25 fs, driven by the narrow energy gap between S₁ and S₂ and the accessibility of the MECI near the FC region. A similar mechanism explains the initial depopulation of S₁ in favor of S₂ when the dynamics starts on the S₁ state, although this occurs to a lesser extent. The geometries at which S₁/S₂ hops occur during the first 25 fs predominantly adopt the (cis-)enol conformation and are structurally similar to the MECI located near the FC point (see Fig. S5), confirming that this intersection is readily accessible upon photoexcitation.

Fig. 5 Electronic state population evolutions upon excitation to the S₁ (top) and S₂ (middle) electronic state. Combined S₁ and S₂ ensembles (bottom).

After the first 50 fs, the population transfer reaches an equilibrium in which the electronic populations are roughly constant in time. No population transfer from the excited states to the ground state is observed within the simulation time. It is consistent with the experimentally observed later radiative decay from the excited state to the ground state through fluorescence measurements.^64,65

3.3.2. Typical ESIPT trajectory and key degrees of freedom. The key bond lengths and S₁/S₀ energy difference along a typical reactive trajectory are depicted in Fig. 6. At early times, the molecule evolves within the reactant well, as evidenced by the rapid fluctuations of the donor oxygen–hydrogen distance (blue curve) at a frequency characteristic of the oxygen–hydrogen stretching mode (ca. 3579 cm⁻¹, corresponding to a period of ca. 9.3 fs). In this trajectory, proton transfer, defined by the proton–acceptor oxygen distance (red curve) becoming shorter than the donor oxygen–proton distance, occurs around 50 fs after photo-excitation. This transfer is accompanied by a marked decrease in the ground/excited state energy gap (grey curve), reflecting the large Stokes shift typical of ESIPT chromophores. After the proton transfer, the system resides in the product tautomer well, where the hydrogen–acceptor oxygen distance exhibits the characteristic oxygen–hydrogen bond stretching frequency. The low frequency oscillation observed in the donor oxygen–hydrogen distance likely arises from the in-plane bending motions, involving the carbonyl and hydroxyl functional groups, which periodically drive the hydrogen atom back toward the donor oxygen. The normal modes (calculated at the keto tautomer S₁ equilibrium geometry) associated with this motion have frequencies of 263 and 316 cm⁻¹, corresponding to periods of approximately 106 to 126 fs, in qualitative agreement with the oscillation period observed following proton transfer.

Fig. 6 Key bond lengths and S₁/S₀ energy difference along a typical reactive trajectory observed in the simulations.

Next, to gain statistical insight into the ESIPT process and possible competing pathways, the key internal coordinates are examined for the ensemble of mixed quantum-classical surface-hopping trajectories. Fig. 7 shows the evolution along all trajectories of the donor oxygen–hydrogen bond distance, d(O_d–H), and the hydrogen–acceptor oxygen bond distance, d(H–O_a), which directly describe the proton transfer reaction. The d(O_d–H) distances are initially distributed around the equilibrium value of 0.98 Å at the FC point, but after only 25 fs part of the trajectories exhibits the elongation of the bond distance in favour of the formation of the proton–acceptor oxygen bond, testified by the swarm of trajectories showing the shortening of the d(H–O_a) distance from the average value of 2.00 Å to 0.98 Å. In addition to the shortening of d(H–O_a), we also observe an increase after 50 fs for another portion of the analysed trajectories caused by the torsion of the hydroxyl group around the dihedral angle ϕ₁.

Fig. 7 Time evolution of the key oxygen–hydrogen distances involved in the description of the ESIPT process along all trajectories. The color code indicate the active state.

In nearly all cases where d(O_d–H) increases while d(H–O_a) decreases, leading to the proton remaining bound to the acceptor oxygen, the system resides in the S₁ electronic state. This is consistent with the lower barrier height along the ESIPT coordinate on S₁ (about 0.42 eV lower than on S₂). Consequently, upon excitation to S₂, the system relaxes efficiently to S₁ before proton transfer takes place, in agreement with Nag et al.⁴³ This finding rules out ESIPT occurring directly on the S₂ surface as the origin of the second (slower) ESIPT time constant, at least within the time window accessible to our simulations.

3.3.3. Kinetic modelling. The yield of the enol and keto forms was monitored in time using the key distances described above. At each time step, the molecular structure visited along the trajectory is classified as the enol tautomer if d(O_d–H) ≤ d(H–O_a) and as the keto tautomer if d(O_d–H) > d(H–O_a) distance. Fig. 8 shows the population dynamics of both tautomers obtained from the combined analysis of the two ensembles of trajectories. Proton migration was observed in 69.8% of the trajectories within the simulation time. The convergence of the tautomeric yields over time for each ensemble individually is reported in Fig. S6.

Fig. 8 Time-dependent yields of the reactant (enol, in blue) and product (keto, in red) forms. The black dashed line represents the bi-exponential fit of the reactant and product populations in time.

The time evolution of the enol population was fitted using a bi-exponential function with latency of the form

where δ₁ and δ₂ are the fractions of population decaying with time constants τ₁ and τ₂, respectively, t₀ is the latency time, C is the asymptotic population. The optimised parameters of the bi-exponential fit are reported in Table 1. A reliable description of the reactant decay requires a bi-exponential function rather than a mono-exponential decay. A comparison between mono- and bi-exponential fits is shown in Fig. S7 and Table S1.

Table 1 Optimised parameters from the fit of the enol yield time evolution with a bi-exponential decay function. In parentheses one standard deviation errors on the fit parameters are reported

Fit parameters
t 0 (fs)	δ 1	δ 2	τ 1 (fs)	τ 2 (fs)	C
24.97 (0.20)	0.42 (0.01)	0.58 (0.14)	24.82 (0.84)	510.92 (162.79)	0.00 (0.13)

---

From the bi-exponential fit, two time constants for the ESIPT process were extracted: a fast component on the sub-100 fs timescale (24.82 fs), consistent with previous non-adiabatic dynamics studies,⁴³ and a slower component of approximately 0.51 ps. Experimentally, a second ESIPT time-constant of ca. 5.5 ps was previously found for the same system in methylcyclohexane.²⁴ To obtain reliable estimate of the uncertainty on the predicted slower component, which is beyond the simulated time scale, a bootstrap procedure was employed. From the set of all the trajectories, 200 set made by 100 trajectories were sampled with replacement and from each set we obtained the τ₂ from the fit of the enol and keto yields. The 3rd and 97th percentiles are 0.11 and 1.15 ps, respectively. Despite the large noise, our identification of an additional picosecond phenomenon is in qualitative agreement with previous experiments.²⁴ The results from the non-adiabatic dynamics simulations in the gas phase are compared with measurements in methylcyclohexane. Due to the aprotic non-polar nature of the solvent, only a limited effect on the time scales of the investigated processes is expected. This is confirmed by a comparison between the gas-phase LIIC and the one in methylcyclohexane (Fig. S8). However, subtle solvation effects could explain a quantitative discrepancy on the picosecond time scale.

3.3.4. Origin of the second ESIPT time constant and competing pathway. The enol and keto yields in time from the ensembles of trajectories that started from the excited-states S₁ and S₂ were analysed separately (Fig. S9 and S10 in the SI). For both ensembles, two ESIPT time scales are observed and they are indistinguishable within the error (see Tables S2 and S3), indicating that the two ESIPT time scales do not originate from the different initial electronic state upon excitation.

To understand the origin of the two time constants, we inspected the distributions of the donor oxygen–hydrogen and hydrogen–acceptor oxygen distances, and of the dihedral angles ϕ₁ and ϕ₂ (in Fig. 9 and 10, respectively). The hopping geometries are highlighted with a different colour in both the internal coordinate subspaces. The tails of the distance distributions (i.e., for d(O_d–H), d(H–O_a) > 3.4 Å) in Fig. 9 suggest the formation of the trans isomer of both enol and keto forms. It is also possible to clearly distinguish the hopping geometries into enol-type and keto-type structures from the distance subspace. From Fig. 10 it is possible to appreciate and confirm the exploration of the configurational space associated with the trans form of the enol (i.e., when |ϕ₁| > 130.0 degrees) and of the keto (i.e., when |ϕ₂| > 130.0 degrees). The 130.0 degrees threshold value was selected based on the results of the scan along the hydrogen torsion coordinate (see Fig. 4). At ϕ₁ ≃ 130.0 degrees the system fully overcame the energy barrier to reach the trans-enol conformer.

Fig. 9 Scatter plot of the d(O_d–H) and d(H–O_a) distances visited by the 3-HC along the surface hopping trajectories.

Fig. 10 Scatter plot of the ϕ₁ and ϕ₂ dihedral angles visited by the 3-HC along the surface hopping trajectories.

Statistically, 54.4% of trajectories originating from the (cis-)enol conformer reach the region of the unreactive S₁ torsional minimum (with nπ* character and |ϕ₁| around 21.6 degrees). Among these, 10.2% proceed further to undergo the full (cis-)enol to trans-enol transition (|ϕ₁| > 130.0 degrees), i.e., crossing the barrier connecting the torsional minimum to the trans-enol conformation. These observations indicate that, after crossing the S₁/S₂ seam near the FC point, the system extensively explores the configurational space associated with out-of-plane hydrogen torsion. In some cases, this leads to formation of the trans-enol conformer, representing a competitive pathway to ESIPT. Notably, the amplitude associated with the slower component of the bi-exponential fit is ca. 58.0% which is comparable to the fraction of the trajectories (54.4%) that visit the region of the S₁ torsional minimum, supporting the assignment of the longer time constant to this torsional dynamics.

To test this hypothesis, we analyzed the time-dependent yields of the ESIPT reactant (enol) and product (keto) for two subsets of trajectories. The first subset comprises trajectories in which the |ϕ₁| torsion never exceeds 21.6 degrees, while the second subset includes trajectories in which |ϕ₁| > 21.6 degrees in at least one frame. As expected, the subset without significant torsion along ϕ₁ exhibits rapid enol-to-keto conversion, whereas the subset visiting the S₁ torsional minimum shows a delayed proton transfer (Fig. 11). These results support the idea that out-of-plane hydrogen torsion, readily accessible after crossing the S₁/S₂ seam near the FC region, modulates the ESIPT timescales and gives rise to the second slower proton transfer time constant. To assess whether the observed differences between the two trajectory subsets arise from structural features inherent to the initial positions or momenta, we performed a linear discriminant analysis (LDA)^66,67 on their respective initial conditions. The LDA results indicate a partial separation of the initial conditions along a discriminant axis dominated by the hydroxyl group's out-of-plane motion (Fig. S11).

Fig. 11 Time-dependent yields of the reactant (enol, in blue) and product (keto, in red) forms obtained from (top) the subset of trajectories that do not display a torsion along the |ϕ₁| coordinate of more than 21.6 degrees and (bottom) the subset of trajectories that meet at least in one frame the condition |ϕ₁| ≥ 21.6 degrees.

3.3.5. Global reaction network. The analyses linked the slower ESIPT time constant to out-of-plane hydrogen torsion and visits to the S₁ torsional minimum. The system's excited-state dynamics involve a broad network of conformational changes and tautomer interconversions. To check the dominant collective motions driving the excited-state evolution, we performed a principal component analysis (PCA)⁶² on the full ensemble of trajectories (see SI, Section S11). It confirms that the previously identified critical structures represent the main basins/critical points of the excited-state dynamics. Accordingly, all structures visited along the trajectories were clustered into four categories based on their geometry:

• cis-enol (d(O_d–H) ≤ d(H–O_a) and |ϕ₁| < 21.6 degrees)

• tor-enol (d(O_d–H) ≤ d(H–O_a) and 21.6 degrees ≤|ϕ₁| < 130.0 degrees)

• trans-enol (d(O_d–H) ≤ d(H–O_a) and |ϕ₁| ≥ 130.0 degrees)

• keto (d(O_d–H) > d(H–O_a))

Since each conformer/tautomer can be populated in either the S₁ or S₂ electronic state, this leads to a total of eight distinct states. Every frame of the ensemble of trajectories was assigned to one of these states, and the transitions between them were counted. The resulting transition matrix (with normalised rows) is reported in Fig. S14. From this state-to-state connectivity, we constructed the complete reaction network shown in Fig. 12.

Fig. 12 Reaction network.

It is possible to observe a population transfer between the electronic states S₁ and S₂ when the system reaches the proton transfer product, i.e., the keto species. Indeed, we found several hopping structures and two non-planar MECI that share structural features with the keto tautomer (see Fig. S15). The first keto-type MECI is located ca. 0.18 eV above the keto minimum on the S₁ surface and it results accessible to the system after the proton is transferred. The second keto-type MECI is characterised by the torsion of the hydrogen that was previously transferred and it is located ca. 0.52 eV above the keto minimum on the S₁ surface. Even if less accessible with respect to the previous one, this observation is consistent with the finding of several crossing geometries from the mixed quantum classical simulations adopting a trans-keto conformation.

Finally, this extended picture of the dynamics of the 3-HC system upon excitation to S₁ or S₂ corroborates the presence of a competitive pathway to ESIPT, associated with the out-of-plane motion of the hydrogen atom rather than its direct transfer between the two oxygen atoms.

4. Conclusions

The present work investigated the photodynamics of the 3-HC molecule, with a focus on the excited-state intramolecular proton transfer and the possible competitive relaxation pathways. The coupled electron–nuclear dynamics of 3-HC following excitation to the S₁ and S₂ states were simulated using surface hopping with TD-PBE0 for the electronic structure. Consistent with previous theoretical studies, our simulations indicate that ESIPT takes place on a sub-100 fs timescale. In addition, we directly observed a slower component of the proton transfer dynamics, giving rise to a second time constant, which was not reported in earlier simulations. The second ESIPT timescale predicted by our simulations is in qualitative agreement with experimental observations from time-resolved UV-vis spectroscopy in organic aprotic apolar solvents.

Our analysis of the critical points on the potential energy surfaces, together with the ensemble of mixed quantum-classical trajectories, reveals that the slower ESIPT component originates from competition with the out-of-plane torsional motion of the transferring proton, accessible after crossing the S₁/S₂ seam near the Franck–Condon region. Ultrafast transient spectroscopic techniques may detect the competitive torsional out-of-plane motion corroborating the findings of our non-adiabatic dynamics simulations. Finally, by clustering the visited geometries into representative conformational/tautomeric states and monitoring their interconversion, we constructed an explicit reaction network for 3-HC in the excited state. This network highlights the coexistence of the ESIPT pathway with torsion-mediated channels, offering a unified mechanistic picture that reconciles the ultrafast and slower components of the ESIPT dynamics.

Conflicts of interest

There are no conflicts to declare.

Note added after first publication

This article replaces the previous version published on 5th February 2026 which contained an error in the title. The text of the article remains unchanged.

Data availability

The code for SHARC can be found at https://sharc-md.org/. The version of the code employed for this study is version 3.0.1.

The data supporting this article have been included as part of the supplementary information (SI). The supplementary information provides: the UV-vis absorption spectrum decomposed in terms of the diabatic states; basis set dependence of the LIIC along the ESIPT coordinate; convergence of the electronic state populations with respect to the number of trajectories; diabatic electronic state populations; representation of the enol-type MECI and hopping geometries; convergence of the enol and keto time-dependent yields with respect to the number of trajectories; comparison between the mono- and bi-exponential fit of the yields; influence of the methylcyclohexane solvent on the LIIC along ESIPT coordinate; enol and keto time-dependent yields from S₁ and S₂ ensembles separately; LDA on the initial conditions; PCA on the ensemble of trajectories; transition matrix between conformers/tautomers in their excited states; representation of the keto-type MECIs and hopping geometries. See DOI: https://doi.org/10.1039/d5cp04236d.

Acknowledgements

The authors acknowledge Benjamin Lasorne for reading the original draft of the manuscript and for the stimulating discussion. The simulations in this work were performed using HPC resources from CCIPL (Le center de calcul intensif des Pays de la Loire) and from GENCI-IDRIS (grant 2021-101353). The project is funded by the European Union (ERC, 101040356, ATTOP). Views and opinions expressed are those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

References

1. A. Weller, Naturwissenschaften, 1955, 42, 175–176.
2. A. Weller, Z. Elektrochem., 1956, 60, 1144–1147.
3. J.-E. A. Otterstedt, J. Chem. Phys., 1973, 58, 5716–5725.
4. P. K. Sengupta and M. Kasha, Chem. Phys. Lett., 1979, 68, 382–385.
5. J. Zhao, S. Ji, Y. Chen, H. Guo and P. Yang, Phys. Chem. Chem. Phys., 2012, 14, 8803–8817.
6. A. C. Sedgwick, L. Wu, H.-H. Han, S. D. Bull, X.-P. He, T. D. James, J. L. Sessler, B. Z. Tang, H. Tian and J. Yoon, Chem. Soc. Rev., 2018, 47, 8842–8880.
7. Y. Li, D. Dahal, C. S. Abeywickrama and Y. Pang, ACS Omega, 2021, 6, 6547–6553.
8. M. J. Paterson, M. A. Robb, L. Blancafort and A. D. DeBellis, J. Phys. Chem. A, 2005, 109, 7527–7537.
9. Y. Gong, Z. Wang, S. Zhang, Z. Luo, F. Gao and H. Li, Ind. Eng. Chem. Res., 2016, 55, 5223–5230.
10. J. E. Kwon and S. Y. Park, Adv. Mater., 2011, 23, 3615–3642.
11. K.-C. Tang, M.-J. Chang, T.-Y. Lin, H.-A. Pan, T.-C. Fang, K.-Y. Chen, W.-Y. Hung, Y.-H. Hsu and P.-T. Chou, J. Am. Chem. Soc., 2011, 133, 17738–17745.
12. Z. Zhang, Y.-A. Chen, W.-Y. Hung, W.-F. Tang, Y.-H. Hsu, C.-L. Chen, F.-Y. Meng and P.-T. Chou, Chem. Mater., 2016, 28, 8815–8824.
13. M. Mamada, K. Inada, T. Komino, W. J. Potscavage Jr, H. Nakanotani and C. Adachi, ACS Cent. Sci., 2017, 3, 769–777.
14. B. Li, L. Zhou, H. Cheng, Q. Huang, J. Lan, L. Zhou and J. You, Chem. Sci., 2018, 9, 1213–1220.
15. K. Wu, T. Zhang, Z. Wang, L. Wang, L. Zhan, S. Gong, C. Zhong, Z.-H. Lu, S. Zhang and C. Yang, J. Am. Chem. Soc., 2018, 140, 8877–8886.
16. O. A. Sytina, D. J. Heyes, C. N. Hunter, M. T. Alexandre, I. H. Van Stokkum, R. Van Grondelle and M. L. Groot, Nature, 2008, 456, 1001–1004.
17. P. J. Tonge and S. R. Meech, J. Photochem. Photobiol., A, 2009, 205, 1–11.
18. D. R. Weinberg, C. J. Gagliardi, J. F. Hull, C. F. Murphy, C. A. Kent, B. C. Westlake, A. Paul, D. H. Ess, D. G. McCafferty and T. J. Meyer, Chem. Rev., 2012, 112, 4016–4093.
19. D. Jacquemin, J. Zuniga, A. Requena and J. P. Céron-Carrasco, Acc. Chem. Res., 2014, 47, 2467–2474.
20. H. C. Joshi and L. Antonov, Molecules, 2021, 26, 1475.
21. P. Zhou and K. Han, Acc. Chem. Res., 2018, 51, 1681–1690.
22. C. Fang, R. R. Frontiera, R. Tran and R. A. Mathies, Nature, 2009, 462, 200–204.
23. F. Han, W. Liu and C. Fang, Chem. Phys., 2013, 422, 204–219.
24. K. Chevalier, A. Gruün, A. Stamm, Y. Schmitt, M. Gerhards and R. Diller, J. Phys. Chem. A, 2013, 117, 11233–11245.
25. A. D. Laurent, C. Adamo and D. Jacquemin, Phys. Chem. Chem. Phys., 2014, 16, 14334–14356.
26. Y. Houari, A. Charaf-Eddin, A. D. Laurent, J. Massue, R. Ziessel, G. Ulrich and D. Jacquemin, Phys. Chem. Chem. Phys., 2014, 16, 1319–1321.
27. A. Chrayteh, C. Ewels and D. Jacquemin, Phys. Chem. Chem. Phys., 2020, 22, 854–863.
28. A. J. Aquino, H. Lischka and C. Hättig, J. Phys. Chem. A, 2005, 109, 3201–3208.
29. A. J. Aquino, F. Plasser, M. Barbatti and H. Lischka, Croat. Chem. Acta, 2009, 82, 105–114.
30. A. L. Sobolewski, D. Shemesh and W. Domcke, J. Phys. Chem. A, 2009, 113, 542–550.
31. P. M. Vérité, S. Hédé and D. Jacquemin, Phys. Chem. Chem. Phys., 2019, 21, 17400–17409.
32. A. Chrayteh, C. P. Ewels and D. Jacquemin, Phys. Chem. Chem. Phys., 2020, 22, 25066–25074.
33. A. L. Sobolewski and W. Domcke, Phys. Chem. Chem. Phys., 1999, 1, 3065–3072.
34. Y. Li, L. Wang, X. Guo and J. Zhang, J. Comput. Chem., 2015, 36, 2374–2380.
35. X.-P. Chang, F.-R. Fan, G. Zhao, X. Ma, T.-S. Zhang and B.-B. Xie, Chem. Phys., 2023, 575, 112056.
36. J. M. Ortiz-Sánchez, R. Gelabert, M. Moreno and J. M. Lluch, J. Chem. Phys., 2007, 127, 084318.
37. A. Perveaux, M. Lorphelin, B. Lasorne and D. Lauvergnat, Phys. Chem. Chem. Phys., 2017, 19, 6579–6593.
38. N. Anand, S. V. K. Isukapalli and S. R. Vennapusa, J. Comput. Chem., 2020, 41, 1068–1080.
39. N. Anand, P. Nag, R. K. Kanaparthi and S. R. Vennapusa, Phys. Chem. Chem. Phys., 2020, 22, 8745–8756.
40. S. Pijeau, D. Foster and E. G. Hohenstein, J. Phys. Chem. A, 2017, 121, 4595–4605.
41. S. Pijeau, D. Foster and E. G. Hohenstein, J. Phys. Chem. A, 2018, 122, 5555–5562.
42. L. Spörkel, G. Cui and W. Thiel, J. Phys. Chem. A, 2013, 117, 4574–4583.
43. P. Nag, N. Anand and S. R. Vennapusa, J. Chem. Phys., 2021, 155, 094301.
44. P. Nag and S. R. Vennapusa, J. Photochem. Photobiol., A, 2022, 427, 113767.
45. L. Zhao, X. Geng, J. Wang, Y. Liu, W. Yan, Z. Xu and J. Chen, Phys. Chem. Chem. Phys., 2024, 26, 20490–20499.
46. G. Estiú, J. Rama, A. Pereira, R. E. Cachau and O. N. Ventura, J. Mol. Struct. THEOCHEM, 1999, 487, 221–230.
47. S. Ash, S. P. De, H. Beg and A. Misra, Mol. Simul., 2011, 37, 914–922.
48. C. A. Kenfack, A. S. Klymchenko, G. Duportail, A. Burger and Y. Mély, Phys. Chem. Chem. Phys., 2012, 14, 8910–8918.
49. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
50. T. H. Dunning Jr, J. Chem. Phys., 1989, 90, 1007–1023.
51. F. Lipparini, G. Scalmani, B. Mennucci, E. Cancès, M. Caricato and M. J. Frisch, J. Chem. Phys., 2010, 133, 014106.
52. 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 A.03, Gaussian Inc., Wallingford CT, 2016.
53. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2025, 15, e70019.
54. S. Maeda, K. Ohno and K. Morokuma, J. Chem. Theory Comput., 2010, 6, 1538–1545.
55. S. Mai, P. Marquetand and L. González, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2018, 8, e1370.
56. R. Crespo-Otero and M. Barbatti, Marco Antonio Chaer Nascimento: A Festschrift from Theoretical Chemistry Accounts, Springer, 2012, pp. 89–102.
57. T. Pitesa, S. Polonius, L. González and S. Mai, J. Chem. Theory Comput., 2024, 20, 5609–5634.
58. G. Granucci, M. Persico and A. Toniolo, J. Chem. Phys., 2001, 114, 10608–10615.
59. G. Granucci and M. Persico, J. Chem. Phys., 2007, 126, 134114.
60. B. G. Levine, C. Ko, J. Quenneville and T. J. Martínez, Mol. Phys., 2006, 104, 1039–1051.
61. R. T. McGibbon, K. A. Beauchamp, M. P. Harrigan, C. Klein, J. M. Swails, C. X. Hernández, C. R. Schwantes, L.-P. Wang, T. J. Lane and V. S. Pande, Biophys. J., 2015, 109, 1528–1532.
62. K. Pearson, The London, Edinburgh, and Dublin philosophical magazine and journal of science, 1901, vol. 2, pp. 559–572.
63. I. B. Bersuker, Chem. Rev., 2013, 113, 1351–1390.
64. A. S. Klymchenko and A. P. Demchenko, New J. Chem., 2004, 28, 687–692.
65. L. Giordano, V. V. Shvadchak, N. Arrupe, L. J. F. Lockhart, V. M. Sanchez and T. M. Jovin, Dyes Pigm., 2023, 218, 111479.
66. R. A. Fisher, Ann. Hum. Genet., 1936, 7, 179–188.
67. T. Hastie, R. Tibshirani, J. Friedman, et al., The Elements of Statistical Learning, 2009.

---