Beyond the three-state picture: when higher-lying excited states become quantitatively indispensable

Abstract

Higher-lying excited states beyond S₁ and T₁ are widely recognized in many photophysical systems, including thermally activated delayed fluorescence (TADF). However, their explicit and quantitative impact on photophysical observables such as photoluminescence quantum yields (PLQY) and lifetimes is difficult to be attained experimentally and it has not been systematically assessed within a fully ab initio kinetic modeling framework. To address this gap, we developed KinLuv, a multistate excited state kinetic model that includes higher-lying excited states (S₂, T₂) and all possible monomolecular interconversion processes between all the electronic states, whose rate constants were computed using Fermi's golden rule (FGR) explicitly including the Herzberg–Teller (HT) vibronic coupling effect. We applied KinLuv to prototypical multi-resonance TADF (MR-TADF) emitters and their derivatives, as well as other representative organic chromophores, demonstrating its broad applicability across diverse photophysical playgrounds beyond TADF. The resulting simulations quantitatively reproduce key experimental observables, including PLQY and prompt/delayed fluorescence lifetimes. Beyond its predictive power, the present results establish clear criteria for identifying when higher-lying excited states influence the excited-state decay and when simplified models remain adequate. This framework enables rational selection of minimal kinetic models that balance physical insight with numerical robustness, with direct implications for the in silico design of high-performance organic emitters.

---

1 Introduction

Understanding the excited-state dynamics of organic chromophores has traditionally relied on a simplified three-state picture involving the ground state (S₀) and the lowest singlet (S₁) and triplet (T₁) excited states. This approximation is generally valid when S₁ and T₁ are well separated in energy from higher-lying excited states (S_n>1 and T_n>1). However, advances in quantum-chemical methods and ultrafast spectroscopies have increasingly demonstrated the protagonist role of higher-lying excited states in a broad variety of photophysical and photochemical scenarios. Representative examples include emission^1,2 and photochemical reactions³ originating from higher-lying excited states (i.e., anti-Kasha photochemistry), triplet formation through S₁ → T_n>1 intersystem crossing (ISC) processes,^4–7 and internal conversion (IC) mediated by higher-lying excited states that enable access to pathways inaccessible from S₁.^8,9 In the above scenarios, a conventional three-state model becomes inadequate to model the excited-state dynamics of organic chromophores, as the explicit inclusion of higher-lying states may qualitatively and quantitatively alter the photophysical outcomes. From an experimental perspective, unambiguously establishing their involvement remains challenging because it requires excitation-wavelength-dependent time-resolved spectroscopy combined with global kinetic analysis across multiple excited-state models and timescales.^9–11 Consequently, such tremendous efforts are rarely reported in the literature. For instance, employing coherence spectroscopy and wavepackets simulations in the context of the ultrafast excited state dynamics of dinuclear Pt(II) complexes, Chen et al.^12–14 provided an experimental and computational evidence of the spin-vibronic mechanism being at play and demonstrated the role of higher-lying excited state to accelerating ISC processes.

The above considerations are particularly relevant for thermally activated delayed fluorescence (TADF) materials,¹⁵ for which a three-state picture is often assumed. However, increasing evidence indicates that more complex excited-state kinetics beyond this picture are not uncommon in TADF emitters.^16–19 In these systems, delayed fluorescence arises from reverse intersystem crossing (rISC) that upconverts triplet excited states (T_n) to singlet excited states (S_n), rendering observables such as photoluminescence quantum yield (PLQY) and prompt/delayed fluorescence lifetimes highly sensitive to the underlying excited-state decay pathways.²⁰ More recently, multi-resonance TADF (MR-TADF) emitters²¹ have attracted considerable attention owing to their narrowband emission and efficient rISC enabled by short-range charge-transfer (SRCT) character and small singlet-triplet energy gaps.^22–24 As mentioned previously, increasing experimental and theoretical studies have suggested that higher-lying excited states may play a key mechanistic role in mediating ISC and rISC processes in these systems. Specifically, spin-vibronic coupling involving higher-lying triplet states has been shown to strongly enhance ISC and rISC processes especially when direct S₁ → T₁ or T₁ → S₁ coupling is symmetry-forbidden or intrinsically weak.^25–28 Computational studies by Olivier et al.^16,17 and Marian et al.¹⁸ have further demonstrated that “dark” triplet states with distinct orbital characters can act as intermediate spin-mixing channels, facilitating efficient population transfer between S₁ and T₁.

Despite these mechanistic insights, most quantitative kinetic analyses of TADF photophysics continue to rely predominantly on simplified three-state models. Experimentally, the transient photoluminescence measurements of TADF emitters are typically fitted by a biexponential decay profile, where the fast component is assigned to the prompt fluorescence and the slow one to the delayed fluorescence. Early kinetic models, such as those by Adachi et al.²⁹ and Dias et al.²⁰ estimated rISC rates under simplifying assumptions, e.g., neglecting non-radiative decay or phosphorescence, whereas later approaches by Monkman et al.³⁰ and Tsuchiya et al.³¹ have progressively relaxed these constraints and enabled direct extraction of all relevant kinetic parameters. Nevertheless, these models generally assume that ISC and rISC occur solely between S₁ and T₁. From a mathematical perspective, the limited number of experimentally measurable quantities (i.e., PLQY and fluorescence lifetimes) is insufficient to uniquely determine all unknowns in the excited-state decay kinetics. Consequently, the roles of higher-lying states cannot be unambiguously proved from experiment alone, making computational simulations indispensable for achieving a complete and quantitative description of the excited-state dynamics. Unlike the above-mentioned kinetic models that extract rate constants by fitting photoluminescence transient data, several theoretical studies have instead attempted to first compute rate constants and then use kinetic models to predict PLQY and fluorescence lifetimes.³² Shizu et al.^33,34 have made significant contributions to the quantitative modeling of excited-state kinetics. They applied their method to molecules like DABNA-1, BNOO, BNSS, and BNSeSe. As with any theoretical framework, their methodology involves certain approximations that define its range of applicability. For example, they neglect Herzberg–Teller (HT) vibronic coupling effects in the rate constant calculations, leading to an underestimation of (r)ISC rate constants. This deficiency is evident from the calculated rISC T₁ → S₁ rate constants in DABNA-1 (2.5 s⁻¹ versus an experimental value of 9.9 × 10³ s⁻¹) and BNOO (8.5 s⁻¹ versus 4.3 × 10⁴ s⁻¹).

As a result, a fundamental gap remains in our understanding of the excited-state dynamics of organic chromophores, including TADF emitters: although higher-lying excited states are frequently invoked qualitatively, their explicit and quantitative impact on PLQY and excited-state lifetimes has not been systematically assessed within a fully ab initio kinetic modelling framework.

In this study, we address this gap by developing multistate kinetic models explicitly including higher-lying singlet and triplet states (S₂ and T₂) and incorporating HT vibronic coupling in the rate constant calculations. Our framework provides the community with a quantitative basis for determining when explicit inclusion of higher-lying excited states becomes essential and when a simplified kinetic model (i.e., three-state picture) remains adequate. We first present a computational protocol that enables the calculation of all radiative and non-radiative rate constants involved in the interconversion between the electronically excited states. We then introduce KinLuv, a Python-based kinetic modeling program that computes prompt and delayed fluorescence lifetimes and PLQY from the calculated rate constants. To validate the proposed framework, we adopt a hierarchical strategy. DOBNA and DiKTa serve as prototypical MR-TADF systems, while DABNA-1 (a DOBNA derivative) and DQAO (a DiKTa derivative) assess the applicability of our approach to related derivatives. To further evaluate the generality of our methods, we additionally consider two non-TADF organic chromophores, namely dibenzothiophene (DBT) and bis-benzo-fused pyrano[3,2-b]pyran-2,6-dione (PPDs-1). DBT³⁵ and PPDs-1 ³⁶ were chosen as representative organic luminophores whose photophysical properties are strongly influenced by higher-lying excited states. In DBT, triplet formation through S₁ → T_n>1 ISC processes has been proved,^37,38 while preliminary investigations in our group show IC mediated by higher-lying excited states in PPDs-1. By systematically comparing multistate kinetic models of increasing complexity, this unified framework enables us to identify when higher-lying excited states play a significant role in their excited-state dynamics, to delineate the limitations of simplified kinetic models, and to quantify the trade-off between mechanistic completeness (inclusion of all physically accessible excited states and pathways) and quantitative necessity (inclusion required to reproduce experimentally observable quantities) with a fully ab initio kinetic modeling framework. This framework offers a principled strategy for constructing minimal yet robust kinetic models, thereby guiding the in silico design of high-performance organic emitters.

