Understanding contrasting S₂ → S₁ internal conversion rates in boron-dipyrromethene derivatives via multi-configuration time-dependent hartree method

Abstract

Internal conversion (IC) dynamics from higher-lying to lower electronic states following photoexcitation are often described as ultrafast, occurring in much less than 1 ps. However, IC processes can exhibit a range of time scales, depending critically on the energetic landscape of excited-state manifolds and the strength of vibronic couplings that drive nonadiabatic transitions. These dynamics play a fundamental role in most photochemical and photophysical applications. A previous time-resolved fluorescence study revealed that two structurally analogous boron-dipyrromethene (BODIPY) molecules, PM650 and PM597, exhibit markedly different IC rates, despite both undergoing ultrafast IC in <100 fs from the S₃/S₂ to the S₁ state. Nuclear wave packets persisting in the S₁ state after the IC were also observed by time-resolved fluorescence. To elucidate the origin of these divergent IC rates and the nature of vibronic interactions among excited states, we performed theoretical simulations using the multiconfiguration time-dependent Hartree (MCTDH) method. Our results reveal nonadiabatic decay pathways mediated by vibronically coupled S₁, S₂, and S₃ potential energy surfaces, with multiple conical intersections (CIs) enabling the IC processes. The IC rates obtained from the MCTDH simulations are in good agreement with the experimental observations, including the contrasting rates for PM650 and PM597. Importantly, the proximity of CIs to the Franck–Condon region was found to significantly influence IC efficiency. As more vibrational modes were incorporated into the model, a consistent acceleration of the IC dynamics was observed, underscoring the role of multimode effects in nonadiabatic transitions through CIs. Additionally, coherent vibrational spectra of the S₁ state, generated from nuclear densities computed via MCTDH following excitation to higher states, were found to match experimental results closely, further supporting the conclusions of this study. Overall, these findings advance our understanding of the intricate excited-state dynamics and highlight the critical role of vibronic coupling and CIs in ultrafast IC.

---

1. Introduction

A molecule excited to electronic states higher than S₁ typically undergoes ultrafast relaxation to the S₁ state in less than 1 ps via internal conversion (IC).^1–3 As a result, photoexcitation to a higher-lying state generally leads to emission exclusively from the S₁ state, a phenomenon known as Kasha's rule.⁴ The details of these relaxation processes, such as IC rates and associated reaction coordinates, have significant implications for the photophysical and photochemical properties of molecules and may offer pathways for actively manipulating these properties.

Although IC from S_n (n ≥ 2) to S₁ is generally considered ultrafast, the actual rates can vary significantly, ranging from a few femtoseconds to several picoseconds, depending on the specific molecular electronic couplings.^3,5–16 This process is facilitated by conical intersections (CIs), which provide highly efficient, non-radiative pathways for electronic relaxation.^17,18 It is important to note that nonadiabatic processes via CIs in polyatomic molecules are generally not simple one-dimensional phenomena involving a single reaction coordinate. Rather, they require the inclusion of multiple vibrational modes to accurately model the vibronic interactions.^19–22 Thus, a detailed understanding of the molecular mechanisms underlying ultrafast IC is crucial for both fundamental photophysics and applications involving light–matter interaction.

The IC process is nonadiabatic and inherently involves vibronic interactions, where electronic and vibrational degrees of freedom are strongly coupled. Upon impulsive photoexcitation to an S_n (n ≥ 2) state, coherent nuclear wave packets (NWPs) are formed in the Franck–Condon (FC) region. These NWPs may survive the IC process and persist in the S₁ state.^10,23–27 During the passage through CIs, the amplitudes and phases of the NWPs can be modified. Vibrational modes that are strongly coupled to the IC process may exhibit particularly significant modulation. Therefore, measuring NWPs in excited states can yield valuable information about the reaction coordinates of IC and the vibronic couplings between excited electronic states.

Previously, by employing time-resolved fluorescence (TF) with a time resolution high enough to resolve the motion of NWPs, we reported IC rates and the formation and evolution of NWPs in the S₁ state following photoexcitation to the S₂ and S₃ states for two Boron-dipyrromethene (BODIPY) molecules, PM597 and PM650.¹² The molecular structures are shown in Fig. 1. Despite the structural similarity of the two molecules, the S₂ → S₁ IC rates were vastly different; in PM650 the IC rate is nearly instantaneous, ≫(20 fs)⁻¹, whereas in PM597 it occurs at (51 fs)⁻¹. These two molecules, which share similar molecular structures and electronic characters but exhibit contrasting IC dynamics, provide a valuable opportunity to investigate how subtle changes in molecular and electronic structures and delicate vibronic interactions influence nonadiabatic energy relaxation pathways and dynamics. BODIPY is also an excellent model system owing to its favourable spectroscopic properties, including a large extinction coefficient, high fluorescence quantum yield, narrow absorption and emission linewidths in solution, structural tunability, low sensitivity to solvent polarity, and resistance to aggregation.

Fig. 1 Molecular structures of (a) PM650 and (b) PM597. The bond lengths and bond angles are indicated.

