Effective models for TADF: the role of the medium polarizability

Abstract

A novel approch to estimate intersystem and reverse intersystem crossing rates (ISC and RISC rates, respectively) is proposed. We build on an effective model Hamiltonian recently parametrized ab initio and validated experimentally for a prototypical dye for thermally activated delayed fluorescence (TADF). The model describes the relevant physics in terms of a few diabatic states coupled to an effective vibrational coordinate and an effective conformational mode. A complete, numerically exact, non-adiabatic solution of the problem opens the way to the calculation of ISC and RISC rates fully accounting for the anharmonic and non-adiabatic nuclear dynamics. The model is further extended to address the role of the environmental polarizability, as described by the medium refractive index. The marginal variability of the refractive index in organic media results in marginal effects on the rates in different media. However, large variations of the rates are predicted when moving from the gas phase to an organic medium, suggesting that a meaningful analysis of experimental data must rely on computational analysis properly accounting for the dielectric properties of the surrounding medium.

---

1 Introduction

Thermally-activated delayed fluorescence (TADF) is a delicate process driven by a complex interplay of several interactions.^1–3 Electronic excited states with different nature, singlets and triplets, on one side, and charge transfer (CT) and local excited (LE) states, on the other side, are involved in a process that is also affected by vibrational and conformational motion.^4–8 Tiny spin–orbit coupling (SOC) interactions must be understood in this complex scenario and a wealth of experimental and theoretical work is devoted to unveiling the relevant physics.^9–19 From a theoretical perspective, the detailed understanding of TADF, as needed to actually simulate the time-dependent phenomena, is highly non-trivial. Indeed the non-adiabatic coupling between electronic and vibrational/conformational degrees of freedom must be accounted for in a system characterized by large anharmonicity.⁹ Moreover, a careful analysis of the different timescales of the competing phenomena (internal relaxation, conformational motion, non-radiative and radiative relaxations, inter-system crossing (ISC) and reverse inter-system crossing (RISC)) must be carefully analyzed and accounted for.

The comparison between theory ad experiment is further hindered by the subtleties introduced in the system by environmental effects. TADF as relevant to real devices in fact occurs in matrices, typically amorphous (glassy) matrices composed of small molecules. The dielectric properties of the matrix, its polarity and polarizability affect in different ways LE and CT states, leading to qualitatively different TADF properties in different media and making the comparison with theoretical gas-phase results unreliable.^{15,16,20–26} The rigidity and viscosity of the medium also play a role, affecting in turn the molecular conformation and flexibility.²⁵ Different theoretical approaches have been exploited to account for environmental effects in TADF processes.^{15,16,23–26} They can be classified in two broad families. Continuum models (typically conductor-like screening model (COSMO) and polarizable continuum model (PCM)) describe the environment as a continuum dielectric medium, then loosing all information about the atomistic nature of the matrix, its rigidity and viscosity, but leading to a simple and economic picture for dielectric effects.^27–29 In atomistic models, a detailed and dynamical model for the environment is defined, typically in terms of molecular mechanics or molecular dynamics models, that, at a larger computational cost than continuum models, allows for a more realistic description of environmental effects.^30,31 However, both families of approaches fail to properly account for the effect due to the medium polarizability.^32–34 Indeed, in both approaches the molecular Hamiltonian is defined as responding to the potential due to the charges on the surrounding medium, in the hypothesis that these charges are equilibrated to the specific molecular state. This adiabatic treatment of environmental effects of course works well for environmental charges related to slow degrees of freedom (vibrational, conformational and orientational motions in the matrix), but it is untenable for the fast degrees of freedom associated with the matrix polarizability. The electronic degrees of freedom of the matrix molecules (typically in the UV region) are characterized by faster timescales than the relevant degrees of freedom of the TADF dyes (typically in the visible) and should rather be treated in the antiadiabatic approximation.³²

Here a model, recently parametrized ab initio and validated experimentally for a prototypical TADF dye (Fig. 1), 9,9-dimethyl-9,10-dihydroacridin-4,6-triphenyl-1,3,5-triazine (DMAC–TRZ),³⁵ is used to propose a novel approach to RISC, ISC, and radiative decay rates, fully accounting for non-adiabatic and anharmonic effects. Gas phase results will be presented first. Then we will address the effect of the matrix polarizability, showing how it affects the behavior of the system. Results will be compared with experimental data collected in a non-polar environment.²² The discussion of medium polarity effects will be deferred to a subsequent publication.

Fig. 1 Top: Chemical structure of DMAC–TRZ. Left panel: A schematic view of the four diabatic states and of the relevant mixing matrix elements. Right panel: The four electronic adiabatic states, the arrows mark the relaxation processes of interest.

2 Setting up the model: TADF in the gas phase

2.1 The model

DMAC–TRZ is a typical TADF dye, made up by an electron acceptor (A, TRZ) and an electron donating (D, DMAC) group joined by a poorly conjugating bridge.³⁶ The dye has been extensively investigated experimentally and we recently proposed a minimal model to describe its low-energy physics.^22,35,36 The model, parametrized against extensive time-dependent density functional theory (TD-DFT) calculations, was validated against spectroscopic data. It accounts for four diabatic states: two singlet states, corresponding to a neutral |N〉 state, that may be represented as DA, and a zwitterionic state |Z〉, D⁺A⁻. While these two states are enough to address optical spectra of DA dyes, TADF also requires triplet states: we will account for the zwitterionic triplet |T〉 and a local excited triplet |L〉. The energy of the diabatic zwitterionic states is set to 2z (the N state sets the energy zero) and 2k is the energy of the local triplet state. The mixing matrix elements of the two singlet (triplet) states is τ. Finally, we neglect the Z–T and N–L SOC (El-Sayed rule), and set the N–T and Z–L couplings to V_SOC and W_SOC, respectively.³⁷ The relevant electronic Hamiltonian reads (N, Z, T, L basis):