2 Results and discussion

2.1 Excited state kinetic models

If a molecule in its excited state is regarded as a distinct chemical species from that in its ground state, the interconversion between its involved electronic states can be formulated as a chemical reaction chain and analyzed with conventional kinetic methods. Within this framework, we systematically examined a hierarchy of excited state kinetic models with increasing complexity to evaluate the mechanistic relevance of higher-lying excited states in the quantitative calculation of PLQY and excited state lifetimes for both TADF and non-TADF organic chromophores.

Fig. 1 schematically summarizes the four kinetic models considered in this study, which differ in the number of electronic states and the possible interconversion processes between them. The two-state model considers the ground state (S₀) and the first singlet excited state (S₁) (Fig. 1a) and is therefore applicable only to model prompt fluorescence in the simplest dyes. The three-state model shown in Fig. 1b additionally includes the first triplet excited state (T₁), allowing for an approximate description of TADF. However, as mentioned above, higher-lying excited states cannot be ruled out in many molecular systems,¹⁹ necessitating more complex excited-state kinetic model. This is illustrated here by the four- and five-state models, which include a higher triplet excited state (T₂) and concomitantly a higher singlet excited state (S₂), respectively. To the best of our knowledge, the latter two models have not yet been widely applied to quantitatively study the excited-state dynamics of functional chromophores. Here, they are deliberately employed to assess the impact of higher-lying excited states on PLQY and prompt/delayed fluorescence lifetimes, across a diverse set of molecular playgrounds, including TADF emitters.

Fig. 1 Schematic diagram of the examined excited state kinetic models for the investigated chromophores: (a) two-state (S₀, S₁), (b) three-state (S₀, S₁, T₁), (c) four-state (S₀, S₁, T₁, T₂), (d) five-state (S₀, S₁, S₂, T₁, T₂) model. The considered interconversion processes between the involved electronic states are highlighted for each model: absorption (ABS), fluorescence (FL), phosphorescence (PH), internal conversion (IC), reverse internal conversion (rIC), intersystem crossing (ISC), and reverse intersystem crossing (rISC).

As shown in Fig. 1, all possible interconversion processes between adjacent electronic states, including absorption, fluorescence (both Kasha and anti-Kasha), phosphorescence, (r)IC and (r)ISC, are explicitly considered. By applying the law of mass action and the rate equation to each pathway, a set of ordinary differential equations (ODEs) is obtained in a general form (eqn (1)):

where R_i is the rate-constant matrix encompassing all possible pathways in the i-state model, N_i(t) is the time-dependent population vector of i states. Detailed ODEs for each model are provided in Section 2–3 of the SI. Once the rate constants are determined, the ODEs can be solved analytically or numerically to yield the population evolution of ground and excited states, from which the prompt and delayed fluorescence lifetime can be extracted.

While the time-dependent solutions of the above ODEs allow for simulating excited-state lifetime, PLQY is conventionally measured by steady-state conditions. In this limit, the system approaches a dynamic equilibrium in which all states populations become time independent. Applying this steady state approximation (SSA), the ODEs are reduced to ordinary algebraic equations (eqn (2)), from which the PLQY is directly derived from its definition (see Section 4 of the SI).

All relevant rate constants were calculated using our computational protocols (see theoretical formalisms in Section 1 of the SI), and the corresponding ODEs were automatically solved using our Python-based code package, KinLuv. The resulting PLQY and prompt/delayed fluorescence lifetimes of the studied emitters (Fig. 2) were subsequently compared with experimental data (Table 1). For the studied systems, anti-Kasha emission from the higher-lying singlet state S₂ is not considered, as no experimental evidence for such anomalous emission has been reported. In addition, rISC processes are not considered for DBT and PPDs-1, as delayed fluorescence has not been experimentally observed for these non-TADF compounds. This workflow enables a quantitative assessment of the validity and predictive capability of fully ab initio multistate kinetics models for functional organic chromophores, with particular emphasis on TADF systems, where both PLQY and prompt/delayed fluorescence lifetimes are of primary interest.

Fig. 2 Molecular structure of DOBNA (a), DiKTa (b), DABNA-1 (c), DQAO (d), DBT (e) and PPDs-1 (f).

Table 1 Experimental photophysical properties of selected molecules from previous studies

Molecule	ΦF (%)	τPF (ns)^g	τDF (µs)^h	ΔEST (S1-T1)^i
a Dispersed in PMMA (1 wt%) at 300 K.b In degassed dilute toluene at 300 K.c In mCBP films (1 wt%) at 300 K.d In mCP films (8 wt%) at 298 K.e In deaerated diluted cyclohexane at 298 K.f In dichloromethane (DCM) at 298 K.g Prompt fluorescence lifetime.h Delayed fluorescence lifetime.i Energy gap between S1 and T1 measured from emission spectra at 77 K.
DOBNA^a	57	11.6	66.1	0.18
DiKTa^b	26	5.1	23	0.18
DABNA-1^c	88	8.8	93.7	0.2
DQAO^d	59.3	10.09	110.58	0.19
DBT^e	9	3	—	—
PPDs-1^f	0.9	0.066	—	—

---

2.2 Excited-state energetics and decay rate constants

Six functional chromophores were deliberately selected as representative systems to examine the putative role of higher-lying excited states (Fig. 2). The selection follows a hierarchical design. We first focus on two prototypical MR-TADF emitters, DOBNA (Fig. 2a) and DiKTa (Fig. 2b), which serve as stringent benchmark systems for in-depth analysis. DOBNA features a central boron atom as the electron acceptor and oxygen atoms as electron donors, while DiKTa uses carbonyl groups as electron acceptors instead of boron. To assess the generality and robustness of the proposed methodology, DABNA-1 (Fig. 2c) and DQAO (Fig. 2d) were subsequently selected as structurally related derivatives of these prototype MR-TADF emitters. These systems allow us to examine whether the conclusions drawn for DOBNA and DiKTa persist upon molecular modification, and to further evaluate the transferability of our approach. Finally, DBT (Fig. 2e) and PPDs-1 (Fig. 2f) were included as representative non-TADF organic chromophores, as previously mentioned to be inadequately described by simplified three-state models. Their incorporation enables the assessment of our approach beyond TADF and provides a critical control for assessing the role of higher-lying excited states. Table 1 summarizes the previously reported experimental results for all the investigated compounds.^{22,24,35,36,39,40}

The optimized ground- and excited-state geometries of all systems are provided in Section 6 of SI: DOBNA and DiKTa (Fig. S3), DABNA-1 and DQAO (Fig. S4), and DBT and PPDs-1 (Fig. S5). In all cases, the structural changes associated with excitations are minimal (Root-Mean-Square Deviation well below 0.2 Å), thereby satisfying one of the prerequisites to apply the harmonic approximation and Fermi's golden rule (FGR) to calculate the excited state decay rate constants. The natural transition orbitals (NTOs) calculated with TD(A)-CAM-B3LYP involved in the excited states of selected compounds, i.e., DOBNA and DiKTa, are shown in Section 5 of the SI. These NTOs were compared with those obtained with SCS-ADC(2). As shown in Fig. S1 for DOBNA, the NTOs calculated using the CAM-B3LYP functional closely resemble those obtained from SCS-ADC(2). Similarly, for DiKTa (Fig. S2), the CAM-B3LYP NTOs show good agreement with those from SCS-ADC(2) for most of the excited states, except for one hole NTO in the T₂ state.