NWPs in the S₁ state following ICs have also been experimentally observed and analysed.¹² In that study, we estimated the displacement between the FC region and the S₁ state, and compared the calculated values with experimental NWPs. However, our analysis was based on the assumption of an impulsive IC process, in which the transition from S₂ to S₁ was treated as instantaneous.²⁸ This approach was partially succeeded in reproducing the experimental observations. A more rigorous approach to calculating NWPs and coherent vibrational spectra (CVS) following photoexcitation and subsequent molecular processes was recently introduced using Born–Oppenheimer molecular dynamics (BOMD) simulations.²⁹ In this method, nuclear dynamics are obtained by projecting the BOMD trajectories onto the normal modes of the product state. However, this molecular-dynamics-based method, despite incorporating on-the-fly quantum calculations of structure and forces, still treats nuclear motions classically. Therefore, a more rigorous quantum mechanical treatment of nuclear motion is essential for a complete understanding of quantum dynamics, particularly for processes involving CIs such as ultrafast IC.

In this study, we adopt a quantum mechanical approach to describe ultrafast IC, explicitly accounting for the vibronic coupling effects. We employed the multi-configuration time-dependent Hartree (MCTDH) method to perform quantum wave packet simulations, enabling tracing of IC pathways and calculation of the NWPs that persist in the S₁ state after IC.^30–35 This quantum nuclear wave packet simulation method not only allows for accurate determination of the IC rates and surviving NWPs, but also provides deeper insight into the role of vibrational modes that are strongly coupled to the IC processes. By capturing the intricate interplay between electronic and nuclear dynamics, our study advances the theoretical framework for understanding ultrafast IC and establishes a foundation for future experimental and computational investigations in photophysics and photochemistry.

2. Methods

2.1. Linear vibronic coupling Hamiltonian

The quantum wave packet propagation, initiated from the electronically excited S_n state, is performed using the MCTDH method within the framework of the linear vibronic coupling (LVC) model. The details of the method have been described extensively in the literature. In brief, a diabatic vibronic Hamiltonian is constructed within the LVC framework to avoid singularities at points where potential energy surfaces (PESs) intersect. The Hamiltonian is expressed as^19,32,36

H = H₀(Q)I₃ + W⁽¹⁾, (1)

where H₀ is the ground-state Hamiltonian, I₃ is the 3 × 3 unit matrix since three electronic states, S₁ to S₃, are considered. Q denotes the mass-weighted dimensionless normal coordinates of the ground state, and W⁽¹⁾ represents the coupling elements between electronic states. The zeroth-order ground-state Hamiltonian includes the diagonal nuclear kinetic energy operator T_N and the harmonic approximation of the diabatic potential energy matrix:where ω_i is the harmonic frequency of the ith vibrational mode. Within the LVC framework, matrix elements of W⁽¹⁾ are given by:where E_m is the vertical excitation energy of the mth state, and κ^m_i, λ^mn_i are the intrastate and interstate vibronic coupling constants, respectively. These are defined as:where V_m and V_n are the adiabatic potential energies of mth and nth electronic states respectively and Q₀ corresponds to the ground-state equilibrium geometry. For nondegenerate electronic state m, κ^m_i is nonzero for totally symmetric vibrational modes. Since the BODIPY molecules studied in this work possess C_s symmetry, κ^m_i is nonzero for the vibrational modes with a′ symmetry across all electronic states. For such modes, λ^mn_i is nonzero when both electronic states m and n share the same symmetry, whereas for modes with a′′ symmetry, λ^mn_i is nonzero only when the electronic states have different symmetries.

The stationary points on the PESs and the positions of minimum energy conical intersections (MECIs) are located within the LVC model using standard expressions reported in the literature. The energy and position of MECI along the coupling modes is given by:¹⁹

andrespectively, whereand f is the number of totally symmetric vibrational modes.

2.2. The MCTDH wave packet propagation method

The details of the MCTDH simulation method are well described in the literature.^30,31,37 In the MCTDH approach, the wavefunction used to solve the time-dependent Schrödinger equation for a system with f degrees of freedom (DOFs) is represented as a sum of products of single-particle functions (SPFs):where Q₁,…,Q_f are the nuclear coordinates, n_k is the number of SPFs used to represent the motion along the kth DOF, and is the SPF for the kth DOF. The SPFs themselves are expanded in terms of time-independent primitive basis functions:where N_k is the number of grid points (or primitive basis functions) for the kth DOF. In this study, the harmonic oscillator discrete variable representation was employed as the primitive basis.

For the MCTDH simulations, a subset of 33 vibrational modes was selected out of 114 and 174 modes for PM650 and PM597, respectively, based on their vibronic coupling strength (κ²/2ω²). The initial ground-state vibrational wave packet was vertically promoted to the S₂ or S₃ electronic state and propagated for 1 ps. To analyse the wave packet dynamics, nuclear density distributions and diabatic electronic state populations were evaluated. Based on the nuclear density projected onto the S₁ state following excitation to S₂/S₃, CVS were computed.

Wave packet simulations were performed using the Heidelberg MCTDH package (version 8.4, Revision 8).³⁸ All simulation parameters were obtained from quantum chemical calculations in the gas phase using the Gaussian 16 software package.³⁹ Ground-state geometry optimizations and frequency calculations for PM650 and PM597 were conducted under C_s symmetry using density functional theory (DFT) with the CAM-B3LYP functional and the 6-311++G(d,p) basis set. Excited-state properties, including vertical excitation energies and oscillator strengths of low-lying electronic states, were computed using time-dependent DFT (TDDFT) with the same functional and basis set. Transition electron density analysis was carried out using Multiwfn software, version 3.8.^40,41 LVC parameters and details of the normal modes for the MCTDH simulations are listed in the SI Tables S2–S11.