A schematic picture of the diabatic energy levels and their mixing matrix elements is reported in the left panel of Fig. 1. The diagonalization of the electronic model leads to the four adiabatic states in the right panel of Fig. 1. These states are named following the standard quantum chemistry notation, S₀ is the lowest singlet (the ground state), S₁ the lowest excited singlet. Triplet states are named T₁, T₂, etc. The arrows mark the radiative and non-radiative relaxation paths of interest.

In eqn (1) we have introduced the dependence of some terms on a conformational and a vibrational coordinate (δ and Q, respectively). The conformational motion describes the torsion around the DA bond. The equilibrium ground state geometry is characterized by an orthogonal configuration of the DMAC and TRZ planes. The torsional coordinate δ measures the deviations of the torsional angle from the orthogonal position. The δ-dependence of the mixing matrix elements is simulated as τ(δ) = τ₀|sin δ| and β(δ) = β₀|sin δ|. An anharmonic (quartic) restoring potential is introduced to mimic the ground and excited state potential energy surfaces, calculated via TD-DFT. Finally, the effective molecular vibration Q accounts for the stabilization of the zwitterionic diabatic states upon geometry relaxation, so that where ε_v is the relaxation energy and ω_v the vibrational frequency.³⁸ The complete Hamiltonian reads:

where ω_c/v are the frequencies associated with the conformational/vibrational coordinates, the momenta associated with the dimensionless conformational and vibrational coordinate and , respectively.

Details of the model parametrization can be found in ref. 35, shortly, TD-DFT (Tamm Dancoff approximation, M062x functional, 6-31G(d) basis) was exploited to calculate the δ-dependence of the first few excited states in the singlet and triplet manifold. The first eight parameters in Table 1 were fixed to best reproduce the resulting PES, as well as the δ-dependence of the SOC matrix elements (ZORA). The vibrational relaxation energy ε_v, measuring the energy gained upon geometry rearrangement when going from the neutral to the zwitterionic state, was estimated summing up the relaxation energies associated with the ionization of the separated D and A fragments. The effective vibrational frequency ω_v was fixed to best reproduce the shape of absorption and emission spectra. As discussed in ref. 35, V_SOC and W_SOC were fixed as to best simulate the δ dependence of the calculated S₁–T₁ and S₀–T₁ SOC. We explicitly checked (see ESI,† Fig. S1) that the best fit of calculated couplings is obtained when V_SOC and W_SOC have the same sign (of course multiplying both terms by −1 does not change the results, but the relative sign of the two terms has important implications, as discussed below).

Table 1 Model parameters (eV) from ref. 35

z	τ 0	k	β 0	Vsoc	Wsoc	ω c	a	ω v	ε v
1.72	0.75	1.96	0.85	3.84 × 10^−4	1.74 × 10^−4	2.40 × 10^−4	1.43 × 10^−7	0.18	0.17

---

Finally, important information is obtained from the evolution of the energy of the ground and the lowest excited states on an electric field applied along the D–A direction: on one side the distinctively different F-dependence of the energies of the states allows to unambiguously assign them as local or CT states. Moreover, the slope of the quasi-linear dependence of the energy of the CT states on the field, offers a reliable estimate of the dipole moment associated with the zwitterionic state as μ₀ = 22.7 D. This quantity enters into the definition of the dipole moment operator, , needed to calculate radiative rates and to account for the interaction of the dye with the dieletric environment. The electric dipole is an electronic operator, diagonal over the vibrational and conformational states. Neglecting all matrix elements of the dipole moment operator on the electronic basis, but μ₀, the dipole moment operator reads:

In ref. 35 we discussed the solution of the Hamiltonian in eqn (2) treating the vibrational coordinate in a non-adiabatic approach, while the conformational coordinate is treated in the adiabatic approximation as a classical coordinate. In other terms, we neglected the kinetic energy associated with the δ coordinate and solved the coupled electronic-vibrational Hamiltonian for fixed δ values. This approximation is perfectly adequate to address optical spectra, since relevant electronic energies are orders of magnitude larger than conformational energies. However, the tiny singlet–triplet gap of TADF dyes makes the adiabatic treatment of the conformational degree of freedom untenable to address RISC and ISC processes.

Here we set up a full non-adiabatic calculation, writing the Hamiltonian in eqn (2) on the basis obtained as the direct product of the four electronic diabatic states times the eigenstates of the harmonic oscillator associated with the vibrational coordinate times the eigenstates of the harmonic part of the potential associated with the conformational motion. Of course, the infinite harmonic oscillator basis must be truncated to a large enough number of states as to get convergence on the quantities of interest. For vibrational states, the relevant number of states is reasonably small (typically 10–20), but the very small frequency associated with the conformational mode requires using a very large number of states, 500–1000, leading to a large overall dimension. Moreover, we need a polynomial expansion of sin δ, and, in line with the quartic expansion of the potential, we truncate the expansion to the third order (see ESI† for technical details).