While geometries, Hessian matrix, and associated properties at the TD(A)-DFT level have been reported to be sufficiently accurate for excited state decay rate constant calculations, excitation energies are often insufficiently accurate.^41,42 To obtain refined adiabatic excitation energies and more accurate estimations of adiabatic energy gaps, we performed single-point SCS-ADC(2) calculations at the TD(A)-DFT optimized geometries. The comparison of the adiabatic excitation energies for the excited states calculated at the TD(A)-CAM-B3LYP and at the SCS-ADC(2) levels is shown in Section 7 of the SI. For instance, for the MR-TADF emitters, the adiabatic TD(A)-DFT excitation energies differ from the SCS-ADC(2) results by 0.1–0.8 eV (Fig. S6). Specifically, TD(A)-DFT tends to yield larger excitation energies for singlet states (S_n) and smaller energies for triplet states (T_n) compared to SCS-ADC(2), resulting in larger singlet-triplet energy gaps at the TD(A)-DFT level. Notably, this discrepancy is particularly pronounced for the T₂ state of DOBNA and DABNA-1, highlighting the necessity of using a higher-level method such as SCS-ADC(2) to refine the excited-state energy gaps. In the case of the non-TADF systems DBT and PPDs-1, deviations between TD(A)-DFT and SCS-ADC(2) excitation energies remain below 0.5 eV. Table 2 lists the adiabatic excitation energies calculated with SCS-ADC(2) for all the compounds, obtained by subtracting the ground-state energy at its optimized geometry from the excited-state energies at their respective minima. For DOBNA, the calculated S₁–T₁ energy gap, i.e., ΔE_ST(S₁–T₁), is 0.18 eV, which shows excellent agreement with the experimental value (Table 1). DiKTa exhibits a slightly larger calculated ΔE_ST(S₁–T₁) of 0.23 eV, which is within 0.05 eV of the experimental value, and comparable accuracy is obtained for DABNA-1 and DQAO. These results highlight the capability of the SCS-ADC(2) method to accurately determine the ΔE_ST(S₁–T₁) values of MR-TADF emitters. In contrast, the calculated ΔE_ST(S₁–T₁) values for the non-TADF emitters are sufficiently large (>0.7 eV), indicating a slowed down rISC process and consistent with the experimental observed absence of delayed fluorescence in these two compounds. Moreover, the ΔE_ST(S₁–T₂) values for DiKTa, DQAO, DBT and PPDs-1, in contrast to those of DOBNA and DABNA-1, are sufficiently small (within ±0.25 eV) to qualitatively suggest that triplet formation via higher-lying triplet excited states may be energetically feasible in these systems. However, this qualitative analysis is based solely on energetic considerations. In addition, a definitive quantitative assessment requires the approach developed here.

Table 2 Adiabatic excitation energies (eV) of selected molecules from SCS-ADC(2) calculations

Molecule	S1	S2	T1	T2	ΔEST (S1–T1)	ΔEST (S1–T2)
DOBNA	3.55	4.31	3.37	3.99	0.18	−0.44
DiKTa	3.27	3.45	3.04	3.40	0.23	−0.13
DABNA-1	3.15	3.98	3.02	3.69	0.13	−0.54
DQAO	3.26	4.02	3.04	3.41	0.23	−0.14
DBT	4.11	4.76	3.42	3.87	0.69	0.25
PPDs-1	3.66	3.97	2.51	3.68	1.15	−0.02

---

Thus, in the following section we focus on a detailed analysis of the calculated rate constants for DOBNA and DiKTa, which serve as representative prototype examples of MR-TADF emitters with expected different behaviors. For completeness, the corresponding data for the other two MR-TADF derivatives, i.e., DABNA-1 and DQAO, as well as for the two non-TADF emitters, are provided in Section 9 of SI. Fig. 3a and b present the calculated rate constants of (r)ISC between singlet states (S_n) and triplet states (T_n) for DOBNA (red) and DiKTa (green), respectively. The calculated ISC rate constants for the individual sublevels are provided in Tables S1–S2 of the SI. As seen in Fig. 3a and b, DiKTa exhibits larger (r)ISC rate constants than DOBNA for most transitions, with the exception of ISC T₁ → S₀. In particular, the larger S₁ → T₁ ISC rate constant calculated for DiKTa than for DOBNA cannot be explained by energetic considerations alone, as DOBNA has a slightly smaller ΔE_ST(S₁–T₁) value of 0.18 eV compared to 0.23 eV for DiKTa (Table 2). According to Marcus theory for ISC rate constants, two main factors account for this behavior: differences in spin–orbit coupling matrix elements (SOCMEs) or variations in the Marcus regime governing the ISC process (i.e., normal or inverted regime). As shown in Tables S4–S5, both compounds feature a very small reorganization energy (ca. 0.02 eV) compared to their singlet-triplet gaps, indicating that the S₁ → T₁ ISC occurs within the same Marcus regime for both compounds. In contrast, the calculated SOCMEs at the reference FC geometry (Tables S6–S7) indicate that DiKTa exhibits larger S₁ → T₁ SOCMEs than DOBNA. These results identify enhanced SOCMEs as the primary origin of the faster S₁ → T₁ ISC rate in DiKTa. Still, the calculated ISC rate constants for both DiKTa and DOBNA are exceptionally large, roughly ranging from 10⁵ s⁻¹ to 10¹¹ s⁻¹. Such high values cannot be rationalized solely based on calculated SOCMEs at FC geometry, nor by invoking simple empirical rules such as the heavy-atom effect⁴³ or El-Sayed rule.⁴⁴ The latter rule predicts faster ISC processes between states of different orbital characters (e.g., nπ* → ππ*). However, analysis of the calculated NTOs (Fig. S1–S2) reveals that all relevant transitions are of ππ* character, with negligible nπ* contributions. Instead, the large (r)ISC rate constants stem from the vibronic coupling, as the HT effect is explicitly considered in the rate constant calculations. As shown in Tables S8–S9, the calculated (r)ISC rate constants are dominated by HT contributions for nearly all transitions. Unlike the FC approximation, the HT effect enables electronic transitions to couple with vibrational motions via a breakdown of the Condon approximation.⁴⁵ The vibrational motions of nuclei effectively increase the SOCMEs, thereby facilitating rapid ISC processes.

Fig. 3 Calculated rate constants of (a) intersystem crossing (ISC), (b) reverse intersystem crossing (rISC), (c) internal conversion (IC), (d) reverse internal conversion (rIC) and (e) emission for molecule DOBNA (red) and DiKTa (green). Calculations at SCS-ADC(2)/def2-TZVP//TDA-CAM-B3LYP/6-311G(d,p) level with the adiabatic Hessian, Franck–Condon/Herzberg–Teller (FC/HT) vibronic model and a 10 cm⁻¹ Lorentzian broadening.

Fig. S7–S8 present the Huang-Rhys factors and reorganization energies per vibrational normal mode for key ISC processes of DOBNA and DiKTa, revealing a strong mode dependence. For DOBNA, in-plane vibrations contribute the most to r(ISC) processes, as they have both the largest Huang-Rhys factors and the largest reorganization energies. Conversely, out-of-plane modes exhibit the largest Huang-Rhys factors for DiKTa, but in-plane modes still contribute the most to the reorganization energies. Overall, these results highlight that out-of-plane vibrational modes have more relevance for DiKTa than for DOBNA. To further elucidate the origins of the HT effects, Fig. S9 and S10 show the SOCMEs calculated at the reference FC geometry and at displaced geometries along specific vibrational modes corresponding to the HT regime for the key (r)ISC channels of DOBNA and DiKTa, respectively. Notably, the SOCMEs in the HT regime are significantly enhanced compared to those in the FC regime particularly for the S₁ → T₁ (Fig. S9a and S10a) ISC and the T₁ → S₁ (Fig. S9e and S10e) rISC. Structurally, all atoms in DOBNA lie in a strictly planar configuration, leading to predominantly in-plane vibrational modes. In contrast, DiKTa exhibits a slight torsion between its molecular fragments, allowing out-of-plane vibrations that more effectively break molecular symmetry and enhance vibronic coupling.⁴⁶ This structural distinction likely contributes to the stronger HT effect observed in DiKTa (Fig. S9 and S10). This analysis demonstrates that vibronic couplings effectively enhance SOCMEs via the HT effect. All in all, vibrational contributions to HT-driven (r)ISC exhibit a strong dependence on molecular structure. For planar systems (e.g., DOBNA), in-plane vibrational modes dominate (Fig. S7 and S9), whereas in slightly twisted systems (e.g., DiKTa), out-of-plane modes play a more significant role, as they more effectively break symmetry and enhance vibronic coupling (Fig. S8 and S10). Notably, the enhancement of SOCMEs arises from the collective contributions of multiple vibrational modes rather than from a single dominant mode.

