Enhancing triplet harvesting in inverted singlet–triplet gap molecules through mechanistic understanding

Abstract

Molecules with an inverted singlet–triplet gap violate Hund's rule and offer promising applications in photoactive materials. In particular, they are expected to improve the efficiency of triplet harvesting via delayed fluorescence. Understanding the mechanistic factors governing electronic transitions is crucial for optimizing their performance. In this work, we employ the nuclear ensemble method combined with high-level electronic structure calculations to investigate the key excited-state processes in four azaphenalene-based molecules. We demonstrate that minimizing the negative singlet–triplet gap is, contrary to the expectations based on the Arrhenius relation, not necessarily a desired strategy for decreasing the delayed fluorescence lifetimes. Nevertheless, due to the vibronic coupling, molecular cores with small negative gaps are a promising starting point for assuring fast reverse intersystem crossing. These findings provide new design principles for developing efficient triplet-harvesting materials, emphasizing the need to control vibrational and symmetry effects to balance radiative and non-radiative transitions effectively.

---

1 Introduction

Singlet–triplet inversion of molecular excited states defies Hund's rule as the lowest singlet state (S₁) is energetically more stable than the lowest triplet state (T₁). This unique property holds promise for the applications of photoactive molecular materials that favor populating singlet states over triplets for maximizing the efficiency. In particular, materials exhibiting thermally activated delayed fluorescence (TADF) have been intensively studied for applications like organic light-emitting diodes, photocatalysis, and bioimaging. The key design parameter for TADF molecules is the minimization of the singlet–triplet gap (ΔE_ST), which has a positive value in the vast majority of cases. While the possibility of inversion via the spin polarization mechanism was theoretically predicted already in 1978,¹ the field remained mostly dormant until 2019 when two independent studies^2,3 confirmed this effect through ab initio calculations. It became clear that inversion results from nearly vanishing exchange interactions and additional stabilization of the singlet state through a higher contribution of double excitations. Since then, singlet–triplet inversion has been confirmed experimentally in molecules based on azaphenalene cores.^4,5

Considering the envisaged applications, it is crucial to understand how inversion affects the photophysical processes, where the forward and reverse intersystem crossing rates are the key parameters. Various computational approaches have been recently used to study the (reverse) intersystem crossing process in ST-inverted molecules. Karak et al.⁶ studied the temperature dependence of rISC rates in inverted molecules using double hybrid time-dependent density functional theory and the correlation-function formalism. They also analyzed the distribution of singlet–triplet gaps in ensembles of geometries generated via Wigner phase space sampling, highlighting that the inversion effect is strongly vibronically coupled. To account for vibrational effects, Dinkelbach et al.⁷ simulated (r)ISC in heptazine and HAP-3MF, incorporating Herzberg–Teller corrections. The rISC rate was found to depend strongly on ΔE_ST, with optimal rates at moderately negative S₁–T₁ energy gaps and a marked decrease as the gap widened. This result underscores the fact that large negative singlet–triplet energy gaps are generally not favorable for accelerating rISC rates. Izu et al.⁸ calculated the rISC rates in triangulene derivatives with a range of electronic structure methods, assuming the spin-vibronic mechanism and a Marcus-like rate expression. This study confirmed that the small SOC between S₁ and T₁ is significantly enhanced via the spin-vibronic mechanism. An explicit exploration of the S₁ and T₁ potential energy surfaces (PES) of azaphenalene cores was done by Valverde et al.⁹ The optimization of conical intersections and minimum energy crossing points (MECPs) allowed for estimating the activation energy barriers for internal conversion and intersystem crossing processes. The Herzberg–Teller expansion of spin–orbit couplings (SOC) showed pronounced vibronic enhancement of their magnitudes. However, the SOC vanishes at the singlet–triplet MECP, which raises questions about whether transitions are more likely to happen around MECPs or somewhere else on the PES.