2.2 RISC and ISC rates

Having a model for DMAC–TRZ and a complete non-adiabatic solution of the relevant Hamiltonian, we now face the problem of calculating relaxation rates. The most elegant way would rely on open quantum systems approaches, coupling the molecular system to a thermal bath, to simulate the energy degradation of the system due to the interaction with the environment. The Redfield model, would nicely do the job, allowing to estimate all relaxation energies in a single shot, as discussed in ref. 39. However, the presence of a vibrational and a conformational coordinate, with distinctively different frequencies would require the introduction of two different coupling channels between the molecule and the bath, one for each coordinate. Accordingly, two different spectral densities should be introduced: the relative strength of the couplings and the specific form of the two spectral densities would affect calculated rates in an uncontrolled way. We therefore proceed in a step by step approach: here we discuss how we estimate RISC and ISC rates and in the next section we will address radiative rates.

RISC and ISC processes are driven by tiny SOC interactions that can be treated perturbatively. The unperturbed states are therefore obtained from the non-adiabatic diagonalization of the molecular Hamiltonian in eqn (2) but setting V_SOC and W_SOC to zero. In these conditions, the singlet and triplet subspaces are decoupled and the two problems can be treated separately. Red and black lines in Fig. 2 show the energy levels relevant to the vibronic states in the T₁ and S₁ manifolds, respectively. Since internal conversion is a very fast process (with typical relaxation times ∼100 fs) we assume that RISC occurs from the thermally equilibrated T₁ states, the relevant distribution being shown in the right panel of Fig. 2. Then we use the Fermi golden rule (FGR) to calculate the rate of transition between states i and j in the two subspaces:

where Ĥ_SOC is the SOC part of the Hamiltonian in eqn (2). A finite width must be associated with each state, and S_ij measures the overlap between the two states. Each state is assigned a Gaussian lineshape, whose width σ is related to the inverse relaxation time τ_r as follows:

Fig. 2 A schematic view of the non-adiabatic calculation of ISC (left panel) and RISC (right panel). In both panels, red and black lines show the energy of the vibronic eigenstates of T₁ and S₁, respectively. In the left panel we show the Gaussian shape (not to scale) assigned to a specific vibronic state in the S₁ manifold and the arrows indicate the transition to specific vibronic eigenstates in the T₁ manifold. The global ISC rate is calculated summing on all S₁ eigenstates, accounting for their thermal population (black dashed line in the left panel). RISC rates (right panel) are calculated in a similar way, but summing on all T₁ eigenstates, accounting for their thermal population (red dashed line in the right panel).

ISC rates are calculated along the same lines, but starting with vibronic eigenstates in the S₁ manifold, whose thermal distribution is shown as a black dashed line in the left panel.

Fig. 3 shows RISC and ISC rates calculated as a function of temperature for three different values of the lifetime of the vibronic eigenstates: 50 fs, 100 fs and 200 fs, spanning the relevant range. Calculated rates increase a little with decreasing relaxation time (increasing the width of the energy interval assigned to each vibronic state), but the observed variation is modest. We therefore stick on results obtained for the intermediate 100 fs relaxation time. As expected, the RISC rate decreases fast upon decreasing temperature, while ISC rate moderately increases. The singlet triplet gap, ΔE_ST, defined as the energy difference between the lowest vibronic state in each subspace, amounts to 0.0513 eV.

Fig. 3 Temperature dependence of calculated RISC and ISC rates. Blue, black and red symbols refer to results obtained imposing a relaxation time for vibronic eigenstates τ_r = 50 fs, 100 fs and 200 fs, respectively. Black open symbols show RISC rates calculated from ISC rates, upon imposing the reversibility condition.

The microscopic reversibility condition relates the RISC and ISC rates: the dashed line in Fig. 2 shows the RISC rate evaluated from the ISC rate upon imposing the microscopic reversibility condition:

where k_B is the Boltzmann constant. Our results marginally deviate from the reversibility condition, suggesting that, while the proposed model is very simple, as it only accounts for a single low-frequency mode, the spacing between vibronic level is small enough to accommodate for thermodynamic behavior.

We can play around with model parameters to understand the effect of specific interactions. Fig. 4 shows the temperature dependence of the RISC and ISC rates calculated with the standard model (black symbols) and imposing hardened vibrational or conformational frequencies (blue and black symbols, respectively). The effect on ISC rates is modest, but a hardened frequency of the conformational mode favors RISC, in good qualitative agreement with recent experimental results.⁴⁰

Fig. 4 Temperature dependence of RISC and ISC rates calculated with model parameters in Table 1 (black curve), and multiplying by a factor of 2 either the vibrational frequency ω_v (blue curve) or the a parameter (red curve).

More interesting and intriguing is the dependence of RISC and ISC rates on the energy of the |L〉 state, as defined by the parameter k. Results in Fig. 5 show a fairly complex and non-monotonic behavior, suggesting that optimizing the rates requires a fine tuning of the position of the local excited state, in a way that may strongly depend on the details of the molecular system at hand. As discussed in ESI† (Fig. S2) this behavior, and more precisely the valley observed in either RISC and ISC rates at k ∼ 2.35 eV originates from the competing contribution from the two SOC channels driven by V_SOC, mixing N and T states and W_SOC, mixing Z and L states.