Fig. 3c and d show the rate constants for IC and rIC; respectively, for DOBNA and DiKTa. The S₁ → S₀ IC rate constants are on the order of 10⁷ s⁻¹. Notably, in DOBNA (Fig. 3d), the calculated rIC rate constant for T₁ → T₂ is comparable in magnitude to that of the T₂ → T₁ IC process and is significantly larger than those of T₂ → S₁ and T₁ → S₁ processes, indicating an efficient and nearly reversible population exchange between the two triplet states prior to population transfer back to the singlet manifold. Fig. 3e collects the values of the calculated radiative rate constants, which are ca. 10⁸ s⁻¹ for fluorescence and ca. 1 s⁻¹ for phosphorescence, aligned to typical values for these compounds.

The calculated excited state decay rate constants for MR-TADF derivatives DABNA-1 and DQAO, as well as for the non-TADF systems (DBT, PPDs-1), are briefly summarized below. Calculated values are provided in Section 9 of the SI. DABNA-1, being closely related to DOBNA shows similar trends in the computed rate constants. Both compounds generally possess smaller (r)ISC rate constants than DQAO, which is closely related to DiKTa and, in turn, also exhibits analogous rate constant trends (Tables S10 and S11). Importantly, by explicitly including HT contributions in the rate constant calculations, the rISC T₁ → S₁ rate constant of DABNA-1 (2.41 × 10⁴ s⁻¹) is obtained with a more reasonable accuracy with respect to the experimental one (9.9 × 10³ s⁻¹) than the one reported by Shizu et al.³³ (2.5 s⁻¹). This comparison highlights the indispensable role of vibronic coupling in accurately describing (r)ISC processes. Interestingly, the calculated rIC rate constant for T₁ → T₂ in DABNA-1 (Table S10) is comparable in magnitude to the T₂ → T₁ IC rate constant. This trend is also found in DOBNA as previously mentioned. These rate constant behaviors observed for the prototype emitters are consistently reproduced in their derivatives, indicating the generality of the underlying trends. In addition, the magnitudes of the ISC rate constants are also substantial for non-TADF compounds DBT and PPDs-1, where the large values arise primarily from HT contributions (Table S15). Both of them show a larger ISC S₁ → T₂ rate constant than that of S₁ → T₁, indicating the non-negligible role of T₂ in mediating ISC. At the same time, PPDs-1 shows a comparably large S₁ → S₂ rIC rate constant compared to S₂ → S₁ IC, indicating a significant population exchange between these two singlet states.

2.3 Excited-state kinetics simulation of MR-TADF

The following discussion focuses on excited-state decay kinetics and the quantitative determination of PLQY and lifetimes. The discussion first focuses on DOBNA and DiKTa, as prototype MR-TADF emitters, and then on their structurally related derivatives (DABNA-1 and DQAO). The remaining systems, PPDs-1 and DBT, where higher-lying excited states are also relevant, are discussed thereafter. The complete excitation and decay profiles directly generated by KinLuv for all compounds are provided in Section 12 of SI. Fig. 4 and 5 collect the simulated decay kinetics of DOBNA and DiKTa using kinetic models of increasing complexity, respectively. A logarithmic time scale is employed to simultaneously capture both prompt and delayed emission regimes within a single plot. Additionally, KinLuv can directly output the PLQY as well as the prompt and delayed fluorescence lifetimes. For the two- and three-state models, the decay of S₁ population can be described by analytical expressions, whereas numerical solutions are employed for the four- and five-state models to ensure computational efficiency. The calculated PLQY, prompt and delayed lifetimes for the four investigated MR-TADF emitters are summarized in Fig. 6.

Fig. 4 Simulated population decay kinetics of DOBNA over a 1 ms timescale based on (a) two-state (S₀, S₁), (b) three-state (S₀, S₁, T₁), (c) four-state (S₀, S₁, T₁, T₂), and (d) five-state (S₀, S₁, S₂, T₁, T₂) kinetic models. The number of absorbed photons is set to 1, with the excitation pulse width of 10 ps. Both time (x-axis) and population (y-axis) are plotted on logarithmic scales.

Fig. 5 Simulated population decay kinetics of DiKTa over a 1 ms timescale based on (a) two-state (S₀, S₁), (b) three-state (S₀, S₁, T₁), (c) four-state (S₀, S₁, T₁, T₂), and (d) five-state (S₀, S₁, S₂, T₁, T₂) kinetic models. The number of absorbed photons is set to 1, with the excitation pulse width of 10 ps. Both time (x-axis) and population (y-axis) are plotted on logarithmic scales.

Fig. 6 Calculated (a) photoluminescence quantum yield (PLQY), (b) prompt fluorescence lifetime, and (c) delayed fluorescence lifetime for DOBNA (red), DABNA-1 (blue), DiKTa (green), and DQAO (purple) compared with experimental results.

DOBNA is first discussed as a representative example (Fig. 4). In the two-state model (Fig. 4a), after initial population of S₁ (green), radiative and nonradiative decay to S₀ (red) occur in the nanosecond regime, giving rise exclusively to prompt fluorescence. Incorporation of the T₁ state in the three-state model (Fig. 4b) introduces an ISC pathway from S₁ to T₁ (blue), followed by rISC back to S₁. Consequently, two distinct regimes showing a marked decrease in S₁ population are observed in Fig. 4b. The first regime corresponds to prompt fluorescence (as in the two-state model) while the second one (ranging between ca. 10² ns and ca. 10² µs) stems from the delayed fluorescence via the T₁ → S₁ → S₀ pathway. However, only a minute fraction of T₁ population is accumulated because the S₁ → T₁ ISC rate constant (1.2 × 10⁶ s⁻¹) is much smaller than the competing fluorescence (6.4 × 10⁷ s⁻¹) and IC (2.1 × 10⁷ s⁻¹) rate constants. Accordingly, the delayed decay of the S₁ population originates from rISC of the T₁ population and is characterized by a small onset amplitude of ca. 10^–6 (Fig. 4b). In the four-state model, inclusion of the higher-lying triplet state T₂ leads to nearly overlapping T₁ (blue) and T₂ (purple) population curves (Fig. 4c and d), indicating synchronized population dynamics. This behavior arises from the comparable rIC rate constant for T₁ → T₂ (3.9 × 10¹⁰ s⁻¹) and the T₂ → T₁ IC rate constant (9.8 × 10¹⁰ s⁻¹), resulting in rapid thermal equilibration between these two states. Although T₂ state introduces an alternative delayed fluorescence pathway (T₁ → T₂ → S₁ → S₀), the limited population accumulated in T₁ restricts the overall impact of this pathway, increasing the onset amplitude of delayed S₁ population only modestly, from ca. 10^–6 (in the three-state model) to ca. 10^–4. In the five-state model, the initial population of the higher-lying singlet state S₂ is followed by ultrafast IC to S₁ (8.6 × 10¹⁰ s⁻¹), which leads to negligible effect on the global excited-state dynamics of the emitter.

Despite the mechanistic differences between the three-, four-and five-state models, the calculated PLQY remains nearly constant at ca. 74% across all models (see red bars in Fig. 6). The calculated prompt fluorescence lifetime reproduces the experimental value of 11.6 ns in the three-state model and shows only minor changes in the four- and five-state model. These results indicate that for DOBNA while the inclusion of higher-lying excited states is mechanistically important to provide a more complete description of the TADF pathways, it is not essential for accurately predicting experimentally observables such as PLQY, prompt and delayed fluorescence lifetimes. This can be attributed to the relatively slower timescales of (r)ISC processes compared to fluorescence and IC, which limit population accumulation in higher-lying triplets and thereby minimize their impact on the overall decay dynamics.