This work aims to address outstanding questions about the potential advantages of azaphenalene-based molecules with negative singlet–triplet gaps for triplet harvesting by providing detailed insights into the mechanism of intersystem crossing and fluorescence. We focus on comprehensive modeling of vibrational effects via the nuclear ensemble method¹⁰ coupled to a high level of electronic structure theory (EOM-CCSD). A unique feature of the nuclear ensemble approach to calculating electronic transition rate constants is that it is based on direct sampling of the potential energy surfaces. It allows us to decompose the rates into the contributions from the configurations in the ensemble and identify those that contribute the most. Contrary to the harmonic and perturbative approximations, energies and couplings at those geometries are calculated exactly within the chosen electronic structure method. Also, it is not assumed that the transitions occur at the minimum energy crossing points between potential energy surfaces, but every sampled point has an associated weight. To generate an ensemble of molecular configurations, a large number of molecular geometries is sampled from a distribution for harmonic oscillators at the target temperature. Considering that the primary interest in ST-inverted systems is due to their potential applications in optoelectronic and photocatalytic materials, we also find it necessary to include condensed-phase effects. Coupling the nuclear ensemble method with a nonequilibrium implicit solvation model enables us not only to capture the solvatochromic effect on the energy levels but also to account for the outer-sphere reorganization energy in electronic transition rate constant calculations. We have recently demonstrated the usefulness of this approach to elucidate the thermally activated delayed fluorescence (TADF) mechanism in a carbene–copper emitter¹¹ and to study energy transfer in organic luminescent materials.^12,13 For the theoretical foundations and implementation details of our methodology, including the underlying assumptions and limitations, we refer to ref. 14 and15.

As the model systems for our study, we choose four related molecules illustrated in Fig. 1. The first two compounds are the prototypical inverted singlet–triplet systems, namely cycl[3.3.3]azine (molecule 1) and heptazine (molecule 2). The remaining two compounds are substituted heptazine derivatives, where its high D_3h symmetry is systematically lowered by replacing hydrogen atoms with methyl and amino groups. To this end, we use 2,5-dimethyl-heptazine (molecule 3) and 2-amino-5-methyl-heptazine (molecule 4), belonging to the C_2v and C_s symmetry point groups, respectively.

Fig. 1 Representation of the molecules studied 1 = cyclazine, 2 = heptazine, 3 = 2,5-dimethyl-heptazine and 4 = 2-amino-5-methyl-heptazine.

2 Methods

All geometry optimizations and frequency calculations were done using density functional theory (DFT) for the ground state and time-dependent DFT (TD-DFT) with the Tamm–Dancoff approximation for the excited states. We employed the long-range-separated hybrid functional ωB97-xD¹⁶ with the 6-31G(d,p) basis set as implemented in Gaussian 16.¹⁷

To generate ensembles and calculate electronic transition rate constants, we used the NEMO package¹⁸ coupled with the Q-Chem 6.0 quantum chemical software.¹⁹ For ensemble calculations, a total of 500 geometries were sampled from the Wigner distribution at 300 K for each of the S₀, S₁, and T₁ ensembles. For each geometry, single-point calculations were performed using the equation-of-motion coupled-cluster method (EOM-CCSD) with the cc-pVDZ basis set. Spin–orbit coupling calculations were conducted using the one-electron Breit–Pauli Hamiltonian. The non-equilibrium solvent correction is obtained using the perturbative state-specific (ptSS) solvation model²⁰ calculated at the TD-DFT level as it is unavailable for the EOM-CCSD method. Since the characters of the transitions (charge density differences) are very similar for both methods, using TD-DFT densities in the solvation protocol is justified. Unless otherwise stated, we have used parameters corresponding to toluene (ε = 2.38 and n = 1.497).

3 Results & discussion

For all the systems, the characters of the S₁ and T₁ states are qualitatively similar. The natural transition orbitals (NTOs) calculated at the ground-state geometry are generally spatially well separated (Fig. S1, ESI†). However, the LUNTO of molecule 4 has a pronounced amplitude also on atoms 3 and 7 (see Fig. 1) and on the amino group. Based on symmetry considerations,²¹ the S₀ → S₁ transition is dipole-forbidden for molecules 1 and 2, so we expect the corresponding oscillator strengths in the ensemble to be close to 0. While the C_2v symmetry of molecule 3 allows for a non-zero oscillator strength, the HONTO and LUNTO patterns are identical as in the D_3h symmetry molecules, so S₁ is also dark. On the contrary, a lower orbital symmetry in the molecule 4 should allow the transition to gain more intensity. Histograms of oscillator strength in the ground-state vibrational ensemble are shown in Fig. 2, along with their values at the S₀ minima and averages over the ensembles. The oscillator strength at the equilibrium geometry of molecule 4 is 0.014, while it is 0.0 for the higher-symmetry molecules. When vibrations are accounted for, the average oscillator strength is increased by 0.002 for molecules with a D_3h symmetry point group (i.e., molecules 1 and 2), and by 0.003 for C_2v and C_s point group symmetry molecules (molecules 3 and 4). The NTOs calculated at the first excited state equilibrium geometries exhibit the same pattern as at the ground-state equilibria (Fig. S3a, ESI†). Thus, the distributions of oscillator strengths in the S₁ ensemble are qualitatively similar (Fig. S2, ESI†). Again, vibrations make the transition weakly allowed by inducing the changes in the orbitals and increasing the overlap between HONTO and LUNTO (Fig. S3b, ESI†).