3. Results

3.1. Vertical excitation energies

Vertical excitation energies and oscillator strengths calculated using TDDFT at the CAM-B3LYP/6-311++G(d,p) level are summarized in Table 1. The experimental absorption spectra in cyclohexane exhibit two bands at 2.08 and 2.88 eV for PM650, and at 2.30 and 3.18 eV for PM597. The computed energy gaps between the S₁ and S₂ states are 0.76 eV for PM650 and 0.92 eV for PM597, in good agreement with the experimental values. Alternative DFT methods show less accurate agreement for the S₁–S₂ energy gap, as shown in the SI Table S1, although some functionals yield better estimates of absolute excitation energies.

Table 1 Vertical excitation energies of the singlet states of PM650 and PM597. The numbers in parentheses indicate oscillator strengths

States	PM650 (eV)	PM597 (eV)
S1 (A′′)	2.65 (0.54)	2.89 (0.70)
S2 (A′′)	3.41 (0.08)	3.81 (0.09)
S3 (A′)	3.64 (0.04)	4.05 (0.05)
S4 (A′)	4.99 (0.26)	4.93 (0.00)

---

The transition density isosurfaces corresponding to these excitations are shown in the SI (Fig. S1). For both dyes, the S₀ → S₁ transition is strongly delocalized over the BODIPY π-system, indicating a bright π → π* excitation with a large oscillator strength. In contrast, the S₀ → S₂ transition exhibits more fragmented and localized lobes, suggesting mixed orbital contributions consistent with a less allowed excitation. The S₀ → S₃ transition shows a similar degree of density localization and is less delocalized than S₁. Although both dyes display qualitatively similar features in their transition density isosurfaces, the transition density of PM597 is more strongly aligned and appears more coherent, thereby strengthening the transition dipole and yielding the larger oscillator strength values reported. Because S₂ and S₃ are closely spaced in energy and have comparable oscillator strengths, both states are expected to be populated upon photoexcitation. Consequently, three excited states, S₁ to S₃ were included in the MCTDH simulations.

3.2. Potential energy surfaces

To identify the vibrational modes that play a dominant role in IC dynamics, we first plotted the adiabatic potential energy curves of S₁, S₂, and S₃ along all totally symmetric vibrational modes for both PM650 and PM597. From these curves, we selected a few modes that exhibited curve crossing near to the Franck–Condon geometry. Such crossings are indicative of accessible CI that can strongly facilitate ultrafast IC. We have identified several normal modes (promoting modes) that are likely to play significant roles in the ultrafast IC processes. They are depicted in Fig. 2, which are all in-plane vibrations. These modes fall within the 1200–1600 cm⁻¹ range, indicating that they primarily involve skeletal stretching vibrations and C–H rocking and scissoring motions. In PM650, the dominant displacements involve skeletal stretching of the cyclic core, whereas in PM597, they include C–H scissoring motions of the tert-butyl substituents as well. Fig. 3 and 4 show one-dimensional cuts of the PESs of the S₁, S₂, and S₃ states along these vibrational coordinates. Distortions along these normal modes cause the PESs to converge, leading to curve crossings. The coordinates of the crossings, denoted as (Q_d, V_d), are provided in parentheses.

Fig. 2 Normal modes of (a)–(c) PM650 and (d)–(e) PM597.

Fig. 3 Adiabatic PESs of the S₁ (red), S₂ (green), and S₃ (blue) states along the (a) Q₆₇, (b) Q₇₅, and (c) Q₉₅ vibrational modes of PM650. The small circles represent the calculated ab initio energies at displaced geometries, and the large open circles indicate curve crossing points.

Fig. 4 Adiabatic PESs of the S₁ (red), S₂ (green), and S₃ (blue) states along the (a) Q₉₆ and (b) Q₁₁₀ vibrational modes of PM597. The small circles represent the calculated ab initio energies at displaced geometries, and the large open circles indicate curve crossing points.

Curve crossings between the S₂ and S₁ states in PM650 are observed at (−4.68, 4.67 eV) and (5.18, 5.00 eV) along the Q₆₇ and Q₇₅ vibrational modes, respectively. Additionally, an S₃/S₂ crossing occurs at (−5.90, 6.02 eV) along the Q₉₅ mode. In PM597, S₂/S₁ curve crossings are located at (6.20, 6.27 eV) and (8.75, 9.68 eV) along the Q₉₆ and Q₁₁₀ modes, respectively, with an S₃/S₂ crossing at (−6.80, 8.07 eV) also along Q₉₆. These curve crossings give rise to the well-known multimode effect, which leads to the formation of multidimensional CIs.

A more detailed analysis of CIs was conducted in the multimode context using the LVC potentials. This includes identification of electronic state minima and the locations of MECIs in the space of promoting modes, as shown in Fig. 5. Corresponding values are summarized in Table 2. In PM650, the FC point on the S₂ state lies at 3.41 eV with a stabilization energy of 0.44 eV, and the S₂/S₁ MECI is located at 2.97 eV. The energy difference between the S₂ minimum (S_2,min) and the MECI is practically zero, indicating that relaxation along the S₂ surface inherently involves passing through the MECI. For PM597, the energy difference between the S_2,min and the MECI is 0.05 eV (400 cm⁻¹), suggesting that the system should overcome a modest barrier to reach the CI.