As a derivative of DOBNA, DABNA-1 exhibits excited-state dynamics that are very similar to those of DOBNA (Fig. S11). The calculated PLQY (82% versus experimental 88%) and prompt fluorescence lifetime (8.4 ns versus 8.8 ns) remain essentially unchanged across the three-, four-, and five-state models (see blue bars in Fig. 6). By contrast, the delayed fluorescence lifetime decreases markedly from 41 µs in the three-state model to 2.3 µs in the four and five state model. Consistent with the behavior observed for DOBNA, higher-lying states may impact the TADF pathways mechanistically, yet the slow (r)ISC rates ensure a negligible effect on PLQY, prompt and delayed fluorescence. The consistency between DOBNA and its derivative DABNA-1 demonstrates that the conclusions drawn from the prototype system extend beyond a single molecule, thereby underscoring the robustness of our approach. Taken together, these results demonstrate that the adequacy of simplified three-state kinetic models is governed by both the magnitude of (r)ISC rate constants relative to fluorescence and IC and whether higher-lying excited states provide kinetically competitive rISC pathways. In systems such as DOBNA and DABNA-1, where (r)ISC processes are slow and the population of higher-lying states remains kinetically limited, their contribution to rISC is secondary, and three-state models remain sufficient to reproduce experimentally observable quantities, including PLQY and fluorescence lifetimes.

DiKTa (Fig. 5) exhibits significantly faster (r)ISC rate constants than DOBNA (Fig. 3a and b), leading to markedly distinct excited-state decay kinetics. In the three-state model (Fig. 5b), the T₁ population becomes significantly larger than that of DOBNA (Fig. 4b) because of the much larger S₁ to T₁ ISC rate constant (1.1 × 10¹⁰ s⁻¹) than that of DOBNA (1.2 × 10⁶ s⁻¹). Moreover, fast S₁ → T₁ ISC suppresses direct decay from S₁ to S₀, causing a pronounced lag in S₀ population growth. As S₁ is gradually repopulated through rISC from T₁, the S₀ population increases again, producing a characteristic sigmoidal profile which was not observed for DOBNA. When moving to the four-state model (Fig. 5c), the inclusion of a higher-lying triplet state T₂ does not introduce an efficient delayed fluorescence pathway analogous to T₁ → T₂ → S₁ → S₀ observed for DOBNA. This is because IC rate constant for T₂ → T₁ (1.4 × 10¹⁴ s⁻¹) is substantially faster than the rIC T₁ → T₂ rate constant (3.3 × 10¹⁰ s⁻¹). Instead, T₂ primarily provides a competing non-radiative pathway (S₁ → T₂ → T₁ → S₀) that effectively reduces the PLQY. In the five-state model, the higher-lying state S₂ introduces an additional delayed fluorescence pathway (T₁ → S₂ → S₁ → S₀). However, this turns out to be a minor contribution that slightly increases the PLQY value, as the rISC rate constant for T₁ → S₂ (2.9 × 10⁵ s⁻¹) is smaller than that for T₁ → S₁ rISC rate constant (1.4 × 10⁶ s⁻¹).

Unlike DOBNA, the higher-lying states in DiKTa play a quantitatively important role in the excited-state kinetics. The calculated fluorescence PLQY decreases progressively from 83.1% in the two-state model to 28.6% in the three-state model and 25.4% and 26.5% in the four- and five-state model, the latter two models yielding PLQY values in excellent agreement with the experimental value of 26% (see the green bars in Fig. 6). Mechanistically, the computed decrease in the PLQY value when going from the three-state to the four-state and five-state models originates from the introduction of T₂ as an additional non-radiative decay pathway, whereas S₂ partially compensates for this loss by providing an extra delayed fluorescence pathway. Overall, while the predominant delayed fluorescence pathway in DiKTa remains that of the three-state model, the inclusion of the higher-lying excited states improves the quantitative agreement with respect to the experiment. In addition, the predicted delayed fluorescence lifetime increases from 20.57 µs in the three-state model to 21.50 µs and 21.02 µs in the four and five state models, the latter two in excellent agreement with the experimental value of 23 µs. All in all, these results show that in DiKTa the inclusion of T₂ is not only mechanistically important but also essential for quantitatively reproducing the experimental PLQY, owing to fast (r)ISC and substantial triplet state populations. In contrast, S₂ primarily improves the mechanistic description of TADF pathways but it is not required to accurately predict the experimentally observable.

As a derivative of DiKTa, DQAO exhibits similar excited-state dynamics (Fig. S12). The calculated fluorescence PLQY decreases progressively from 77.5% in the two-state model to 66.1% in the five-state model. The comparison of the PLQY values for the three-state (49.3%) and the four-state model (60.7%) reveals that only after the inclusion of T₂ in the kinetic models can we recover the experimental value of 59% (see the purple bars in Fig. 6). In addition, the delayed fluorescence lifetime decreases from 38.75 µs in the three-state model to 22.92 µs in the four-state model, and to 13.67 µs in the five-state model. Consistent with the results obtained for DiKTa, the inclusion of T₂ in DQAO is therefore essential for quantitatively reproducing the experimental PLQY, reflecting the fast (r)ISC processes and the dominant triplet state populations. Unlike in DiKTa, T₂ in DQAO additionally mediates a minor alternative delayed fluorescence pathway (T₁ → T₂ → S₁ → S₀), as reflected in the changes in delayed fluorescence lifetime, rather than acting solely as a competing non-radiative decay pathway. The inclusion of S₂ further refines the mechanistic description of the TADF process but has only a minor impact on the predicted experimentally observables. Overall, these results demonstrate that the adequacy of four-state kinetic models is likewise determined by the magnitude of (r)ISC rate constants relative to fluorescence and IC and by whether higher-lying triplet states provide kinetically competitive rISC pathways. In systems such as DiKTa and DQAO, where fast (r)ISC leads to substantial triplet state populations, the explicit inclusion of T₂ is required to quantitatively reproduce experimentally observable quantities.

2.4 Excited-state kinetics simulation beyond TADF emitters

As discussed above, higher-lying excited states may not only modify delayed fluorescence pathways but can also introduce additional nonradiative decay channels. Consequently, for certain non-TADF organic emitters, a three-state kinetic model is likewise insufficient. This behavior is illustrated below using DBT (see its excited state population dynamics in Fig. S13) and PPDs-1 (see its excited state population dynamics in Fig. S14).

For DBT, the three-state model already predicts a substantial population of T₁ (Fig. S13b) from fast ISC S₁ → T₁ (2.13 × 10⁸ s⁻¹), giving rise to a nonradiative decay pathway S₁ → T₁ → S₀. Consequently, the S₀ population exhibits a sigmoidal temporal profile, similar to the one observed by DiKTa. When moving to the four-state model, the inclusion of T₂ introduces an additional major non-radiative pathway (S₁ → T₂ → T₁ → S₀), because the ISC rate constant for S₁ → T₂ is also very large (2.65 × 10⁸ s⁻¹). In the five-state model, the higher-lying state S₂ provides another nonradiative pathway (S₁ → S₂ → S₀). However, its impact on the overall kinetics is of very minor relevance compared to the dominant ISC-mediated processes. As a result, the calculated PLQY progressively decreases from 73.3% in the two-state model to 17.1% in the three-state model and further to 8.8% in the four-state model and 8.5% in the five-state model (Fig. 7a). The PLQY of the four-state model is thus in good agreement with the experimental value of 9%. This highlights the necessity of including a higher-lying triplet excited state (T₂) to attain quantitatively accurate PLQY estimates. In parallel, the calculated prompt fluorescence lifetime decreases from 15.4 ns in the two-state model to 1.8 ns in the four-state model and 2.6 ns in the five-states model, closely matching the experimental value of ca. 3 ns. All in all, the higher-lying T₂ plays a crucial role by introducing competing nonradiative decay pathways (S₁ → T₂ → T₁ → S₀), which are essential for quantitatively reproducing both PLQY and prompt fluorescence lifetimes in systems dominated by fast ISC processes involving higher-lying triplets. Accordingly, the four-state model, explicit inclusion of T₂ represents the minimal extension needed to achieve quantitative agreement for non-TADF organic emitters in which ISC processes involving higher-lying triplets are significant. We note, however, that the present analysis does not exclude the possible involvement of higher-lying triplet states beyond T₂ (e.g., T₃) in other molecular systems. In such cases, extensions of the present framework to include additional excited states may be required to achieve quantitative accuracy, which is beyond the scope of the present work.