Fig. 2 Distribution of oscillator strengths for molecules 1 to 4 in the S₀ ensemble. Vertical dashed lines show the ensemble average whereas solid vertical lines show values for optimized structures.

Vibronic coupling alleviates the symmetry constraints in cyclazine and heptazine and makes the S₀ → S₁ transition weakly allowed. Indeed, as can be seen in (Fig. 3(a) and (b)), which compare experimental^22,23 and simulated S₀ → S₁ absorption bands in solution, the transition is now possible with the inclusion of vibrational effects. Though the ensemble simulations do not resolve the vibrational structure and the peak for heptazine is off by approximately 0.12 eV, the broadening of the absorption band is generally well captured.

Fig. 3 Simulated S₀ → S₁ absorption and S₁ → S₀ emission spectra of molecule 1 (a) and (c) and molecule 2 (b) and (d). Shaded areas in the simulated spectra depict statistical uncertainty due to the sampling of molecular configurations. Black dashed lines correspond to available experimental spectra of molecules 1²² and 2²³ and are shown for comparison.

The emission spectra of cyclazine and heptazine are shown in (Fig. 3(c) and (d)). The hypothetical emission from the S₁ state of cyclazine is predicted to peak at λ_em = 997 nm in hexane solution. In practice, cyclazine is an anti-Kasha molecule;²² therefore, the S₁ excited state decays nonradiatively and no emission is observed experimentally. The emission peak of 2 in acetonitrile is computed at λ_em = 494 nm, which corresponds to an absolute deviation of ΔE = 0.12 eV compared to the experimental value λ_max = 518 nm.²³ The simulated spectrum is broader than its experimental counterpart, including emission contributions below 500 nm. This mismatch is possibly due to the experimental spectrum being measured with an excitation wavelength of 436 nm, at the lower energy portion of the S₀ → S₁ absorption spectrum. It means that the excitation happens to lower vibronic states and, considering its rigidity, the molecule might not thermalize fully before emission. Since the simulations consider a thermalized population in S₁ as the initial state, the calculated spectrum includes higher energy contributions as well.

An overall good agreement between simulated and available experimental spectra reassures that our methodology is adequate for the investigation of electronic transitions in molecules based on azaphenalene cores. In particular, the nuclear ensemble approach allows us to represent the relevant regions of the potential energy surfaces well.

Fig. 4 shows distributions of fluorescence rate contributions from different geometries in the ensemble, as well as ensemble averages and values from equilibrium geometries. We notice that lowering the symmetry from D_3h to C_2v in molecule 3 allows for a non-vanishing fluorescence rate even at the equilibrium geometry (k_F = 2.50 × 10⁵ s⁻¹). For molecule 4, the corresponding equilibrium value is an order of magnitude larger (k_F = 7.9 × 10⁶ s⁻¹). Molecular vibrations shift the fluorescence rates to 9.0 × 10⁴ and 8.6 × 10⁵ s⁻¹ for molecules 1 and 2, respectively. Incidentally, the latter value is in an excellent agreement with the experimental radiative rate (5.1 × 10⁵ s⁻¹).²³ For 3, the increase in the rate constant, compared to the equilibrium contribution, is less pronounced as they differ only by a factor of 4. In contrast, for molecule 4, averaging over the ensemble does not cause almost any change in the fluorescence rate as the equilibrium and peak of the distribution coincide. Furthermore, the nuclear ensemble method gives a similar fluorescence rate for 2 (k_F = 8.6 × 10⁵ s⁻¹) as when using the Herzberg–Teller expansion (k_F = 3.0 × 10⁵ s⁻¹).⁷

Fig. 4 Histograms of contributions to the fluorescence rate constant k_F from the {S₁} ensemble geometries.