Fig. 5 Schematic of the energetics of FC point on S₂, PES minima and MECI. Location of S₂/S₁ MECI marked with the dashed line.

Table 2 Energies (in eV) of the various points in S₁, S₂, and S₃ states of PM650 and PM597

	S2,FC	S3,FC	S2/S1 MECI	S3/S2 MECI	S1,min	S2,min	S3,min
PM650	3.41	3.64	2.97	3.22	2.57	2.97	3.22
PM597	3.81	4.05	3.47	4.43	2.82	3.42	3.71

---

The geometries of the MECIs can be approximated by their positions along the promoting modes; Q₆₇ and Q₇₅ in PM650, and Q₉₆ in PM597. For PM650, the coordinate of the S₂/S₁ MECI along Q₆₇ (Q²¹_CI,67) is −0.57 and Q²¹_CI,75 = 0.84. In contrast, for PM597, Q²¹_CI,96 = 0.91, indicating that a larger distortion from the FC geometry is required to access the MECI. Thus, the MECI for PM650 is located closer to the FC point, facilitating more efficient IC, whereas in PM597, the MECI lies farther from the FC region, possibly slowing down the IC process.

The geometric parameter Q^mn_CI,i provides insights into the degree of distortion along the normal modes from the FC point to the MECI. To evaluate the associated nuclear configurations, distortions corresponding to the Q^mn_CI,i values were applied to the vibrational modes using the same procedure employed in determining the vibronic coupling constants. For the vibrational modes Q₆₇ and Q₇₅ in PM650 and Q₉₆ in PM597, we constructed a graph in Fig. 6 that shows the percentage difference of relevant structural parameters. The specific bond lengths and angles used in this analysis are indicated in Fig. 1. In PM650, Q²¹_CI,75 shows substantial deviations (>1%) in r₂, r₈, and a₅ with notable changes also observed in r₆ and a₉. For Q²¹_CI,67 significant deviations are found in r₃ and a₄. For PM597 Q²¹_CI,96 exhibits prominent structural changes in r₂, r₃, and a₄.

Fig. 6 Deviations in the structural parameters of (a) Q²¹_CI,75 and (b) Q²¹_CI,67 in relation to the FC geometry of PM650. (c) Deviations in the structural parameters of Q²¹_CI,96 in relation to the FC geometry of PM597.

3.3. MCTDH calculation of the nuclear wave packet dynamics

We performed MCTDH simulations within the LVC framework to investigate the excited-state dynamics and NWPs persisting in the S₁ state following ICs. In these simulations, the ground-state vibrational wave packet is vertically excited to either the S₂ (A′′) or S₃ (A′) singlet states of PM650 and PM597 including 33 selected vibrational modes.

The diabatic electronic populations of the S₁–S₃ states following excitation to S₂ or S₃ in PM650 are shown in Fig. 7. Upon excitation to S₂, PM650 exhibits ultrafast population transfer to the S₁ state in less than 10 fs, accompanied by a minor population in S₃. By ∼10 fs, over half of the initial S₂ population is redistributed to S₁ and S₃. Exponential fit of the population of the S₁ state gives a rise time constant of 9.6 fs (Fig. S2). Experimental and calculated time constants for the S₂ → S₁ ICs in PM650 and PM597 upon excitation to S₂ are summarized in Table 3. Since the MCTDH method does not incorporate dissipative mechanisms, a dynamic equilibrium is established between S₂ and S₁ populations, leading to a persistent residual population in S₂. In actual systems, however, this remaining S₂ population is expected to decay rapidly to S₁via dissipation. The proximity of the CI to the FC point and the near-degeneracy between the S_2,min and the S₂/S₁ MECI provide favourable conditions for such ultrafast wave packet transfer well within a single vibrational period along the relevant normal mode. Fig. 7(b) displays the population dynamics following excitation to the S₃ state. Due to the negligible energy gap between the S_3,min and the S₃/S₂ MECI, the wave packet on S₃ quickly transfers to S₂ in less than 10 fs and subsequently to S₁ by ∼20 fs. This ultrafast IC from both S₂ and S₃ to S₁ is consistent with the experimental observation.¹²

Fig. 7 Diabatic electronic populations of the S₁, S₂ and S₃ states including 33 vibrational modes of PM650 following initial excitations to (a) S₂ and (b) S₃ states. Panel (c) shows the populations of the S₁, S₂ and S₃ states including 23 vibrational modes following initial excitation to S₂ state.

Table 3 Experimental and calculated time constants for the S₂ → S₁ ICs in PM650 and PM597

PM650		PM597
Exp. (fs)	Calc. (fs)	Exp. (fs)	Calc. (fs)
<20	9.6	51	85

---

To investigate the collective role of vibrational modes, we performed additional MCTDH simulations with various combinations of vibrational modes. One such case is shown in Fig. 7(c), where nine specific vibrational modes at frequencies of 31.5, 35.5, 40.3, 294.4, 1228.4, 1266.3, 1451.0, 1522.0, and 1629.2 cm⁻¹ are selectively excluded from the simulation. A marked slowdown in population decay was observed compared to that in Fig. 7(a). Notably, the exclusion of the 1522.0 cm⁻¹ mode, along with the others, led to a significantly reduced rate of population transfer. This gradual decay highlights the critical role of the collective contribution of multiple vibrational modes to the ultrafast IC in polyatomic molecules.