Fig. 5 Dependence of the RISC and ISC rates on the energy difference between the |L〉 and |T〉 states.

2.3 Radiative rates

Fluorescence is an allowed radiative process, occurring from the relaxed excited singlet, S₁ towards vibronic states in the S₀ manifold. The probability of the fluorescence process from state i to state f is calculated again using the FGR, as follows:⁴¹where ω_fi and μ_fi are the transition frequency and dipole moment, respectively.

The above equation can be exploited in two different approaches to the fluorescence rate. In the first approach, the complete non-adiabatic diagonalization of the Hamiltonian in eqn (2) identifies i and f into vibronic states (accounting for both vibrational and conformational motions) of the S₁ and S₀ manifold, respectively. The overall fluorescence rate from state i is then obtained summing over all decay rates towards the f states in the vibronic manifold of S₀. Finally, the overall fluorescence rate k_fluo is obtained accounting for the thermal distribution of the i states in the equilibrated S₁ manifold. This approach is fairly expensive, since the full non-adiabatic Hamiltonian must be diagonalized and, to get convergence, it requires including ∼10 vibrational states and ∼600 conformational states. Of course calculations may be limited to the two electronic states in the singlet subspace, but the total basis is fairly large, with ∼1 × 10⁴ states. In an easier approach, the conformational mode is dealt with in the adiabatic approximation: neglecting the kinetic energy of the conformational motion, the vibronic Hamiltonian accounting for the coupled electronic and vibrational motions is diagonalized for fixed values of the conformational coordinate δ. The calculation implies the diagonalization of the vibronic Hamiltonian on a basis with dimension ∼20 on a grid of δ values. For each δ, only the lowest eigenstate of the S₁ vibronic manifold, the fluorescent state, is populated (the vibrational energy is very large versus thermal energy). For each δ value, eqn (7) is summed over the f vibronic states in the S₀ manifold to obtain the δ-dependent fluorescence rate, that is finally averaged over the thermal distribution on the δ-dependent energy of the fluorescent state. We explicitly verified that the two approaches to the k_fluo estimate give exactly the same result.

The temperature-dependence of the calculated fluorescence rates is reported in the bottom panel of Fig. 6: the calculated rate, of the order of 1 × 10⁶ s⁻¹, is fairly small, as expected for S₁ state with a dominant CT character. Its decrease with decreasing temperature can be rationalized since the S₁ equilibrium geometry has δ = 0, where the transition dipole moment and hence the radiative rates vanish. Finite radiative rates are expected at finite δ:^15,42 as temperature increases, states with finite δ are progressively populated, leading to a progressive increase of the fluorescence rate.

Fig. 6 Temperature dependence of calculated ISC, RISC and fluorescence rates. Black symbols refer to gas phase results (ε_opt = 1), blue, and red symbols refer to results for matrices with ε_opt = 2.6 and 3.0, respectively. RISC and ISC rates are obtained setting τ_r = 100 fs.

3 Environmental effects: the role of medium refractive index

TADF phenomena involve molecular excitations in the visible portion of the electromagnetic spectrum. The dielectric response of a generic matrix can therefore be partitioned in two contributions with distinctively different dynamics. The electronic degrees of freedom of the matrix are characterized by timescale in in UV region, faster than the molecular degrees of freedom of interest. On the opposite, the vibrational and orientational motions of the matrix are much slower. These slow motions, on the other hand, only contribute to the dielectric response of the medium in polar matrices. Here we will only address the effect of the fast dielectric response, an will defer to a subsequent publication a detailed analysis of the dielectric response in polar matrices.

The Onsager model offers the simplest description of the interaction between a molecular species and an embedding dielectric medium.⁴³ Specifically, the solute–solvent interaction is treated in the dipolar approximation, so that the molecule feels an environmental field, F, usually called the reaction field.^44–46 If we only account for the electronic (fast) component of the dielectric response, the Hamiltonian that describes the dye in the matrix is:

where Ĥ_mol is the molecular Hamiltonian in eqn (2), F is the reaction field in energy units, amounting to the product of the electric field times μ₀, the dipole moment associated withe the zwitterionic Z and T diabatic states, T_el is the kinetic energy associated with the reaction field, and the last term is the corresponding potential energy, where ε_el measures the matrix polarizability. In the hypothesis that the solute occupies a spherical cavity of radius r₀ inside the dielectric medium, an explicit expression for ε_el can be written, in terms of the matrix dielectric constant at optical frequencies, ε_opt = η² (where η is the matrix refractive index):

Since the electronic degrees of freedom of the matrix are faster than the molecular degrees of freedom of interest, the antiadiabatic approximation can be adopted to renormalize out the medium degrees of freedom, leading to a very simple and inspiring effective molecular Hamiltonian that has precisely the same expression as in eqn (2), but with a renormalized z value:^32,34

Matrices of interest for organic light-emitting devices applications have ε_opt spanning a very narrow range comprised between 2.6 and 3.0. Setting the cavity radius to the Onsager radius (increased by 0.5 Å), r₀ = 6.44 Å, the renormalized z varies between 1.57 eV and 1.55 eV.