After discussing the absorption and fluorescence of molecules 1–4, we move to the intersystem crossing dynamics, for which ΔE_ST is the key parameter. The singlet–triplet gaps at the optimized S₁ geometries of cyclazine (1) and heptazine (2) are −0.11 eV and −0.21 eV, respectively, in close agreement with the gas-phase estimates available in the literature.^2,3,24 Dimethyl heptazine has ΔE_ST = −0.21 eV, i.e. the same as heptazine, while an asymmetric substitution with methyl and amino groups in molecule 4 raises the gap to ΔE_ST = −0.14 eV. These gaps are displayed in Fig. 5(a) along with ther ensemble distributions. When averaged over the {S₁} ensemble, ΔE_ST increases to −0.02 eV for molecule 1, with a similar effect being seen for molecules 2, 3, and 4, with average gaps moving to less negative values.

Fig. 5 Histograms of (a) ΔE_ST, (b) 〈S₁|H_SO|T₁〉, and (c) log k_rISC for molecules 1 to 4 in the {S₁} ensemble.

The second parameter determining the intersystem crossing dynamics is the spin–orbit coupling (SOC) between the S₁ and T₁ states. Analogously to the oscillator strengths, this transition appears forbidden in the D_3h symmetry, resulting in vanishing SOC at the equilibrium structures of 1 and 2 and in the meV range for 3 and 4. Molecular vibrations induce an increase in the ensemble-averaged SOCs to 0.030 meV in all four systems (Fig. 5(b)). A study using linear vibronic model at the TDA-DFT level found a very similar value for heptazine, but an order of magnitude lower for cyclazine.⁹ The discrepancy could result from the nonlinear vibrational effects and different electronic structure method. In any case, this increase in SOC, in conjunction with the variations in ΔE_ST, enables appreciable intersystem crossing rates in the order of 10⁷ s⁻¹ (Fig. 5(c)). Again, the rate constant for heptazine (k_ISC = 2.2 × 10⁷ s⁻¹) is in good agreement with the value obtained using Herzberg–Teller expansion (k_ISC = 1 × 10⁷ s⁻¹).⁷

Fig. 6 shows the relationship between the ΔE_ST, SOC, and the contribution to the S₁ → T₁ ISC rate for every configuration in the ensemble. Unsurprisingly, the configurations that have the highest ISC rates (marked with a red cross) are those with ΔE_ST ≈ 0. However, such configurations do not necessarily exhibit the highest SOCs within the ensemble. Interestingly, the distribution of ISC rate contributions versus the energy gap varies across the molecules. While molecule 1 displays a symmetrical and evenly distributed pattern, molecules 2 and 3 show a broader distribution with only a few configurations having a positive energy gap. In contrast, molecule 4 exhibits a distribution more similar to molecule 1.

Fig. 6 Distribution of the k_ISC rates versus the ΔE_ST gap and the SOC between the considered excited states displayed on the color scale. The red “X” corresponds to the geometry that contributes the most to the rate. Its properties are specified on the top left for each molecule.

Since each configuration in the ensemble is generated by sampling the normal modes of the system, we can analyze the relationships between vibrational modes and electronic properties. Correlating the absolute displacement of nuclei with ΔE_ST and SOC for each vibrational mode reveals that only one mode for each molecule is significantly correlated (Fig. S5, ESI†). Specifically, the vibrational modes ω₅₁ = 2250 cm⁻¹, ω₃₉ = 1767 cm⁻¹, ω₅₃ = 1783 cm⁻¹, and ω₅₁ = 1737 cm⁻¹ for 1, 2, 3, and 4, respectively, all correspond to in-plane anti-symmetric stretching vibrations of peripheral carbon atoms (for 1) and nitrogen atoms (for the others) (see Fig. S6, ESI†). The correlation weakens as the molecular symmetry decreases, suggesting that for more complex systems, additional factors must be considered alongside those presented here.

We have employed the same procedure to analyze the rISC mechanism, which requires repeating calculations in the T₁ ensemble as it is the initial state. Despite rISC being a thermodynamically downhill process in inverted-ST molecules, it is still thermally activated due to an energy barrier on the T₁ potential energy surface. The distributions of key parameters are displayed in Fig. 7. Among the four molecules, cyclazine exhibits the same singlet–triplet gap for both S₁ and T₁ states at their respective equilibrium geometries (Fig. 7(a)). In contrast, heptazine and molecule 3 show an increase in ΔE_ST to −0.18 eV and −0.15 eV, respectively. The most pronounced variation occurs in molecule 4, where the optimized gap approaches zero. The ensemble-averaged gap of cyclazine (ΔE_ST = 0.01 eV) shows the largest deviations from the optimized-geometry value. Molecules 2 and 3 retain the negative gaps of ΔE_ST = −0.11 eV, while molecule 4 once again exhibits the smallest variation, with ΔE_ST = 0.03 eV. The distribution of spin–orbit coupling constant between the T₁ and S₁ states is depicted in Fig. 7(b). The qualitative picture is the same as in the S₁ ensemble, but the ensemble averages are a bit higher (by 0.001–0.011 meV).