The population dynamics shown in Fig. 7(a) and (b) exhibit oscillatory features indicative of coherent nuclear motions. The population traces were fit using a tri-exponential function, and the residual was Fourier transformed to obtain the spectra shown in Fig. 8. Distinct peaks appear at 1242 and 1438 cm⁻¹, corresponding to the vibrational modes of the ν₆₇ and ν₇₅ modes, respectively. This suggests that these modes play a significant role in the IC reaction coordinates and contribute strongly to the vibronic coupling between the S₁ and S₂ states.

Fig. 8 Fourier transform (FT) power spectra of the oscillation part of the populations of (a) S₁ and (b) S₂ states obtained from the MCTDH simulation following excitation to S₂ state.

The analysis of specific vibrational modes involved in the early stages of excited-state dynamics can be performed by tracking the temporal evolution of nuclear density distributions obtained from MCTDH simulations. As an example, Fig. S3 shows the nuclear distributions on the S₁–S₃ states along the Q₆₇ and Q₇₅ vibrational modes at various time intervals following vertical excitation to the S₂ state. As early as 5 fs, substantial nuclear density is already observed in the S₁ state, along with some nuclear density in the S₃ state. By 20 fs, the nuclear density is distributed mostly in S₁ state, as reflected in the population dynamics in Fig. 7. Since the vibrational periods of the Q₆₇ and Q₇₅ modes are 26 and 23 fs, respectively, the initial nuclear wave packet has traversed nearly a full circle in this 20 fs timeframe.

The diabatic electronic populations of the S₁–S₃ states following excitation to S₂ or S₃ state in PM597 are shown in Fig. 9. In contrast to the nearly instantaneous depopulation of the S₂ state in PM650, the decay of the wave packet in PM597 is substantially slower. The rise time of the S₁ population and the decay time of the S₂ population are both approximately 85 fs (Fig. S2), which is much slower than that in PM650. As discussed in Section 3.2, this difference may be attributed to variations in the energetics and spatial location of the S₂/S₁ MECI within the multidimensional coordinate space. The S₂/S₁ MECI in PM650 lies 0.44 eV below the FC point, while in PM597, it lies 0.3362 eV below. In addition, in PM597, the S₂/S₁ MECI lies slightly above the S_2,min by 0.05 eV, which could also contribute to the slower IC.

Fig. 9 Diabatic electronic populations of S₁, S₂ and S₃ states including 33 vibrational modes of PM597 after excitation to (a) S₂ and (b) S₃ states. Panel (c) shows the populations of the S₁, S₂ and S₃ states including 28 vibrational modes following initial excitation to S₂ state.

When PM597 is excited to the S₃ state, the S₃ population undergoes ultrafast decay to S₂ within ∼10 fs, which is similar to the S₃ → S₂ IC time in PM650, followed by a much slower transition to S₁ on a timescale similar to that of direct S₂ excitation. Interestingly, the S₃ → S₂ IC rates of PM650 and PM597 are comparable, even though the S₃/S₂ MECI lies above both the S₃ FC point and the S_3,min, whereas it lies below in PM650.

We also performed additional MCTDH simulations for PM597 with various combinations of vibrational modes. One such example is shown in Fig. 9(c), where five vibrational modes with frequencies of 1007, 1147.8, 1271.9, 1427, and 1442.1 cm⁻¹ are omitted deliberately. As a result, the IC dynamics for both S₂ → S₁ and S₃ → S₂ transitions slowed greatly, with little progress observed even after 1 ps.

The temporal evolution of nuclear density distributions for PM597 was also obtained from MCTDH simulations. Fig. S4 presents a series of snapshots of the wave packet density distributions along the Q₉₆ and Q₁₁₀ vibrational modes following excitation to the S₂ state. At 25 fs, the nuclear density in the S₂ state shows minimal change, reflecting the much slower IC dynamics in PM597. The density gradually shifts to the S₁ state over approximately 100 fs, as indicated by the population dynamics in Fig. 9. Since the vibrational periods of the Q₉₆ and Q₁₁₀ modes are 26 and 23 fs, respectively, the nuclear wave packet in the S₂ state at 25 fs is returned close to its initial position (center) of this two-dimensional space.

3.4. Coherent vibrational spectra

NWPs in the S₁ state following photoexcitation to the S₂ and S₃ states of PM597 and PM650 have been previously reported using TF spectroscopy.¹² The frequency range was limited to approximately 600 cm⁻¹ in that work due to the finite time resolution of the experiment. Since the initial NWPs in S₂ or S₃ state evolve on excited-state PESs and survive the ICs through CIs, the consequential NWPs in S₁ carry valuable information about the PESs and may provide insights into vibrational dynamics, relaxation pathways, and vibronic interactions in the excited state.^{10,23–27,42}

We calculated the CVS of the S₁ state following excitation to the S₂ or S₃ states, using nuclear densities obtained from MCTDH simulations. The oscillation amplitude of a vibrational mode in a TF signal is be proportional to the amplitude of the NWP oscillation multiplied by the vibrational reorganization energy of the corresponding mode in S₁.^28,43 To compare with the experimental CVS from TF spectroscopy, we first calculated the first moment of the nuclear density in the S₁ state as a function of time for each normal mode i, denoted as Q̄_i. The amplitude of the ith mode in the CVS is then given by:

I_i ∝ (Q̄_i,max − Q̄_i,min)·ω_id̃_i = (Q̄_i,max − Q̄_i,min)·κ¹_i, (12)

where d̃_i is the displacement of the ith mode between the ground and S₁ states.

The experimental and theoretical CVS results are shown in Fig. 10 and 11 for PM650 and PM597, respectively, with the vibrational assignments provided in Tables 4 and 5. Calculated CVS up to 1700 cm⁻¹ are shown in Fig. S5 and S6 for PM650 and PM597, respectively. For PM650, all major peaks in the experimental CVS are well reproduced by the MCTDH-calculated spectra, whereas the agreement is less satisfactory for PM597. Nonetheless, the overall good agreement between the calculated and experimental CVS corroborates the MCTDH simulations presented in this work.

Fig. 10 (a) Experimental CVS of PM650 acquired by linear prediction singular value decomposition method of the TF of PM650 in cyclohexane following photoexcitation to S₂ (and S₃) states. For comparison with the calculated CVS, oscillation amplitudes were corrected for the finite experimental time-resolution. Each peak is a Lorentzian with a linewidth of 5 cm⁻¹. (b) Calculated CVS from MCTDH simulation following excitation to S₂ state. Each peak is a Gaussian with a linewidth of 5 cm⁻¹.

Fig. 11 (a) Experimental CVS of PM597 acquired by linear prediction singular value decomposition method of the TF of PM597 in cyclohexane following photoexcitation to S₂ (and S₃) states. For comparison with the calculated CVS, oscillation amplitudes were corrected for the finite experimental time-resolution. Each peak is a Lorentzian with a linewidth of 5 cm⁻¹. (b) Calculated CVS from MCTDH simulation following excitation to S₂ state. Each peak is a Gaussian with a linewidth of 5 cm⁻¹.

Table 4 Peak assignments of the CVS of PM650

Amplitude	Exp. ω (cm^−1)	Calc. ω (cm^−1)	Mode
2.89	153.3	151.6	ν 13
2.62	243.2	250.8	ν 18
2.09	292.0	292.0	ν 22
0.97	392.6	386.3	ν 28
7.06	412.2	427.5	ν 29
2.82	547.7	548.5	ν 33

---

Table 5 Peak assignments of the CVS of PM597

Amplitude	Exp. ω (cm^−1)	Calc. ω (cm^−1)	Mode
1.14	94.5	76.6	ν 5
0.88	140.5	141.1	ν 9
7.06	257.1	271.8	ν 23
2.82	489.1	510.6	ν 47

---

4. Discussion

Although both BODIPY molecules exhibit ultrafast IC from the S₂/S₃ states to S₁ in much less than 100 fs, their IC rates differ significantly. The MCTDH simulations performed in this work successfully reproduce the experimentally observed IC rates. In general, IC rates are determined by the vibronic coupling constants between electronic states and the energetic landscape of the excited-state manifolds. The average magnitudes of the intrastate and interstate coupling constants for the two molecules are summarized in Table 6. The coupling constants for PM650 are consistently, though only slightly, larger than those of PM597. For example, the S₁/S₂ interstate coupling constants are 0.0969 eV and 0.0907 eV for PM650 and PM597, respectively. While stronger electronic couplings in PM650 may contribute to faster nonadiabatic transitions, the differences in coupling alone may not be sufficient to account for the large disparity in the IC rates.

Table 6 Average values of the intrastate (κ) and interstate (λ) coupling constants of PM650 and PM597. The numbers in superscripts indicates electronic states

Couplings	PM650 (eV)	PM597 (eV)
κ ^1	0.0099 ± 0.014	0.0069 ± 0.008
κ ^2	0.0257 ± 0.035	0.0214 ± 0.026
κ ^3	0.0247 ± 0.032	0.0202 ± 0.023
λ ^12	0.0969 ± 0.066	0.0907 ± 0.068
λ ^13	0.0498 ± 0.046	0.0364 ± 0.034
λ ^23	0.0296 ± 0.032	0.0193 ± 0.023

---

One of the main reasons for the large difference in IC rates may originate from the different energetic landscapes along the vibrational modes, in particular the promoting modes. The relative energy of S₂/S₁ MECI with respect to the S_2,min in PM650 is nearly degenerate, whereas S₂/S₁ MECI is located 400 cm⁻¹ higher than S_2,min in PM597. For PM650, the negligible energy barrier allows the wave packet to access the MECI shortly after the photoexcitation is initiated. In addition, the S₃/S₂ MECI and S_3,min in PM650 are nearly degenerate, whereas S₃/S₂ MECI is located 5800 cm⁻¹ higher than S_3,min in PM597. In particular, the S₃/S₂ MECI is located even higher than the S₃ FC point. Considering that the diabatic populations in S₃ are significant and that S₃/S₁ interstate couplings are significant, the S₃/S₂ MECI may also affect the overall IC rates to the S₁ state.

Although the coupling constants and the energetics contributes to some extent, the substantial structural deformation inherent to the Q^mn_CI,i configuration should also contribute to the difference in IC rates for the two molecules. In PM650, the S₂/S₁ MECI is located close to the FC point, while in PM597, it is more distant. PM650 requires minimum structural reorganization to access MECI due to the close positioning of MECI to FC point. In contrast, PM597 requires more extensive distortion as the MECI located farther away from the FC point. It is important to note that in this structural relaxation, the multimode effect is critical to reach the CI in this nonadiabatic process. The multimode effect is clearly demonstrated in the simulation with reduced DOFs, where the IC is decelerated significantly in both molecules.