For the sake of comparison, Table S1 (ESI†) lists vertical transition energies (at the equilibrium ground state geometry) calculated in gas-phase and for ε_opt = 2.6 and 3 eV, in our effective model, and via TDF-DFT treating the solvent according the three flavors of PCM available in Gaussian 16.⁴⁷ As expected,³² linear response (LR) fails to reproduce the stabilization of CT states, either singlet or triplets: LR corrections are in fact proportional to transition densities, that are negligible for states like S₁ and T₁ with an almost pure CT character. Indeed, the marginal increase of the transition frequencies in LR vs the gas-phase is ascribable to the stabilization of the ground state energy. The external iteration approach fails in the opposite direction, largely overestimating the dielectric effects on transition energies. The perturbative corrected linear response (cLR) approach leads to transition energies that compare very well with those obtained in our effective model. However, in cLR the excited state wavefunctions are the same as the LR wavefunctions, and are therefore only marginally different from the gas-phase wavefunctions. Accordingly, cLR will not offer reliable estimates of SOC matrix elements or of transition dipole moment elements involving S₁ or T₁ states.

The singlet–triplet gap, calculated in the effective model as the energy difference between the lowest non-adiabatic eigentates in the singlet and triplet subspaces, decreases from the gas phase value, 0.0513 eV, to 0.0260 eV and 0.0238 eV in a dielectric medium with ε_el = 2.6 and 3, respectively. Comparing these values, obtained from non-adiabatic calculation and hence accounting for zero-point corrections, with the adiabatic gaps in TD-DFT is of little use. The comparison with experimental data is delicate as well, due to the large experimental uncertainties. Specifically, estimates for DMAC–TRZ in ref. 22 are obtained from the energy difference between the onset of fluorescence and phosphorescence bands. The very broad structure of both bands makes the estimate uncertain. Even worst, fluorescence spectra are measured at room temperature, while phosphorescence spectra are collected at 20 K.²² The dielectric response (both the refractive index and the dielectric constant) of any material show a large dependence on temperature, and, in turn, this variation largely affects the position of the singlet and triplet states.

Fig. 6 shows the temperature-dependent rates (ISC, RISC, radiative) calculated in the gas phase (z = 1.72 eV) and for the two limiting values of the effective z = 1.57 eV and 1.55 eV. Radiative rates are marginally affected by the matrix polarizability, while both RISC and ISC rates decrease considerably when going from gas phase to the matrix, suggesting that relying on gas-phase computational results may affect the quality of the analysis. On the other hand, the variability of ε_opt in common matrix is very limited, so that addressing the effects of the matrix polarizability from the experimental perspective is difficult.

RISC and ISC rates are calculated based on ORCA results for SOC matrix elements.^48,49 Out of the three spin states ORCA gives finite SOC for only one state, so that a single channel is available for either ISC or RISC. Calculated rates refer to this channel. However, in the hypothesis that the population of the three triplet states is instantaneously equilibrated, the effective RISC rate, to be compared with experiment, must be divided by a factor of 3.⁵⁰ Accordingly, we estimate the effective RISC rate for DMAC–TRZ in a non-polar matrix in the 7 × 10⁴−2 × 10⁵ s⁻¹ range, in line with experimental results falling in the range 2–5 × 10⁵ s⁻¹ in non polar matrices (Zeonex and UGH) as well as in toluene solution.²² ISC rates are estimated of the order of 2 × 10⁶ s⁻¹ approximately one order of magnitude smaller than experimental results (∼2 × 10⁷ s⁻¹). Calculated radiative rates ∼2 × 10⁶ s⁻¹ are again roughly one order of magnitude smaller than experimental fluorescence rates, which are however also affected by the non-radiative decay. In any case, we feel that the overall agreement is very good considering that we are exploiting the same model that was previously defined and parametrized against TD-DFT results and validated against experiment.³⁴

4 Discussion and conclusions

The Marcus equation is often adopted to calculate RISC rates.^12,16,51 However, attention must be payed to the basic approximation underlying the Marcus model. Specifically, the original Marcus model^52,53 or its generalization to include high frequency molecular vibrations (Marcus–Levich–Jortner, MLJ model^54,55), apply to the calculation of transfer rates between two diabatic states, whose energy depends on one or more vibrational and/or conformational and/or solvation coordinates. Relevant potential energy surfaces are usually set as harmonic, with the same frequency but displaced minima (spin-boson model with linear coupling). But the most stringent approximation is that the matrix element that mixes the two diabatic states is constant.⁵³ However, when applied to RISC rates calculations, the Marcus model is typically exploited with reference to adiabatic singlet and triplet states as obtained from quantum chemical calculations. In the adiabatic picture the state-to-state FGR, at the heart of the Marcus model, requires the calculation of the following matrix elements:

〈ψ_i(r,Q)χ_i,v(Q)|Ĥ^SOC|ψ_f(r,Q)χ_f,u(Q)〉 = 〈χ_i,v(Q)|H^SOC_if(Q)|χ_f,u(Q)〉 (11)

where ψ_i/f(r,Q) is the electronic wavefunction relevant to the initial/final state and χ_i/f,v/u is the v/uth vibrational wavefunction relevant to the i/f electronic manifold. The electronic wavefunction describes the motion of electronic coordinates r, and parametrically depends on the nuclear coordinates Q. The vibrational wavefunctions only depend on Q. In the second term of the above equation an integral over the electronic coordinates allows to define the relevant matrix element of the interaction Hamiltonian:

H^SOC_if(Q) = 〈ψ_i(r,Q)|H^SOC|ψ_f(r,Q)〉_r (12)

where 〈…〉_r stands for the integral on the electronic coordinates. If H^SOC_if(Q) is weakly dependent on Q it can be Taylor-expanded around the equilibrium geometry. The zeroth order term (the Condon term) brings the problem back to the Q-independent interactions characteristic of the Marcus model, so that the Marcus or MLJ equations can be safely applied. However if the Q dependence of H^SOC_if(Q) cannot be disregarded, the factorization of the electronic and vibrational problems is more delicate and leads to additional terms (Herzberg–Teller, etc.) well beyond the Marcus model.

In TADF systems, and specifically in DMAC–TRZ (see Fig. S1, ESI†), H^SOC_if(Q) shows a quite impressive Q-dependence as it vanishes at δ = 0, where the singlet–triplet gap closes, and becomes sizeable at intermediate angles and specifically at the triplet equilibrium geometry. Applying the Marcus or MLJ expressions in these conditions represents a strong approximation, explaining the observation of wildly different results for RISC rates estimates. Indeed, fixing the H^SOC_if(Q) to the value relevant to δ = 0 (or in any case to the value obtained where the singlet–triplet gap vanishes) leads to vanishing RISC rates.⁵¹ On the opposite, if the SOC is set to the value relevant to the equilibrium geometry for the triplet state as in ref. 12 and 15 sizable RISC rates are obtained, that for the specific case of DMAC–TRZ are roughly one order of magnitude larger than in our non-adiabatic calculation (Fig. S4, ESI†). Since, as explained above, the Marcus factorization does not apply when the Q-dependence of the interaction matrix element is non-negligible, augmenting the Marcus (or MLJ) expression just introducing a Q-dependent interaction¹⁸ is formally incorrect and leads to uncontrolled results.

The importance of properly accounting for spin-vibronic terms, well beyond the Marcus or MLJ models, has been very clearly expressed by Penfold et al.¹⁴ In their approach, using a diabatic basis, with Q independent SOC among basis states, they regain the Q-dependence of the adiabatic SOC matrix elements accounting for the non-adiabatic mixing of different diabatic states, getting closed formulas thanks to a perturbative treatment. In our essential state model³⁴ the mixing between the CT and the LE triplets, is parametrized to reproduce the adiabatic Q-dependent SOC interaction and singlet–triplet gap, as obtained from TD-DFT. Then RISC and ISC rates are obtained via a direct non-adiabatic diagonalization of the adopted effective Hamiltonian in an approach that fully accounts for non-adiabatic phenomena and the anharmonicity of the system, without the need to rely on perturbative expansions for the spin-vibronic matrix elements. Moreover, in the Penfold model a single channel for SOC is considered (corresponding to the |Z〉–|L〉 coupling in our model), so that only the mixing, driven by vibronic effects, of the |L〉 and |T〉 into T₁ states is responsible for a sizable 〈S₁|H_SOC|T₁〉 interaction. In our model, we also account for the (vibronically induced) mixing of |N〉 and |Z〉 states into the S₁ state. This mixing on one side explains the observation of a finite absorption intensity of the CT transition, on the other side, it opens a second SOC channel, related to the SOC interaction between the diabatic |N〉 and |T〉 states.

The effective model developed for DMAC–TRZ in ref. 34 opens the way to a nominally exact calculation of RISC, ISC and radiative rates fully accounting for anharmonicity and non-adiabaticity. It also allows to investigate the effects of a fine-tuning of model parameters, showing, e.g., that a softening of the conformational modes has detrimental effects on TADF. More importantly, we are in the position to address environmental effects on TADF systems. TADF is governed by a subtle interplay of singlet and triplets, LE and CT states: environmental effects, and specifically the medium polarity and polarizability have different effects on the different states affecting the singlet–triplet gap as well as the mixing between the LE and CT triplets with large effects on the SOC. A rationalization of these highly non trivial effects is however difficult as several protocols are presently available leading to contrasting results, as recently well illustrated by Mewes.¹⁵ The origin of the contrasting protocols for environmental effects on TADF (and more generally on the definition of spectral properties of molecules in condensed media) can be traced back to the improper treatment of the medium polarizability. Specifically, either continuum solvation models (including PCM, COSMO, etc.) or atomistic solvation models treat the environmental polarizability in an adiabatic approach, solving the molecular Hamiltonian for fixed environmental charges.^32,56 This approximation is well suited for polar solvation, related to very slow degrees of freedom, but fails on the environmental polarizability, related to the very fast electronic motion in the medium. This leads to a variety of solvation models, none of which is able to properly address the effect of the medium polarizability. The problem is overcome here adopting an antiadiabatic approach to address the role of the medium polarizability.³²