Fig. 7 Histograms of (a) ΔE_ST, (b) 〈T₁|H_SO|S₁〉 and (c) the distribution of log k_rISC in which the rate is expressed in s⁻¹ for molecules 1 to 4 from the {T₁} ensemble geometries.

The distribution of rISC rate contributions is shown in Fig. 7(c) for all four systems. As in the case of ISC, the contribution from the equilibrium geometry is negligible for 1 and 2; however, it corresponds to k_rISC = 3.4 × 10⁶ s⁻¹ for molecule 4. The latter results from a combination of a smaller |ΔE_ST| and nonvanishing SOC at the T₁ equilibrium geometry. Vibrations induce a strong increase in the rate constants, which are in the order of 3.5–5.6 × 10⁷ s⁻¹ for all systems. Contrary to the S₁ ensemble case, in the T₁ ensemble, only molecule 1 displays a direct correlation between electronic parameters and the displacement along any of the vibrational modes.

Analyzing the relationship between ΔE_ST, SOC and k_rISC depicted in Fig. 8, it is interesting to note that the SOC tends to increase when ΔE_ST is the most negative. Also, the configuration that contributes the most to the rate is not the one with the highest SOC nor the one with the lowest singlet–triplet gap. Indeed, the latter is close to 0 or slightly positive for 4 while high SOC geometries have gaps between 0.06 meV to 0.18 meV.

Fig. 8 Distribution of the k_rISC rates versus the ΔE_ST gap and the SOC between the considered excited states displayed on the color scale. The red “X” corresponds to the geometry that contributes the most to the rate. Its properties are specified on the top left for each molecule.

From Fig. 6 and 8, it is clear that there are multiple configurations in the ensemble with nearly degenerate S₁ and T₁ states and significant contributions to the ISC and rISC rate constants. Analysis of the configurations that account for 95% of the rate constant shows that there is no correlation between their potential energy and the magnitude of contribution. It is often assumed that electronic transitions are most likely to occur in the vicinity of minimum energy crossing points, whose energy is the effective activation energy. This assumption is particularly questionable for azaphenales due to the vanishing SOC at the MECP,⁹ and raises a question of what the effective activation energy for the transition is. The potential energies of configurations with significant contributions to the rate (Fig. S7, ESI†) can be as high as a few eV, an order of magnitude higher than the reported MECP energies (∼0.1–0.4 eV), which indicates that the majority of (r)ISC transition events happen far away from the MECP.

For each of the studied molecules, we have identified the configuration in the S₁/T₁ ensemble that contributes the most to the ISC/rISC rate constant. For each configuration, we calculated the zero-point vibrational energy (E_ZPV) at the corresponding equilibrium structure, potential energy under the harmonic approximation (V_harmonic), ab initio potential energy (V_abinitio), and the minimal thermal energy needed to reach this configuration using the harmonic approximation (E_thermal) (Table 1). The latter quantity is calculated as the total excess energy above the zero-point energy for each of the normal modes: , where Q_i is the displacement along the ith vibrational mode and μ_i the corresponding reduced mass. As such, E_thermal is the minimum thermal energy that needs to be deposited in the vibrational motion to go beyond the classical turning points in the ground vibrational state and reach the analyzed configuration.

Table 1 Ab initio and harmonic potential energies, zero-point vibrational energies, and lower bounds for the thermal energy in the harmonic approximation for configurations with the highest contributions to the ISC and rISC rates for molecule 1–4. All values in eV

Molecule	V abinitio	V harmonic	E ZPV	E thermal
ISC
1	3.23	3.44	4.85	1.30
2	1.76	1.88	2.85	0.77
3	3.21	3.12	4.35	0.97
4	3.42	1.69	4.06	0.29
rISC
1	1.69	1.87	4.78	0.34
2	1.93	1.75	2.79	0.75
3	4.57	3.03	4.28	0.93
4	2.57	2.27	4.01	0.60