All vibrational modes in PM650 that show large amplitude of motion in S₁ state have strong intrastate coupling to S₂ except the 250.8 (ν₁₈) cm⁻¹ mode. This is natural as these modes are excited strongly by the initial photoexcitation to the S₂ state. The ν₁₈ mode, however, is not activated by the photoexcitation due to its negligible intrastate coupling, but appear strongly in the experimental CVS indicating that it is activated in the course of the IC process from S₂ to S₁. Because the IC process is essentially impulsive for all the vibrational modes except the high frequency CH stretching modes, any vibrational modes that has large displacement between the ground (FC) state and S₁ states should be excited strongly. The ν₁₈ mode indeed has a large displacement (large κ¹) between the ground and S₁ states. Therefore, the CVS calculated simply from the displacements matches well with the experiment in PM650. In contrast, for PM597, where the IC process is slow and therefore not impulsive, the CVS calculated from the displacement does not match well, because details of the IC process influence the formation of the NWPs in S₁.

Interestingly, the S₁ and S₂ state populations in PM650 exhibit oscillatory behaviour driven by the in-plane promoting modes ν₆₇ and ν₇₅, which are associated with the S₂/S₁ CIs. Due to strong vibronic coupling along these coordinates, the wave packets initially populated in the S₂ state propagate to the S₁ state through the CIs, and a portion of them subsequently return to the S₂ state. This clearly indicates that the S₁ and S₂ PESs are strongly coupled along these two modes, which play a significant role in the ultrafast IC. In contrast, population oscillations are barely observed in PM597, suggesting that there are no dominant vibrational modes facilitating the fast IC. Instead, the IC is likely driven by small contributions from many modes.

5. Conclusions

In this work, MCTDH quantum dynamics simulations were used to investigate the evolution of nuclear wave packets and ultrafast IC in BODIPY molecules following excitation to electronic states higher than S₁. In PM650, rapid population transfer occurs through CIs along a few key vibrational modes, whereas PM597 exhibits significantly slower population transfer. Notably, nuclear wave packet coherence persists in the S₁ state after ultrafast IC, as reflected in the coherent vibrational spectra. This computational study provides a comprehensive understanding of the IC dynamics of BODIPY derivatives PM650 and PM597, with a focus on mechanistic details including MECIs, nuclear wave packet dynamics, and other factors that influence IC rates.

The striking difference in IC rates between these structurally similar derivatives is attributed to several computationally identified factors. In PM650, the MECI lies close to the Franck–Condon region, with an almost negligible energy gap between the S₂ minimum and the S₂/S₁ MECI. This allows the nuclear wave packet to rapidly access the MECI following photoexcitation, resulting in the ultrafast IC from S₂ to S₁. Most of the initial population in S₂ transfers to the S₁ state within 10 fs. This ultrafast IC process involves only small structural rearrangement. Furthermore, PM650 exhibits strong vibronic coupling along a couple of dominant vibrational modes, providing efficient channels for population transfer. In contrast, PM597 shows a much slower relaxation process from S₂ to S₁ in approximately 86 fs, substantially slower than that of PM650. This slower decay is attributed to the remote location of the MECI from the FC point and a larger energy gap between the S₂ minimum and the S₂/S₁ MECI. Accessing the MECI in PM597 requires significant structural distortion, and the molecule exhibits weaker vibronic coupling between excited states.

These results extend our understanding of the photophysics of molecules excited to electronic states above S₁. Identifying vibrational modes that play significant roles in IC may enable fine-tuning of the IC process by modulating vibronic coupling along these key modes.

Author contributions