Here we only address the role of the matrix polarizability, as measured by the refractive index at the optical frequencies. We postpone to another publication the discussion of polarity effects that require a careful consideration of the relative rates of molecular processes (RISC, ISC, radiative and non-radiative decays) and the medium (solvent, matrix, etc.) relaxation. The marginal variability of the refractive index in common organic solvents and matrices suggests that it is impossible to optimize the TADF behavior of a dye acting on the medium refractive index. However it is extremely important to recognize that when going from gas phase to an organic medium, the large variation of the refractive index has enormous effects on the energies of the excited state and more generally on TADF. It is therefore of little use to compare results obtained from gas-phase calculations with experimental data in condensed phases. As already recognized by Mewes, indeed, the RISC rates may change by orders of magnitude when going from gas phase to a non-polar organic medium. Quite interestingly, for DMAC–TRZ, RISC and ISC rates decrease by approximately one order of magnitude when going from gas phase to the organic medium, even if the singlet triplet gap decreases. The reason for less effective RISC can be traced back to the stabilisation in the medium of the CT states, that leads to a decoupling of the CT and LE triplet and hence to a decrease of the SOC. Quite interestingly, the opposite behavior is observed by Mewes, with RISC rates increasing by a few order of magnitude from gas phase to condensed media. Apart from the use of a different solvation model, we believe that this contrasting result is indeed related to the different nature of excited states in the molecules investigated by Mewes with respect to DMAC–TRZ. In Mewes molecules, in fact, the LE triplet lays in the gas phase lower in energy than the CT triplet, but the order is reversed in non-polar matrices, leading to a large amplification of RISC.

In conclusion, we presented an original approach to the calculation of ISC and RISC rates in TADF molecules. The approach relies on the definition of a reliable essential state model for the dye, based on first-principles calculations. The model accounts for few electronic diabatic states, and for effective vibrational and conformational coordinates. Not relying on further approximations, the model fully addresses anharmonic and non-adiabatic phenomena. Quite interestingly, the model lends itself quite naturally to describe environmental effects: here we discuss the role of the environmental polarizability, the more subtle polar effects are deferred to a subsequent publication.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This project received funding from the European Union Horizon 2020 Research and Innovation Programme under Grant Agreement No. 812872 (TADFlife), and benefited from the equipment and support of the COMP-HUB Initiative, funded by the “Departments of Excellence” program of the Italian Ministry for Education, University and Research (MIUR, 2018–2022), and benefited from the equipment and the support from the HPC (High Performance Computing) facility of the University of Parma, Italy.

References