---

There is an overall good agreement between harmonic and ab initio potential energies, except for two configurations (molecule 4 in the S₁ ensemble and molecule 3 in the T₁ ensemble), for which anharmonicity or mode coupling appear to be significant. In all cases, total potential energy is well below the zero-point vibrational energy. It does not mean, of course, that the configurations are accessible at 0 K, as the energy cannot be exchanged between different modes in the ground state. Indeed, E_thermal, which can be loosely associated with an effective activation energy for the specific configuration, varies between 0.3 and 1.3 eV. It is notably higher than the MECP energy, confirming that significant transitions are happening far away from it. More importantly, heptazine at MECP is planar and has a vanishing SOC, while the analyzed configuration is distorted (Fig. S8, ESI†), and simultaneous thermal activation of multiple modes is necessary.

The simulations show that it is possible to have an inverted singlet–triplet molecule with an appreciable radiative rate constant and rISC faster than ISC. The magnitudes of these rate constants are still too low for practical applications as TADF emitters but serve as a proof of principle for further development. At the same time, k_rISC > k_ISC together with a slow radiative rate constant is a promising prerequisite for photocatalysts, as it offers a pathway for long-lived excitations with suppressed triplet quenching by molecular oxygen. As was already pointed out by Marian and collaborators,⁷ large negative gaps are not necessarily favorable for speeding up rISC. We see it also from our simulations gathered in Table 2, where the comparison of the molecules 2 and 4 show that the latter has faster rISC despite a less negative gap. Nevertheless, larger negative gaps may still be favorable for maintaining higher populations of singlet excitons, higher quantum yields, and lower radiative lifetimes, which would depend on the interplay between all the possible radiative and nonradiative transitions.

Table 2 Comparison of k_ISC and k_rISC for the four molecules

Molecule	k ISC (× 10^7 s^−1)	k rISC (× 10^7 s^−1)
1	4.1 ± 0.4	3.5 ± 0.3
2	2.2 ± 0.5	3.5 ± 0.6
3	1.6 ± 0.4	5 ± 1
4	5.1 ± 0.9	5.6 ± 0.5

---

To analyze the overall population dynamics after an excitation, we computed the time evolution of the S₁ and T₁ state populations as well as time-resolved photoluminescence (TRPL) spectra under optical excitation (Fig. S9a–h, ESI†) using a kinetic model detailed in the ESI.† By biexponential fitting of the TRPL spectra, we extract the amplitudes and time scales of prompt and delayed fluorescence (Table 3). The prompt fluorescence lifetime of the first three systems is within the range of 13 to 17 ns, while molecule 4 has the shortest lifetime with τ_PF = 9 ns. Nevertheless, all prompt fluorescence lifetimes have the same order of magnitude as they are primarily determined by the intersystem crossing dynamics.

Table 3 Fitting parameters for the prompt and delayed fluorescence lifetimes. Computed with EOM-CCSD/cc-pVDZ in toluene

Molecule	A PF	τ PF (s)	A DF	τ DF (s)
1	56	1.32 × 10^−8	47	2.39 × 10^−5
2	39	1.74 × 10^−8	64	1.86 × 10^−6
3	19	1.50 × 10^−8	82	1.19 × 10^−6
4	51	8.97 × 10^−9	54	2.53 × 10^−7

---

On the other hand, delayed fluorescence lifetimes differ by two orders of magnitude, with the shortest for molecule 4 and the longest for 1. This trend aligns well with the computed fluorescence rates, where k_F is higher for 4 than for 1. This observation means that singlet–triplet gap inversion, and consequently k_rISC > k_ISC, is of little significance for the effectiveness of triplet harvesting through delayed fluorescence in this class of molecules. Similar delayed lifetime results are shown for simulations under an electrical, rather than optical, excitation (Fig. S10a–h, ESI†).

4 Conclusion

In this paper, we calculated and analyzed the main rate constants determining the efficiency of triplet harvesting in four azaphenalene-based inverted singlet–triplet gap molecules. To account for the dielectric condensed-phase effect on the rate constants, we used a state-specific implicit solvation model. Using the nuclear ensemble method, we captured the vibrational effects on rate constants and decomposed them into contributions from different configurations on the relevant potential energy surfaces. The electronic structure of the configurations in the ensembles was calculated at a high level of theory (EOM-CCSD).