Fig. 7 Calculated photoluminescence quantum yield (PLQY) and prompt fluorescence lifetime for (a) DBT (light green) and (b) PPDs-1 (orange), compared with experimental results.

For PPDs-1, the three-state model (Fig. S14b) predicts a low T₁ population due to the relatively slow S₁ → T₁ ISC (4.61 × 10³ s⁻¹) relative to S₁ → S₀ IC (1.67 × 10⁷ s⁻¹) and fluorescence (3.74 × 10⁸ s⁻¹). In the four-state model, a faster ISC from S₁ → T₂ (4.07 × 10⁷ s⁻¹) increases the T₁ population and introduces an additional nonradiative decay pathway (S₁ → T₂ → T₁ → S₀). In the five-state model, the presence of the higher-lying singlet state S₂ opens up a highly efficient nonradiative decay channel, S₁ → S₂ → S₀. This pathway is important because the rIC S₁ → S₂ rate constant (3.15 × 10⁹ s⁻¹) exceeds that of the direct IC S₁ → S₀ (1.67 × 10⁷ s⁻¹), allowing a fraction of the S₁ population to access S₂ and subsequently undergo rapid IC (3.86 × 10¹⁰ s⁻¹) to S₀. As a result, the calculated PLQY decreases progressively from 95.7% in the two-state model to 77% in the five-state model, while the calculated prompt fluorescence lifetime decreases from 2.56 ns in the two-state model to 2.09 ns in the five-state model. This clear trend indicates that T₂ and S₂ do indeed introduce additional nonradiative pathways and therefore play an important role. However, the experimental PLQY is reported to be 0.9%. To investigate the discrepancies between the simulations and the experiments in greater depth, we performed further analyses. First, we excluded the possibility that nonradiative decay pathways involving the higher triplet state T₃ are operative, because both ISC (S₁ → T₃) (7.36 × 10⁴ s⁻¹) and ISC (S₂ → T₃) (9.87 × 10⁵ s⁻¹) are negligible compared to the fast (r)IC processes. Second, IC (S₂ → S₁) and rIC (S₁ → S₂) are extremely sensitive to the adiabatic energy gap ΔE(S₂–S₁). By subtracting 0.2 eV from the calculated ΔE(S₂–S₁) value at the SCS-ADC(2) level (we note that ±0.2 eV is the expected accuracy of SCS-ADC(2)), the calculated (r)IC rate constants significantly change by one to two orders of magnitude (Fig. S17). This leads to a reduction in the calculated PLQY to 4.0% and the lifetime to 0.14 ns (see fifth columns in Fig. 7b), bringing them substantially closer to the experimental values of 0.9% and 0.066 ns. Furthermore, we also calculated ΔE(S₂–S₁) using the CC2 method, which is of comparable accuracy to ADC(2) but has been reported to provide more reliable excitation energies for valence excitations in benchmark studies^47–49 (both S₁ and S₂ of PPDs-1 are of ππ* character, see Fig. S18). CC2 predicts a smaller value of 0.23 eV compared to that obtained with SCS-ADC(2) (0.31 eV). Correspondingly, the calculated rIC S₁ → S₂ rate constant increases by approximately one order of magnitude, from 3.15 × 10⁹ s⁻¹ (SCS-ADC(2)) to 1.27 × 10¹⁰ s⁻¹ (CC2), while the IC S₂ → S₁ rate constant increases more moderately, i.e., from 2.54 × 10¹² s⁻¹ (SCS-ADC(2)) to 4.37 × 10¹² s⁻¹ (CC2). Consequently, the calculated PLQY and prompt fluorescence lifetime further decrease to 68.3% and 1.85 ns, respectively (see fourth columns in Fig. 7b), in closer agreement with the experimental results than those obtained using the calculated ΔE(S₂–S₁) value at the SCS-ADC(2) level. Given the pronounced sensitivity of (r)IC rate constants to ΔE(S₂–S₁) (see Fig. S17) and the likely overestimation of valence excitation energies by ca. 0.2 eV at the SCS-ADC(2) level, we consider the inaccuracy in the calculated ΔE(S₂–S₁) value to be a major contributor to the remaining quantitative discrepancy with experiment (see also the discussion below in Section 2.5). These results demonstrate that higher-lying singlet states S₂ can play a significant role by introducing efficient competing nonradiative decay pathways (S₁ → S₂ → S₀), which are essential for quantitatively reproducing both PLQY and prompt fluorescence lifetimes in systems influenced by fast rIC processes. Accordingly, the five-state model, explicitly including both T₂ and S₂, represents the minimal kinetic description beyond the three-state description required to achieve quantitative agreement for non-TADF emitters dominated by fast ISC and rIC processes.

2.5 Accuracy and limitations of our approach

Including higher-lying excited states in the kinetic models improves the mechanistic completeness of the excited state dynamic description and can yield more reliable predictions of PLQY and excited-state lifetimes by explicitly accounting for all relevant photophysical pathways in the solution of the rate equations. Nevertheless, achieving quantitative agreement with experiment from ab initio calculations remains highly challenging, largely due to uncertainties in the calculated excited-state decay rate constants. For instance, as illustrated in Fig. S15, even small variations in the calculated energy gaps ΔE(S₁–S₀) can lead to orders-of-magnitude changes in the calculated S₁ → S₀ IC rate constants. Even when a very narrow Lorentzian broadening function is employed (HWHM∼0.001 eV), the physical vibrational overlap becomes negligible once the adiabatic energy lies beyond the inflection point. In such cases, the computed IC rate constant may be dominated by the tails of the broadening function, potentially leading to an overestimation. Consequently, extending the kinetic model to include additional excited states and pathways inevitably increases the number of rate constants that must be estimated, which can result in the accumulation and propagation of errors in the predicted observables.

For DOBNA and its derivative DABNA-1, the calculated delayed fluorescence lifetimes obtained from the four- and five-state kinetic models are underestimated by nearly one order of magnitude with respect to the experimental values. Even for the most rigid systems, errors of this size in the computed lifetimes are not uncommon. All in all, the deviations in the predicted PLQY and delayed fluorescence lifetime are to be unambiguously traced back to inaccuracies in some specific rate constant calculations. For DOBNA (Fig. S16a), gradually varying the value of the S₁ → S₀ IC rate constant over a range of 10⁷–10⁸ s⁻¹ while keeping all the other computed rate constants fixed leads to a pronounced change in the calculated PLQY, decreasing from approximately 80% to 40%, whereas the prompt and delayed fluorescence lifetimes remain essentially unaffected. This indicates that uncertainties in the calculated S₁ → S₀ IC rate constants are the primary source of the PLQY deviation. Furthermore, Fig. S16b demonstrates that variations in the T₁ → S₀ ISC rate constant strongly affect the delayed fluorescence lifetime, identifying this process as the dominant contributor to the inaccuracy in the predicted delayed lifetime for DOBNA. For its derivative DABNA-1 (Fig. S16c and S16d), the analysis indicates that uncertainties in the calculation of the S₁ → S₀ IC rate constant predominantly affect the predicted PLQY, whereas variations in the T₁ → S₁ ISC rate constant largely govern the accuracy of the delayed fluorescence lifetime.