1. C. A. Parker and C. G. Hatchard, Trans. Faraday Soc., 1961, 57, 1894–1904.
2. A. Endo, K. Sato, K. Yoshimura, T. Kai, A. Kawada, H. Miyazaki and C. Adachi, Appl. Phys. Lett., 2011, 98, 083302.
3. H. Uoyama, K. Goushi, K. Shizu, H. Nomura and C. Adachi, Nature, 2012, 492, 234–238.
4. M. K. Etherington, J. Gibson, H. F. Higginbotham, T. J. Penfold and A. P. Monkman, Nat. Commun., 2016, 7, 13680.
5. M. K. Etherington, F. Franchello, J. Gibson, T. Northey, J. Santos, J. S. Ward, H. F. Higginbotham, P. Data, A. Kurowska, P. L. D. Santos, D. R. Graves, A. S. Batsanov, F. B. Dias, M. R. Bryce, T. J. Penfold and A. P. Monkman, Nat. Commun., 2017, 8, 14987.
6. C. S. Oh, D. d. S. Pereira, S. H. Han, H.-J. Park, H. F. Higginbotham, A. P. Monkman and J. Y. Lee, ACS Appl. Mater. Interfaces, 2018, 10, 35420–35429.
7. F. B. Dias, J. Santos, D. R. Graves, P. Data, R. S. Nobuyasu, M. A. Fox, A. S. Batsanov, T. Palmeira, M. N. Berberan-Santos, M. R. Bryce and A. P. Monkman, Adv. Sci., 2016, 3, 1600080.
8. T.-T. Huang and E. Y. Li, Org. Electron., 2016, 39, 311–317.
9. X.-K. Chen, S.-F. Zhang, J.-X. Fan and A.-M. Ren, J. Phys. Chem. C, 2015, 119, 9728–9733.
10. J. Gibson, A. P. Monkman and T. J. Penfold, ChemPhysChem, 2016, 17, 2956–2961.
11. I. Lyskov and C. M. Marian, J. Phys. Chem. C, 2017, 121, 21145–21153.
12. P. K. Samanta, D. Kim, V. Coropceanu and J.-L. Brédas, J. Am. Chem. Soc., 2017, 139, 4042–4051.
13. T. J. Penfold, E. Gindensperger, C. Daniel and C. M. Marian, Chem. Rev., 2018, 118, 6975–7025.
14. T. J. Penfold, E. Gindensperger, C. Daniel and C. M. Marian, Chem. Rev., 2018, 118, 6975–7025.
15. J.-M. Mewes, Phys. Chem. Chem. Phys., 2018, 20, 12454–12469.
16. Y. Olivier, J.-C. Sancho-Garcia, L. Muccioli, G. DAvino and D. Beljonne, J. Phys. Chem. Lett., 2018, 9, 6149–6163.
17. P. de Silva, C. A. Kim, T. Zhu and T. Van Voorhis, Chem. Mater., 2019, 31, 6995–7006.
18. L. E. de Sousa and P. de Silva, J. Chem. Theory Comput., 2021, 17, 5816–5824.
19. B. C. Garain, P. K. Samanta and S. K. Pati, J. Phys. Chem. A, 2021, 125, 6674–6680.
20. G. Méhes, K. Goushi, W. J. Potscavage and C. Adachi, Org. Electron., 2014, 15, 2027–2037.
21. P. L. dos Santos, J. S. Ward, M. R. Bryce and A. P. Monkman, J. Phys. Chem. Lett., 2016, 7, 3341–3346.
22. K. Stavrou, L. G. Franca and A. P. Monkman, ACS Appl. Electron. Mater., 2020, 2, 2868–2881.
23. Y. Olivier, B. Yurash, L. Muccioli, G. D'Avino, O. Mikhnenko, J. C. Sancho-García, C. Adachi, T.-Q. Nguyen and D. Beljonne, Phys. Rev. Mater., 2017, 1, 075602.
24. H. Sun, Z. Hu, C. Zhong, X. Chen, Z. Sun and J.-L. Brédas, J. Phys. Chem. Lett., 2017, 8, 2393–2398.
25. T. Northey, J. Stacey and T. J. Penfold, J. Mater. Chem. C, 2017, 5, 11001–11009.
26. L. Kunze, A. Hansen, S. Grimme and J.-M. Mewes, J. Phys. Chem. Lett., 2021, 12, 8470–8480.
27. J. Tomasi, B. Mennucci and R. Cammi, Chem. Rev., 2005, 105, 2999–3094.
28. A. V. Marenich, C. J. Cramer, D. G. Truhlar, C. A. Guido, B. Mennucci, G. Scalmani and M. J. Frisch, Chem. Sci., 2011, 2, 2143–2161.
29. A. Klamt, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2011, 1, 699–709.
30. T. Vreven and K. Morokuma, Hybrid Methods: ONIOM(QM:MM) and QM/MM in Annual Reports in Computational Chemistry, ed. D. C. Spellmeyer, Elsevier, 2006, ch. 3, vol. 2, pp. 35–51.
31. J. Li, G. D'Avino, I. Duchemin, D. Beljonne and X. Blase, Phys. Rev. B, 2018, 97, 035108.
32. D. K. A. Phan Huu, R. Dhali, C. Pieroni, F. Di Maiolo, C. Sissa, F. Terenziani and A. Painelli, Phys. Rev. Lett., 2020, 124, 107401.
33. D. K. A. P. Huu, C. Sissa, F. Terenziani and A. Painelli, Phys. Chem. Chem. Phys., 2020, 22, 25483–25491.
34. R. Dhali, D. K. A. Phan Huu, F. Terenziani, C. Sissa and A. Painelli, J. Chem. Phys., 2021, 154, 134112.
35. R. Dhali, D. K. A. P. Huu, F. Bertocchi, C. Sissa, F. Terenziani and A. Painelli, Phys. Chem. Chem. Phys., 2021, 23, 378–387.
36. W.-L. Tsai, M.-H. Huang, W.-K. Lee, Y.-J. Hsu, K.-C. Pan, Y.-H. Huang, H.-C. Ting, M. Sarma, Y.-Y. Ho, H.-C. Hu, C.-C. Chen, M.-T. Lee, K.-T. Wong and C.-C. Wu, Chem. Commun., 2015, 51, 13662–13665.
37. M. A. El-Sayed, J. Chem. Phys., 1963, 38, 2834–2838.
38. A. Painelli, Chem. Phys. Lett., 1998, 285, 352–358.
39. F. Di Maiolo and A. Painelli, J. Chem. Theory Comput., 2018, 14, 5339–5349.
40. M. Hempe, N. A. Kukhta, A. Danos, M. A. Fox, A. S. Batsanov, A. P. Monkman and M. R. Bryce, Chem. Mater., 2021, 33, 3066–3080.
41. R. Loudon, The Quantum Theory of Light, Oxford University Press, Oxford New York, 2000.
42. Y. Olivier, M. Moral, L. Muccioli and J.-C. Sancho-García, J. Mater. Chem. C, 2017, 5, 5718–5729.
43. L. Onsager, J. Am. Chem. Soc., 1936, 58, 1486–1493.
44. E. G. McRae, J. Phys. Chem., 1957, 61, 1128–1129.
45. S. Di Bella, T. J. Marks and M. A. Ratner, J. Am. Chem. Soc., 1994, 116, 4440–4445.
46. A. Painelli, Chem. Phys., 1999, 245, 185–197.
47. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16, Revision B.01, Gaussian, Inc., Wallingford CT, 2016.
48. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 73–78.
49. F. Neese, J. Chem. Phys., 2005, 122, 034107.
50. C. M. Marian, J. Föller, M. Kleinschmidt and M. Etinski, in Intersystem Crossing Processes in TADF Emitters, John Wiley & Sons, Ltd, 2018, ch. 8, pp. 257–296.
51. N. Aizawa, Y. Harabuchi, S. Maeda and Y.-J. Pu, Nat. Commun., 2020, 11, 3909.
52. R. A. Marcus, Angew. Chem., Int. Ed. Engl., 1993, 32, 1111–1121.
53. A. Nitzan, Chemical dynamics in condensed phases: relaxation, transfer, and reactions in condensed molecular systems, Oxford University Press, Oxford, 2013.
54. N. R. Kestner, J. Logan and J. Jortner, J. Phys. Chem., 1974, 78, 2148–2166.
55. J.-L. Brédas, D. Beljonne, V. Coropceanu and J. Cornil, Chem. Rev., 2004, 104, 4971–5004.
56. R. Dhali, D. K. A. P. Huu, F. Terenziani, C. Sissa and A. Painelli, J. Chem. Phys., 2021, 154, 134112.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d1tc05296a

---