Our work sheds new light on the design of inverted singlet–triplet molecules as triplet harvesters. While Hund's rule violation is fundamentally intriguing, minimizing the (negative) singlet–triplet gap to enhance rISC thermodynamic driving force appears counterproductive from an application standpoint. While it may lead to the rISC rate constant being higher than for ISC, the effect is rather small and can defeat the overall goal by concomitantly decreasing both rates. As long as photon emission is the bottleneck, the focus should be on enhancing the fluorescence rate constant. Indeed, our simulations show that functionalizing heptazine can decrease the delayed fluorescence lifetime by an order of magnitude despite pushing the singlet–triplet gap toward positive values. The latter observation also highlights the fact that non-radiative electronic transitions occur at configurations where the states are quasi-degenerate, rather than when the final state is energetically below the initial.

On the other hand, a promising observation for the inverted gap molecules is that ensemble-averaged gaps and spin–orbit couplings are generally higher than calculated at the equilibrium geometries. Therefore, not only do statistically more configurations get closer to the resonance condition, but also the coupling becomes stronger. The main conclusion is that molecular cores with small negative gaps at equilibrium geometries are indeed interesting starting points for further molecular design. Like in the case of molecule 4, the effect of functional groups and molecular vibrations can lead to the ensemble-averaged singlet–triplet gap being close to zero and many configurations contributing significantly to (r)ISC. We also demonstrated that a tacit assumption that the electronic transitions are most likely to happen around minimum energy crossing points is, in this case, not valid. We show that the most effective transitions happen at geometries that require thermal activation of multiple modes at the same time, and the estimated activation energy is significantly higher than at MECP.

While this study focuses on (r)ISC and emissive rate constants, non-radiative relaxation to the ground state can also significantly affect triplet harvesting efficiency. It has been recently shown that suppression of some vibrations via aggregation is necessary to avoid rapid internal conversion in heptazine derivatives.²⁵ Altogether, these results point to the need for a thorough understanding and control of vibrational effects in inverted singlet–triplet gap materials, both at the level of single molecules and molecular aggregates.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the ESI.†

Acknowledgements

This work was supported by a research grant (VIL50199) from VILLUM FONDEN, grant no. 2032-00144B from the Independent Research Fund Denmark, and grant from the Novo Nordisk Foundation Data Science Research Infrastructure 2022 Grant: A high-performance computing infrastructure for data-driven research on sustainable energy materials, Grant no. NNF22OC0078009.

References