Neethu Anand: conceptualization, methodology, validation, investigation, data curation, writing – original draft, writing – review & editing, visualization. Munnyon Kim: conceptualization, methodology, data curation, writing – review & editing. Changmin Lee: conceptualization, methodology, data curation, writing – review & editing. Jinhyuk Ma: conceptualization, methodology, data curation, writing – review & editing. T. Joo: conceptualization, methodology, writing – review & editing, visualization, supervision, project administration, funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request. Input parameters and operator files for the MCTDH simulations are available in the GitHub repositories (https://github.com/AnandNeethu/BODIPY_MCTDH_DYNAMICS).

Supplementary information (SI): Contains isosurfaces of transition densities for low-lying excited-states, exponential fits of the populations of S₁ state following excitation to S₂ state, nuclear density contours, CVS calculated from MCTDH simulations, vertical excitation energies of PM650 using different DFT methods, and LVC parameters for the MCTDH simulations. See DOI: https://doi.org/10.1039/d5cp02771c.

Acknowledgements

This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (no. RS-2020-NR049542).

Notes and references

1. W. Domcke and D. R. Yarkony, Annu. Rev. Phys. Chem., 2012, 63, 325–352.
2. J. C. del Valle and J. Catalán, Phys. Chem. Chem. Phys., 2019, 21, 10061–10069.
3. W. K. Peters, D. E. Couch, B. Mignolet, X. Shi, Q. L. Nguyen, R. C. Fortenberry, H. B. Schlegel, F. Remacle, H. C. Kapteyn, M. M. Murnane and W. Li, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, E11072–E11081.
4. M. Kasha, Discuss. Faraday Soc., 1950, 9, 14–19.
5. M. Beer and H. C. Longuet-Higgins, J. Chem. Phys., 1955, 23, 1390–1391.
6. S. Murata, C. Iwanaga, T. Toda and H. Kokubun, Chem. Phys. Lett., 1972, 13, 101–104.
7. H. Chosrowjan, S. Tanigichi, T. Okada, S. Takagi, T. Arai and K. Tokumam, Chem. Phys. Lett., 1995, 242, 644–649.
8. G. G. Gurzadyan, T. H. Tran-Thi and T. Gustavsson, J. Chem. Phys., 1998, 108, 385–388.
9. T. Horio, T. Fuji, Y.-I. Suzuki and T. Suzuki, J. Am. Chem. Soc., 2009, 131, 10392–10393.
10. T. S. Kuhlman, T. I. Sølling and K. B. Møller, ChemPhysChem, 2012, 13, 820–827.
11. J. Sung, P. Kim, S. Saga, S.-Y. Hayashi, A. Osuka and D. Kim, Angew. Chem., Int. Ed., 2013, 52, 12632–12635.
12. C. Lee, K. Seo, M. Kim and T. Joo, Phys. Chem. Chem. Phys., 2021, 23, 25200–25209.
13. D. W. Cho, M. Fujitsuka, J. H. Ryu, M. H. Lee, H. K. Kim, T. Majima and C. Im, Chem. Commun., 2012, 48, 3424–3426.
14. G. J. Hedley, A. Ruseckas, A. Harriman and I. D. Samuel, Angew. Chem., Int. Ed., 2011, 50, 6634–6637.
15. P. Toele, H. Zhang, C. Trieflinger, J. Daub and M. Glasbeek, Chem. Phys. Lett., 2003, 368, 66–75.
16. T. Fujino, S. Y. Arzhantsev and T. Tahara, J. Phys. Chem. A, 2001, 105, 8123–8129.
17. G. Herzberg and H. C. Longuet-Higgins, Discuss. Faraday Soc., 1963, 35, 77–82.
18. M. S. Schuurman and A. Stolow, Annu. Rev. Phys. Chem., 2018, 69, 427–450.
19. H. Köppel, W. Domcke and L. S. Cederbaum, in Adv. Chem. Phys., ed. I. Prigogine and S. A. Rice, John Wiley & Sons, Inc., 1984, vol. LVII, pp. 59–246.
20. D. R. Yarkony, Rev. Mod. Phys., 1996, 68, 985–1013.
21. F. Bernardi, M. Olivucci and M. A. Robb, Chem. Soc. Rev., 1996, 25, 321–328.
22. B. G. Levine and T. J. Martínez, Annu. Rev. Phys. Chem., 2007, 58, 613–634.
23. T. Fuji, H. J. Ong and T. Kobayashi, Chem. Phys. Lett., 2003, 380, 135–140.
24. Y.-I. Suzuki, T. Fuji, T. Horio and T. Suzuki, J. Chem. Phys., 2010, 132, 174302–174308.
25. J. L. P. Lustres, A. L. Dobryakov, A. Holzwarth and M. Veiga, Angew. Chem., Int. Ed., 2007, 46, 3758–3761.
26. M. Liebel, C. Schnedermann and P. Kukura, Phys. Rev. Lett., 2014, 112, 198302.
27. S. Y. Kim and T. Joo, J. Phys. Chem. Lett., 2015, 6, 2993–2998.
28. C. H. Kim and T. Joo, Phys. Chem. Chem. Phys., 2009, 11, 10266–10269.
29. W. Heo, C. Lee, S. H. Sohn and T. Joo, J. Chem. Phys., 2025, 162, 154108.
30. H. D. Meyer, U. Manthe and L. S. Cederbaum, Chem. Phys. Lett., 1990, 165, 73–78.
31. M. H. Beck, A. Jäckle, G. A. Worth and H. D. Meyer, Phys. Rep., 2000, 324, 1–105.
32. N. Anand, K. Welke, S. Irle and S. R. Vennapusa, J. Chem. Phys., 2019, 151, 214304.
33. S. Mandal and C. Daniel, J. Phys. Chem. A, 2024, 128, 3126–3136.
34. J. A. Green, D. Brey, L. P. Razgatlioglu, B. Ali, B. Błasiak and I. Burghardt, J. Chem. Theory Comput., 2024, 20, 8127–8139.
35. U. Manthe, H. D. Meyer and L. S. Cederbaum, J. Chem. Phys., 1992, 97, 3199–3213.
36. L. S. Cederbaum, H. Köppel and W. Domcke, Int. J. Quantum Chem., 1981, 20, 251–267.
37. H. Wang and H.-D. Meyer, J. Chem. Phys., 2018, 149, 044119.
38. G. A. Worth, M. H. Beck, A. Jäckle and H. D. Meyer, The MCTDH Package, Version 8.4, 2007.
39. 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 CT, 2016.
40. T. Lu and F. Chen, J. Comput. Chem., 2012, 33, 580–592.
41. T. Lu, J. Chem. Phys., 2024, 161, 082503.
42. D. Picconi, Photochem. Photobiol. Sci., 2021, 20, 1455–1473.
43. G. Lee, J. Kim, S. Y. Kim, D. E. Kim and T. Joo, ChemPhysChem, 2017, 18, 670–676.

---