For DiKTa and its derivative DQAO, the three-, four-, and five-state kinetic models all slightly underestimate the prompt fluorescence lifetime by nearly one order of magnitude. Among them, the four-state model yields the closest overall agreement with experiment when prompt and delayed fluorescence lifetimes and PLQY are considered simultaneously. In contrast to DOBNA and DABNA-1, the deviations observed for DiKTa and DQAO cannot be traced back to a single dominant rate constant but rather reflect the coupled and interdependent contributions of multiple excited-state decay processes. For the non-TADF emitters DBT, the calculated lifetimes and PLQY show deviations from the experimental data within an acceptable range. For PPDs-1, as discussed above, an overestimation of ΔE(S₂–S₁) is the major contributor to the remaining discrepancy with the experiment (Fig. S17). Additionally, it is possible that the (r)IC process between S₁ and S₂ proceeds through a conical intersection (CI), resulting in ultrafast dynamics that cannot be adequately described by Fermi's golden rule within the harmonic approximation. This scenario merits further exploration using e.g., nonadiabatic excited-state dynamics methods.

All in all, the above findings highlight a limitation of our fully ab initio kinetic modeling approach: the quantitative accuracy is ultimately determined by the accuracy in the state-of-the-art approaches to calculate excited state decay rate constants. First, the static approach we used to calculate rate constants relies on the harmonic approximation, which may introduce non-negligible errors in systems with significant anharmonicity. In such cases, more advanced approaches, such as Morse-type potential⁵⁰ and time-dependent density matrix renormalization group theory (TD-DMRG),⁵¹ may provide more reliable approaches to calculate rate constants. Second, nonadiabatic spin-vibronic effects remain challenging to incorporate within this static framework. These effects can be more naturally captured by ab initio excited-state dynamics methods (e.g., trajectory surface hopping^52,53). However, their high computational cost generally restricts their application to short time windows. In contrast, kinetic modeling enables the simulation of population evolution over extended (nanosecond to microsecond) timescales, thereby providing a direct connection to experimental photophysical observables of rigid chromophores. The complementarity of these two approaches represents a promising direction for future studies in computational photochemistry. Specifically, static rate kinetic models can be used to construct a global photophysical picture of the excited-state population evolution over long time scales and to predict observables such as PLQY and lifetimes. In regimes where the static approach becomes inadequate, such as in the presence of strong anharmonicity, near conical intersections (CI), or under pronounced spin-vibronic coupling, ab initio excited-state dynamics simulations can offer highly accurate mechanistic insight into short-time population evolution.^54,55 Overall, this multiscale simulation framework, integrating both long- and short-time dynamics, represents a highly promising paradigm for the future of computational photochemistry. Nevertheless, despite the inherent uncertainties in the calculated rate constants, our investigations unambiguously highlight that a three-state model may be insufficient for some MR-TADF and non-TADF organic emitters, underscoring the need for more comprehensive excited-state kinetic models.

3 Conclusion

In this study, we developed KinLuv, a Python-based kinetic modeling package which explicitly includes higher-lying excited states (T₂ and S₂) into the excited state kinetic analysis. All rate constants were computed using Fermi's golden rule, explicitly including the HT vibronic coupling effect. Our approach was applied to prototypical MR-TADF emitters (DOBNA, DiKTa) and their derivatives (DABNA-1, DQAO), as well as to non-TADF organic emitters (DBT, PPDs-1). Key findings are summarized as follows:

(a) In the investigated organic emitters, fast (r)ISC is primarily governed by HT vibronic coupling, as revealed by the relative contributions of HT and FC mechanisms to the calculated rate constants.

(b) Mechanistically, inclusion of T₂ state in the four-state model additionally mediates an additional delayed fluorescence pathway (T₁ → T₂ → S₁ → S₀) and a competing non-radiative decay pathway (S₁ → T₂ → T₁ → S₀). Inclusion of S₂ state in the five-state model enables an additional delayed fluorescence pathway (T₁ → S₂ → S₁ → S₀) and introduces further efficient competing nonradiative decay pathways (S₁ → S₂ → S₀).

(c) Quantitively, when (r)ISC processes are slow relative to fluorescence and IC, higher-lying excited states are not significantly populated and do not affect the overall decay dynamics; in this regime, a three-state description already reproduces the experimental observables (e.g., DOBNA and DABNA-1).

(d) When (r)ISC processes are fast leading to substantial triplet populations, and when the population of higher-lying triplet states become kinetically relevant; in this regime, explicit inclusion of at least T₂ is required to achieve quantitative agreement (e.g., DiKTa, DQAO and DBT).

(e) When rIC processes dominate the excited-state dynamics, higher-lying singlet states become essential nonradiative decay channels; in this regime, explicit inclusion of S₂ is required to reproduce the experimental observables (e.g., PPDs-1).

Overall, this unified framework clarifies that higher-lying excited states may play a mechanistically important role by introducing additional decay pathways and that their inclusion in the kinetic models may become necessary to attain quantitative accuracy in the calculated values of PLQY and fluorescence lifetimes. At the same time, increased mechanistic completeness does not inherently guarantee improved quantitative accuracy owing to the amplification of errors and uncertainties in the underlying calculation of ab initio rate constants. Further developments in FGR-based theories, including anharmonic effects, will be necessary to improve the quantitative accuracy of PLQY and fluorescence lifetimes. Regarding kinetic model selection, we suggest that higher-lying triplet states should be included explicitly when ISC is sufficiently fast to generate substantial populations in these states. Likewise, higher-lying singlet states become essential when rIC populates them more efficiently than nonradiative decay depletes the excited-state manifold. For more efficient model construction, simple descriptors may serve as preliminary indicators for assessing the importance of higher-lying states in future studies, thereby guiding model selection without requiring computationally intensive rate-constant calculations. Taken together, these insights establish a rational criterion for selecting minimal kinetic models that balance physical insight and numerical robustness, with direct implications for in silico design of high-performance organic emitters.

Author contributions

Y. H. performed the theoretical calculations, and D. E. supervised the project. All authors participated in the manuscript preparation and approved the final version for publication.

Conflicts of interest

The authors declare no competing interests.

Data availability

The source code, KinLuv, used in this study is publicly available on GitHub at: https://github.com/stevenuoa/KinLuv. Supplementary information: computational methods and theoretical background; ODEs and algebraic equations under the SSA; natural transition orbitals; optimized ground- and excited-state geometries; adiabatic excitation energies and rate constant calculations; error diagnosis of computational rate constants; and simulated excitation and decay kinetics generated by KinLuv. See DOI: https://doi.org/10.1039/d6sc01726f.

Acknowledgements

D. E. acknowledges FWO (project numbers G079122N, G022324N) and KU Leuven Internal Funds (NovelTADF, 3E250385). The quantum chemical calculations were performed on the VSC (Flemish Supercomputer Center), funded by the Research Foundation - Flanders (FWO) and the Flemish Government. We acknowledge the support of OpenAI's ChatGPT in polishing the language of this manuscript and assisting with code refinement.

References