1. H. Kollmar and V. Staemmler, Violation of Hund's rule by spin polarization in molecules, Theor. Chim. Acta, 1978, 48, 223–239.
2. P. de Silva, Inverted Singlet-Triplet Gaps and Their Relevance to Thermally Activated Delayed Fluorescence. The, J. Phys. Chem. Lett., 2019, 10, 5674–5679.
3. J. Ehrmaier, E. J. Rabe, S. R. Pristash, K. L. Corp, C. W. Schlenker, A. L. Sobolewski and W. Domcke, Singlet-Triplet Inversion in Heptazine and in Polymeric Carbon Nitrides, J. Phys. Chem. A, 2019, 123, 8099–8108.
4. N. Aizawa, Y.-J. Pu, Y. Harabuchi, A. Nihonyanagi, R. Ibuka, H. Inuzuka, B. Dhara, Y. Koyama, K.-I. Nakayama, S. Maeda, F. Araoka and D. Miyajima, Delayed fluorescence from inverted singlet and triplet excited states, Nature, 2022, 609, 502–506.
5. K. D. Wilson, W. H. Styers, S. A. Wood, R. C. Woods, R. J. McMahon, Z. Liu, Y. Yang and E. Garand, Spectroscopic Quantification of the Inverted Singlet–Triplet Gap in Pentaazaphenalene, J. Am. Chem. Soc., 2024, 146, 15688–15692.
6. P. Karak, P. Manna, A. Banerjee, K. Ruud and S. Chakrabarti, Reverse Intersystem Crossing Dynamics in Vibronically Modulated Inverted Singlet-Triplet Gap System: A Wigner Phase Space Study, J. Phys. Chem. Lett., 2024, 15, 7603–7609.
7. F. Dinkelbach, M. Bracker, M. Kleinschmidt and C. M. Marian, Large Inverted Singlet-Triplet Energy Gaps Are Not Always Favorable for Triplet Harvesting: Vibronic Coupling Drives the (Reverse) Intersystem Crossing in Heptazine Derivatives, J. Phys. Chem. A, 2021, 125, 10044–10051.
8. A. E. Izu, J. M. Matxain and D. Casanova, Reverse intersystem crossing mechanisms in doped triangulenes, Phys. Chem. Chem. Phys., 2024, 26, 11459–11468.
9. D. Valverde, C. T. Ser, G. Ricci, K. Jorner, R. Pollice, A. Aspuru-Guzik and Y. Olivier, Computational Investigations of the Detailed Mechanism of Reverse Intersystem Crossing in Inverted Singlet–Triplet Gap Molecules, ACS Appl. Mater. Interfaces, 2024, 16, 66991–67001.
10. R. Crespo-Otero and M. Barbatti, Spectrum simulation and decomposition with nuclear ensemble: formal derivation and application to benzene, furan and 2-phenylfuran, Theor. Chem. Acc., 2012, 131, 1237.
11. L. de Thieulloy, L. E. de Sousa and P. de Silva, TADF mechanism in a carbene-copper emitter: Insights from the nuclear ensemble simulations, J. Phys. Chem. C, 2024, 128, 14887–14896.
12. L. E. de Sousa, L. dos Santos Born, P. H. de Oliveira Neto and P. de Silva, Triplet-to-singlet exciton transfer in hyperfluorescent OLED materials, J. Mater. Chem. C, 2022, 10, 4914–4922.
13. F. T. Bueno, T. d S. A. Cassiano, P. de Silva, P. H. de Oliveira Neto and L. E. de Sousa, Exploring the triplet-to-singlet conversion mechanism in persistent luminescence: insights from a host–guest system, J. Mater. Chem. C, 2025, 13, 2673–2680.
14. L. E. de Sousa and P. de Silva, Unified framework for photophysical rate calculations in TADF molecules, J. Chem. Theory Comput., 2021, 17, 5816–5824.
15. L. E. de Sousa and P. de Silva, Photophysics of Solvated Molecules: Computational Protocol Combining Nuclear Ensemble and Nonequilibrium State-Specific Solvation Methods, J. Phys. Chem. A, 2023, 127, 8200–8208.
16. J.-D. Chai and M. Head-Gordon, Long-range corrected hybrid density functionals with damped atom-atom dispersion corrections, Phys. Chem. Chem. Phys., 2008, 10, 6615–6620.
17. M. J. Frisch, et al., Gaussian∼16 Revision C.01., Gaussian Inc., Wallingford CT, 2016.
18. Nemo - photophysics with the nuclear ensemble method. https://github.com/M3G-DTU/NEMO, 2023; Accessed on March 1st, 2024.
19. E. Epifanovsky, et al., Software for the frontiers of quantum chemistry: An overview of developments in the Q-Chem 5 package, J. Chem. Phys., 2021, 155, 084801.
20. J.-M. Mewes, Z.-Q. You, M. Wormit, T. Kriesche, J. M. Herbert and A. Dreuw, Experimental benchmark data and systematic evaluation of two a posteriori, polarizable-continuum corrections for vertical excitation energies in solution, J. Phys. Chem. A, 2015, 119, 5446–5464.
21. G. Ricci, J.-C. Sancho-García and Y. Olivier, Establishing design strategies for emissive materials with an inverted singlet–triplet energy gap (INVEST): a computational perspective on how symmetry rules the interplay between triplet harvesting and light emission, J. Mater. Chem. C, 2022, 10, 12680–12698.
22. W. Leupin and J. Wirz, Low-lying electronically excited states of cycl[3.3.3]azine, a bridged 12.pi.-perimeter, J. Am. Chem. Soc., 1980, 102, 6068–6075.
23. A. M. Halpern, M. A. Rossman, R. S. Hosmane and N. J. Leonard, Photophysics of the S₁ ↔ S₀ transition in tri-s-triazine, J. Phys. Chem., 1984, 88, 4324–4326.
24. R. Pollice, P. Friederich, C. Lavigne, G. dos Passos Gomes and A. Aspuru-Guzik, Organic molecules with inverted gaps between first excited singlet and triplet states and appreciable fluorescence rates, Matter, 2021, 4, 1654–1682.
25. L. de Thieulloy, R. S. Oliboni, P. de Silva and L. G. C. Rego, Aggregation Induced Effects on the Nonradiative Recombination Dynamics of Inverted Singlet–Triplet Heptazine-Based Materials, J. Phys. Chem. A, 2025, 129, 5220–5233.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5tc01155h

---