1. T. Itoh, Chem. Rev., 2012, 112, 4541–4568.
2. K. Veys and D. Escudero, Acc. Chem. Res., 2022, 55, 2698–2707.
3. B. Pfund and O. S. Wenger, J. Am. Chem. Soc., 2025, 147, 26477–26485.
4. V. S. Reddy and S. Irle, J. Chem. Theory Comput., 2017, 13, 4944–4949.
5. A. W. Kohn, Z. Lin and T. Van Voorhis, J. Phys. Chem. C, 2019, 123, 15394–15402.
6. K. Shizu and H. Kaji, J. Phys. Chem. A, 2021, 125, 9000–9010.
7. S. Mai, P. Marquetand and L. González, J. Phys. Chem. Lett., 2016, 7, 1978–1983.
8. Z. Lin, A. W. Kohn and T. Van Voorhis, J. Phys. Chem. C, 2020, 124, 3925–3938.
9. Y. Sawada, A. Takada, K. Onda, K. Miyata, Y. Sasaki and N. Kimizuka, Phys. Chem. Chem. Phys., 2025, 27, 19677–19683.
10. S. Lochbrunner, T. Schultz, M. Schmitt, J. P. Shaffer, M. Z. Zgierski and A. Stolow, J. Chem. Phys., 2001, 114, 2519–2522.
11. C. M. Marian, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 187–203.
12. P. Kim, A. J. S. Valentine, S. Roy, A. W. Mills, A. Chakraborty, F. N. Castellano, X. Li and L. X. Chen, J. Phys. Chem. Lett., 2021, 12, 6794–6803.
13. P. Kim, A. J. S. Valentine, S. Roy, A. W. Mills, F. N. Castellano, X. Li and L. X. Chen, Faraday Discuss., 2022, 237, 259–273.
14. S. R. Rather, N. P. Weingartz, S. Kromer, F. N. Castellano and L. X. Chen, Nature, 2023, 620, 776–781.
15. A. Endo, K. Sato, K. Yoshimura, T. Kai, A. Kawada, H. Miyazaki and C. Adachi, Appl. Phys. Lett., 2011, 98, 083302.
16. E. W. Evans, Y. Olivier, Y. Puttisong, W. K. Myers, T. J. H. Hele, S. M. Menke, T. H. Thomas, D. Credgington, D. Beljonne, R. H. Friend and N. C. Greenham, J. Phys. Chem. Lett., 2018, 9, 4053–4058.
17. D. Hall, J. C. Sancho-García, A. Pershin, D. Beljonne, E. Zysman-Colman and Y. Olivier, J. Phys. Chem. A, 2023, 127, 4743–4757.
18. I. Lyskov and C. M. Marian, J. Phys. Chem. C, 2017, 121, 21145–21153.
19. J. M. Dos Santos, D. Hall, B. Basumatary, M. Bryden, D. Chen, P. Choudhary, T. Comerford, E. Crovini, A. Danos, J. De, S. Diesing, M. Fatahi, M. Griffin, A. K. Gupta, H. Hafeez, L. Hämmerling, E. Hanover, J. Haug, T. Heil, D. Karthik, S. Kumar, O. Lee, H. Li, F. Lucas, C. F. R. Mackenzie, A. Mariko, T. Matulaitis, F. Millward, Y. Olivier, Q. Qi, I. D. W. Samuel, N. Sharma, C. Si, L. Spierling, P. Sudhakar, D. Sun, E. Tankelevičiu̅tė, M. Duarte Tonet, J. Wang, T. Wang, S. Wu, Y. Xu, L. Zhang and E. Zysman-Colman, Chem. Rev., 2024, 124, 13736–14110.
20. F. B. Dias, T. J. Penfold and A. P. Monkman, Methods Appl. Fluoresc., 2017, 5, 012001.
21. H. Hirai, K. Nakajima, S. Nakatsuka, K. Shiren, J. Ni, S. Nomura, T. Ikuta and T. Hatakeyama, Angew. Chem., Int. Ed., 2015, 54, 13581–13585.
22. D. Hall, S. M. Suresh, P. L. Dos Santos, E. Duda, S. Bagnich, A. Pershin, P. Rajamalli, D. B. Cordes, A. M. Z. Slawin, D. Beljonne, A. Köhler, I. D. W. Samuel, Y. Olivier and E. Zysman-Colman, Adv. Opt. Mater., 2020, 8, 1901627.
23. S. Madayanad Suresh, D. Hall, D. Beljonne, Y. Olivier and E. Zysman-Colman, Adv. Funct. Mater., 2020, 30, 1908677.
24. T. Hatakeyama, K. Shiren, K. Nakajima, S. Nomura, S. Nakatsuka, K. Kinoshita, J. Ni, Y. Ono and T. Ikuta, Adv. Mater., 2016, 28, 2777–2781.
25. C. M. Marian, J. Phys. Chem. C, 2016, 120, 3715–3721.
26. P. K. Samanta, D. Kim, V. Coropceanu and J.-L. Brédas, J. Am. Chem. Soc., 2017, 139, 4042–4051.
27. Y. Olivier, J.-C. Sancho-Garcia, L. Muccioli, G. D'Avino and D. Beljonne, J. Phys. Chem. Lett., 2018, 9, 6149–6163.
28. T. J. Penfold, E. Gindensperger, C. Daniel and C. M. Marian, Chem. Rev., 2018, 118, 6975–7025.
29. K. Masui, H. Nakanotani and C. Adachi, Org. Electron., 2013, 14, 2721–2726.
30. N. Haase, A. Danos, C. Pflumm, A. Morherr, P. Stachelek, A. Mekic, W. Brütting and A. P. Monkman, J. Phys. Chem. C, 2018, 122, 29173–29179.
31. Y. Tsuchiya, S. Diesing, F. Bencheikh, Y. Wada, P. L. Dos Santos, H. Kaji, E. Zysman-Colman, I. D. W. Samuel and C. Adachi, J. Phys. Chem. A, 2021, 125, 8074–8089.
32. S. Mikulin, K. Bergmann, B. T. Luppi, S. A. Elgadi and Z. M. Hudson, J. Phys. Chem. B, 2025, 129, 13211–13220.
33. K. Shizu and H. Kaji, Commun. Chem., 2022, 5, 53.
34. K. Shizu and H. Kaji, Nat. Commun., 2024, 15, 4723.
35. N. Nijegorodov and R. Mabbs, Spectrochim. Acta, Part A, 2001, 57, 1449–1462.
36. Y. Tani and T. Ogawa, J. Mater. Chem. C, 2021, 9, 4281–4288.
37. K. Veys, M. H. E. Bousquet, D. Jacquemin and D. Escudero, J. Chem. Theory Comput., 2023, 19, 9344–9357.
38. C. A. Griffith, S. E. Krul, S. J. Hoehn, E. Mao, G. Sleyko and C. E. Crespo-Hernández, J. Phys. Chem. B, 2023, 127, 5924–5932.
39. N. Ikeda, S. Oda, R. Matsumoto, M. Yoshioka, D. Fukushima, K. Yoshiura, N. Yasuda and T. Hatakeyama, Adv. Mater., 2020, 32, 2004072.
40. S.-N. Zou, C.-C. Peng, S.-Y. Yang, Y.-K. Qu, Y.-J. Yu, X. Chen, Z.-Q. Jiang and L.-S. Liao, Org. Lett., 2021, 23, 958–962.
41. A. Pershin, D. Hall, V. Lemaur, J.-C. Sancho-Garcia, L. Muccioli, E. Zysman-Colman, D. Beljonne and Y. Olivier, Nat. Commun., 2019, 10, 597.
42. M. T. Do Casal, K. Veys, M. H. E. Bousquet, D. Escudero and D. Jacquemin, J. Phys. Chem. A, 2023, 127, 10033–10053.
43. N. J. Turro, Modern molecular photochemistry, University science books, 1991.
44. M. A. El-Sayed, J. Chem. Phys., 1963, 38, 2834–2838.
45. T. J. Penfold, E. Gindensperger, C. Daniel and C. M. Marian, Chem. Rev., 2018, 118, 6975–7025.
46. Y. Qian, T. Zhang, J. Han, A. R. Harutyunyan, G. Chen, Y. Rao and H. Chen, J. Phys. Chem. A, 2021, 125, 3589–3599.
47. M. Schreiber, M. R. Silva-Junior, S. P. A. Sauer and W. Thiel, J. Chem. Phys., 2008, 128, 134110.
48. A. Tajti and P. G. Szalay, J. Chem. Theory Comput., 2019, 15, 5523–5531.
49. A. Tajti, L. Tulipán and P. G. Szalay, J. Chem. Theory Comput., 2020, 16, 468–474.
50. R. R. Valiev, B. S. Merzlikin, R. T. Nasibullin, A. Kurtzevitch, V. N. Cherepanov, R. R. Ramazanov, D. Sundholm and T. Kurtén, Phys. Chem. Chem. Phys., 2023, 25, 6406–6415.
51. Y. Wang, J. Ren and Z. Shuai, J. Chem. Phys., 2021, 154, 214109.
52. J. C. Tully, J. Chem. Phys., 1990, 93, 1061–1071.
53. M. Barbatti, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2011, 1, 620–633.
54. D. Tang, C. Liao, M. Taub, G. C. Schatz, H. Guo and X. Li, J. Phys. Chem. Lett., 2025, 16, 3816–3821.
55. C. Liao, C. Rosenbaum, A. M. Glaudin, M. Taub, R. Banerjee Ghosh, S. Pristash, C. W. Schlenker and X. Li, J. Am. Chem. Soc., 2025, 147, 22176–22184.

---