Thermally Activated Delayed Fluorescence (TADF)¶

Electronic Structure of OLED Materials¶

OLED devices typically consist of a thin layer of a host material which is doped with an organic (or metalorganic) emitter material. Applying an electric voltage across such semiconductor layers results in free electrons and electron holes, so-called charge carriers, travelling through the material. When these charge carriers recombine, strongly bound electron-hole pairs, so-called excitons are formed, which are entirely localized on a single emitter molecule and thereby correspond to electronic excitations.

As free electrons (and holes) are equally likely to be created with spin up or spin down configurations, 75% of the formed excitons emerge in a triplet state (total spin \(S\) = 1) while 25% are found in a singlet state (\(S\) = 0).

Fig. 27 Mechanism of first generation OLED emitters involving a fluorescent decay path with a 25% maximum theoretical quantum yield.

This spin statistics negatively affects the quantum efficiency of an OLED material, hence on the proportion of electron-hole pairs that successfully undergo a radiative decay in to the ground state \(S_{\mathrm{0}}\) and emit a photon. A direct de-excitation from \(T_{\mathrm{1}}\) to \(S_{\mathrm{0}}\) is spin-forbidden in a non-relativistic theory. In the absence of any spin-orbit coupling effects in the dye, a radiative \(T_{\mathrm{1}} \rightarrow S_{\mathrm{0}}\) emission becomes very unlikely and light is only be emitted from a fluorescent \(S_{\mathrm{1}} \rightarrow S_{\mathrm{0}}\) transition.

To address this problem a second generation of OLED emitters was developed involving OLED dyes containing heavy element atoms. Electrons in the valence shells of such atoms (e.g. Iridium) are subject to significant spin-orbit coupling (SOC). SOC effects in turn, permit \(T_{\mathrm{1}} \rightarrow S_{\mathrm{0}}\) transitions and facilitate so-called intersystem-crossings (ISC) e.g. from \(S_{\mathrm{1}}\) to \(T_{\mathrm{1}}\). As a result the \(T_{\mathrm{1}}\) states become more populated. From there the system undergoes a phosphorescent \(T_{\mathrm{1}} \rightarrow S_{\mathrm{0}}\) decay and emit light.

Fig. 28 Second generation OLED dyes rely on intersystem crossing to convert all excitons into a triplet state from which phosphorescent emissions occur.

While SOC-based emitters can reach quantum efficiencies of nearly 100%, the presence of long-living triplet states can cause chemical degradations. This tendency increases with the energy difference to the ground state (and thus with the frequency of the emitted light). As a result no stable, long-lasting OLED materials to deep blue light have been found as of yet. Apart from that, SOC-based dyes involve relatively rare and expensive heavy elements which significantly increases their costs.

These disadvantages motivated the development of third generation OLED materials: Dye molecules in which the states \(S_{\mathrm{1}}\) and \(T_{\mathrm{1}}\) are close in energy and (to some degree) vibrationally coupled, can exhibit a phenomenon known as thermally activated delayed fluorescence (TADF). Using emitter molecules specifically tailored to maximize TADF, this effect can be exploited to avoid the degradation problems of SOC-based emitters. In addition, TADF occurs in purely organic dyes so that these third generation OLED emitters do not require expensive heavy metal elements.

TADF involves a reverse intersystem crossing (RISC) process, in which \(S_{\mathrm{1}}\) states are populated from \(T_{\mathrm{1}}\) levels. This is followed again by a normal fluorescent transition \(S_{\mathrm{1}} \rightarrow S_{\mathrm{0}}\) and allows for light emissions with nearly 100% quantum yields.

Fig. 29 The thermally assisted reverse intersystem crossing mechanism in third generation OLED materials converts excitons from triplet into singlet states from which a fluorescent radiation process can occur.

RISC competes with other processes such as the \(S_{\mathrm{1}} \rightarrow T_{\mathrm{1}}\) intersystem crossing. However, if non-radiative paths are negligible, \(k_{\mathrm{RISC}}\) represents the rate-determining factor for TADF.

TADF represents a rather complex process and \(k_{\mathrm{RISC}}\) has been shown to be strongly affected by various individual properties of a given dye material:

• singlet-triplet gap \(\Delta E_{\mathrm{ST}} = E(S_{\mathrm{1}}) - E(T_{\mathrm{1}})\)
• spin-orbit interactions
• vibrational couplings between different electronic states
• solvent effects

Each of these aspects will be addressed in the following.

Excited States Geometry Optimizations¶

Following Samanta et al., we use TDDFT calculations in the Tamm-Dancoff approximation (TDA) to describe the TADF process in 1,2,3,5-tetrakis(carbazol-9-yl)-4,6-dicyanobenzene (4CzIPN), a typical TADF emitter.

Fig. 30 1,2,3,5-tetrakis(carbazol-9-yl)-4,6-dicyanobenzene (4CzIPN)

Thereby we use the range-separated LCY-PBE functional, which employs a Yukawa potential to realize the switching function between long- and short-range terms. The sharpness of the Yukawa potential is controlled by the parameter \(\gamma\), for which a value of 0.22 is used in the following.

To model the TADF process of 4CzIPN properly, the geometry of the system needs to be relaxed for all three relevant electronic states, \(S_{\mathrm{0}}\), \(S_{\mathrm{1}}\), and \(T_{\mathrm{1}}\), respectively. We can do this from GUI with ADF with the following settings:

Download the 4CzIPN structure here

Start AMSinput

Load the structure with File → Import Coordinates…

Select the following in the Main panel:

Task → Geometry Optimization

XC Functional → Range Separated → LCY PBE

Click on right of the XC Functional field

Enter Gamma 0.22

Return to the Main panel

Make sure that Relativity → Scalar

Basis set → TZP

Frozen core → None

Numerical quality → Good

Click on the Model panel and select Solvation

Solvation method → COSMO

COSMO solvent → Toluene

Save the input File → Save as… under e.g. 4CzIPN_GSopt.ams

Run the calculation File → Run

A few more additional options are needed to define the excited states geometry optimizations. In the case of the \(S_{\mathrm{1}}\)-optimization:

In the window with the previously created 4CzIPN_GSopt.ams file still open:

Click in the panel bar on Properties and select Excitations (UV/Vis), CD

Type of excitations → SingletOnly

TDA → Yes

Enter the Number of excitations 2

Tick the Calculate → Charge transfer descriptors box

Click on Properties panel and select Excited State Geometry

Enter at the Excitation 1A

Spin type → Singlet

Save the input File → Save as… under e.g. 4CzIPN_S1opt.ams

Run the calculation File → Run

Analogously for the \(T_{\mathrm{1}}\)-state:

Click on the Properties panel and select Excitations (UV/Vis), CD

Type of excitations → TripletOnly

TDA → Yes

Enter the Number of excitations 2

Tick the Calculate → Charge transfer descriptors box

Click on the Properties panel and select Excited State Geometry

Excitation enter 1A

Spin type → Triplet

Save the input File → Save as… under e.g. 4CzIPN_T1opt.ams

Run the calculation File → Run

We are now ready to run all three optimization jobs. Note that, for a system like 4CzIPN an excited states geometry optimization with RSH can take multiple hours. For your convenience, the relaxed structures for the \(S_{\mathrm{0}}\), \(S_{\mathrm{1}}\), and \(T_{\mathrm{1}}\) states, respectively can be downloaded here:

After these optimizations are concluded, the respective minimum energies of the states can be extracted at the end of the logfile. For the ground state,

<Nov12-2020> <03:56:22>   Bond Energy         -33.76473746 a.u.

For excited states the end of the logfile looks like this (the T1 state in this case):

<Nov13-2020> <04:18:25>  Total excited state energy:                 -33.67697839 Hartree

From there we can extract the following energy differences (neglecting vibrational contributions for the moment):

The result is somewhat different from the results (\(\Delta E_{\mathrm{ST}}\) = 0.01 eV) of the original paper due to the different electronic structure description used therein. Compared to the experimental value of \(\Delta E_{\mathrm{ST}}\) = 0.1 eV, the difference is however smaller and we will use our computed results of 0.16 eV in the subsequent calculations of this tutorial.

The output file provides the charge transfer descriptors according to the prescription of Plasser and Lischka, which are obtained from transition density matrices.

For the \(T_{\mathrm{1}}\) state we find a charge transfer (CT) character of 82% with an average electron-electron hole distance (R_HE) of 4.6 Å, while in the case of \(S_{\mathrm{1}}\) the CT character amounts to 91% with an average of 5.0 Å between the excited electron and electron hole. Note, that such high CT characters and electron-hole distances are a direct consequence of the design principles of typical OLED emitters featuring well-separated electron donor and acceptor groups. In consequence, these excited states also exhibit a relatively small exchange energy, which contributes most to the small \(\Delta E_{\mathrm{ST}}\) value computed above from TDDFT excitation energies.

Vertical Absorption¶

Following Samanta et al. and to analyze the system further, we compute vertical excitations at the ground state geometry to obtain the absorption spectrum of 4CzIPN. Note, that these results do not enter the calculation of the RISC rate but rather serve as an additional assessment of the accuracy of the description employed here. Regarding most aspects this is done analogously to the previous calculations. The main differences are that we now treat the excitations in a spin-orbit description which allows for pure spin-states (i.e. no spin contamination). Furthermore, the response of the solvation model is ignored as the solvent molecules cannot rearrange within the short timespan of a vertical excitation; a better approximation of this interaction requires the optical part of the solvents dielectric constant (see the corresponding section in the ADF Manual).

Open a new instance of AMSinput

Load the previously optimized GS structure with File → Import Coordinates…

In the Main panel:

Task → Single Point

XC Functional → Range Separated → LCY PBE

Click on right of the XC Functional field

Enter Gamma 0.22

Return to the Main panel

Make sure that Relativity → Scalar

Basis set → TZP

Frozen core → None

Numerical quality → Good

Click on the Model panel and select Solvation

Solvation method → COSMO

COSMO solvent → Toluene

Click on the Details panel and select COSMO

Untick Include response in TDDFT

Click on the Properties panel and select Excitations (UV/Vis), CD

Type of excitations → Spin-Orbit (Perturbative)

TDA → Yes

Enter the Number of excitations 5

Tick the Calculate → Charge transfer descriptors box

Save the input File → Save as… under e.g. 4CzIPN_GS_VA.ams

Run the calculation File → Run

After completing this calculation, we open the results with AMSoutput and search for the results with Response Properties → All Spin-Orbit Coupling Excitation Energies:

Recalling the spin-forbidden transitions appear with near vanishing oscillator strengths in an SOC treatment, we note the following from these results:

• Excitation no. 8 is the lowest with a relevant oscillator strength \(f\), it corresponds to \(S_{\mathrm{0}} \rightarrow S_{\mathrm{1}}\).
• The \(S_{\mathrm{0}} \rightarrow S_{\mathrm{1}}\) excitation energy of 2.77 eV compares well with the experiment value, 2.85 eV.
• At the GS geometry the \(T_{\mathrm{1}}\) and the \(T_{\mathrm{2}}\) state lie below \(S_{\mathrm{1}}\).

Calculating Spin-Orbit Couplings¶

The reverse intersystem crossing process itself occurs on a too short timescale for the systems geometry to readapt during the \(T_{\mathrm{1}}\rightarrow S_{\mathrm{1}}\) transition. The factor \(|V_{\mathrm{SOC}}|^{2}\) in the above expression for \(k_{\mathrm{RISC}}\) is therefore computed at the optimized geometry of the \(T_{\mathrm{1}}\) state.

Open a new instance of AMSinput

Load the previously optimized T1 structure with File → Import Coordinates…

In the Main panel:

Task → Single Point

XC-Functional → Range Separated → LCY PBE

Click on right of the XC-Functional field

Enter Gamma 0.22

Return to the Main panel

Make sure that Relativity → Scalar

Basis set → TZP

Frozen core → None

Numerical quality → Good

Click on the Model panel and select Solvation

Solvation method → COSMO

COSMO solvent → Toluene

Click on the Details panel and select COSMO

Make sure that the Include response in TDDFT box is ticked

Click on the Properties panel and select Excitations (UV/Vis), CD

Type of excitations → Spin-Orbit (Perturbative)

TDA → Yes

Enter the Number of excitations 5

Tick the Calculate → Charge transfer descriptors box

Save the input File → Save as… under e.g. 4CzIPN_SOCME.ams

After the calculation is completed, open the output file with AMSoutput. In the AMSoutput the SOCME results can be found with Response Properties → Spin-orbit couplings:

The above root mean square coupling element \(|\langle S_{\mathrm{1}}|\hat{H}_{\mathrm{SOC}}|T_{\mathrm{1}}\rangle|\) is then divided by three to represent the average over the three triplet states (see above). From the \(|\langle S_{\mathrm{1}}|\hat{H}_{\mathrm{SOC}}|T_{\mathrm{1}}\rangle|\) value of 0.4 cm^-1 we then obtain \(|V_{\mathrm{SOC}}|^{2} =\) (0.4 cm^-1)² / 3 = 8.2\(\cdot\)10^-10 eV².

Marcus Theory¶

Classical Marcus theory (CMT) provides a simple way to describe this thermokinetic barrier. The potential energy surfaces of both states around their respective minimum structures are thereby described in terms of parabolic functions.

Fig. 31 Potential energy surfaces of the \(S_{\mathrm{1}}\) and \(T_{\mathrm{1}}\) states as quadratic functions centered at \(x_{\mathrm{S1}}\) and \(x_{\mathrm{T1}}\), respectively. In Marcus theory the difference between \(E_{\mathrm{S1/S1}}\) and \(E_{\mathrm{S1/T1}}\) is used to estimate the thermokinetic barrier of the \(T_{\mathrm{1}}\rightarrow S_{\mathrm{1}}\) transition.

Assuming the same curvature of both states, i.e. \(E(T_{1}) = a\cdot x^{2}\) and \(E(S_{1}) = a\cdot (x_{} + x_{\mathrm{S1}} - x_{\mathrm{T1}})^{2} + \Delta E_{\mathrm{ST}}\), one finds the for the \(T_{\mathrm{1}}\rightarrow S_{\mathrm{1}}\) barrier height as the intersection of both parabolas:

Using this barrier height in an Arrhenius-type expression yields then the RISC rate:

The reorganization energy \(\lambda\) can thereby be considered as the energy required to bring the system in the \(S_{\mathrm{1}}\) to the minimum structure of the \(T_{\mathrm{1}}\) state. Designating the energies of the \(S_{\mathrm{1}}\) state at \(x_{\mathrm{S1}}\) and \(x_{\mathrm{T1}}\) as \(E_{\mathrm{S1/S1}}\) and \(E_{\mathrm{S1/T1}}\), respectively, one can use

to calculate \(k_{\mathrm{RISC}}\). We can thereby set up the necessary calculations analogously to the vertical absorption calculations:

Open a new instance of AMSinput

Load the previously optimized S1 structures with File → Import Coordinates…

In the Main panel:

Task → Single Point

XC-Functional → Range Separated → LCY PBE

Click on the right of the XC-Functional field

Enter Gamma: 0.22

Return to the Main panel

Make sure that Relativity → Scalar

Basis set → TZP

Frozen core → None

Numerical quality → Good

Click on the Model panel and select Solvation

Solvation method → COSMO

COSMO solvent → Toluene

Click on the Details panel and select COSMO

Untick Include response in TDDFT

Click on the Properties panel and select Excitations (UV/Vis), CD

Type of excitations → Spin-Orbit (Perturbative)

TDA → Yes

Enter the Number of excitations 5

Tick the Calculate → Charge transfer descriptors box

Save the input File → Save as… under e.g. 4CzIPN_S1_S1.ams

Analogously for the S1/T1 calculation we set:

Load the previously optimized T1 structures with File → Import Coordinates…

Repeat the above steps

Save the input File → Save as… under e.g. 4CzIPN_S1_T1.ams

After running the calculations the required energies can be extracted again with AMSoutput. Note, that the excited states energies of \(E_{\mathrm{S1/S1}}\) and \(E_{\mathrm{T1/S1}}\) are obtained as sum of the ground state energy and the respective excitation energy. These can be retrieved in AMSoutput as Total Bonding Energy under Properties → Bonding Energy Decomposition and as the first excitation with a significant oscillator strength under Response Properties → All Spin-Orbit Coupling Excitation Energies, respectively. So, \(E_{\mathrm{S1/S1}}\) and \(E_{\mathrm{S1/T1}}\) can be calculated by addition of its first singlet excitation energy to its bond energy.

For 4CzIPN, this results in a value of \(\lambda\) = \(E_{\mathrm{S1/T1}} - E_{\mathrm{S1/S1}}\) = -916.153 eV - -916.230 eV = 0.076 eV. This result can now be used along with the previously obtained values \(\Delta E_{\mathrm{ST}}\) = 0.16 eV and \(|V_{\mathrm{SOC}}|^{2} =\) 8.2\(\cdot\)10^-10 eV² in the above expressions for \(\rho_{\mathrm{FCWD}}^{\mathrm{CMT}}\) and \(k_{\mathrm{RISC}}\) to yield at a temperature of \(T\) = 298 K:

and thus a reverse intersystem crossing rate of \(k_{\mathrm{RISC}}\) = 5.24\(\cdot\)10⁴ / s. This result underestimates the crossing rate compared to the experimental reference of \(k_{\mathrm{RISC}}\) = 2.7\(\cdot\)10⁶ / s. The reason for this might be due to small discrepancies in either \(\Delta E_{\mathrm{ST}}\) or \(\lambda\) which get amplified by the exponential term in \(\rho_{\mathrm{FCWD}}^{\mathrm{CMT}}\).

Compared to that, Samanta et al. obtained higher rates for 4CzIPN (4.9\(\cdot\)10⁶ – 1.8\(\cdot\)10⁷/ s). This difference can be rationalized by their usage of different DFT approximations, the significantly smaller value of \(\Delta E_{\mathrm{ST}}\) obtained thereby, and a different handling of the spin-orbit couplings. We also note that Samanta et al. argued that for emitters like 4CzIPN \(\lambda\) can be replaced with an effective value between 0.1 and 0.2 eV. When using \(\lambda\) = 0.1 eV instead of the previously calculated reorganization energy, the RISC rate rises to \(k_{\mathrm{RISC}}\) = 8.8\(\cdot\)10⁵ / s.

Franck-Condon Principle and Marcus-Levich-Jortner Theory¶

While the classical Marcus theory often provides acceptably accurate results, it relies on some strong assumptions and the real situation turns out more intricate: First, the potential energy surfaces around both minima are described by the entirety of vibrational normal modes. Each normal mode defines a separate harmonic oscillator system which, in a quantum mechanical description, leads to a separate set of vibrational eigenfunctions (quanta). The Franck-Condon principle then states that the \(T_{\mathrm{1}}\rightarrow S_{\mathrm{1}}\) transition requires a significant overlap between nuclear wavefunctions of both electronic states.

Fig. 32 The Franck-Condon principle describes \(T_{\mathrm{1}}\rightarrow S_{\mathrm{1}}\) transition rate in terms of the overlap between the nuclear vibronic wavefunctions of the initial and final state, respectively.

One can distinguish between two limiting cases: the quanta low frequency vibrations from shallow modes and from rearranging solvent molecules (see below) typically form an essentially continuous spectrum. Pairs of overlapping quanta from low frequency modes will therefore occur at energies close to the intersection of the two parabolic potential curves, which leads to results equivalent to those of classical Marcus theory. Within classical Marcus theory the effect of these modes can be described in terms of a Marcus reorganization energy \(\lambda_{\mathrm{M}}\). If the transition is dominated by high frequency modes, nuclear tunneling effects and the overlap of individual pairs of quanta are prevalent. The effect of such modes consists mainly in a renormalization of the coupling parameters.

An extension of CMT, Marcus-Levich-Jortner theory, includes both of these limiting cases. The high frequency modes can thereby be projected onto the transition path which summarizes them into a single effective mode with frequency \(\omega_{eff}\) and gives rise to the following expression:

Thereby, \(S_{eff}\) denotes the Huang-Rhys factor of this effective mode, which is a measure for the strength of the electron-phonon coupling. \(S_{eff}\) and the corresponding frequency \(\omega_{eff}\) are obtained from the Huang-Rhys factors and frequencies of the individual high frequency normal modes.

Whereas modes with \(\omega_{i} >\) 1000 cm^-1 are typically included in the summation.

Samanta et al. computed \(S_{i}\) factors from mode resolved reorganization energies \(\lambda_{i} = k_{i}\Delta Q_{i}^{2}/2\), which in turn results from the frequencies \(k_{i} = \omega_{i}^{2}\) and \(\Delta Q_{i}\), the normal mode projection onto the transition path.

For a large emitter molecule like 4CzIPN however, the calculation of the full Hessian matrix turns out to be computationally too expensive.

Effective Modes and Huang-Rhys Factors from DFTB and FCF¶

In some cases, the DFTB method can be used as a more approximate but efficient alternative to compute the normal modes. The FCF program of the Amsterdam Modeling Suite can then be used to calculate the \(S_{i}\) factors which then lead to the \(S_{eff}\) and \(\omega_{eff}\) values. The following procedure illustrates how this can be achieved with the Amsterdam Modeling Suite.

Open a new instance of AMSinput

Load the previously optimized \(S_{\mathrm{1}}\) structure via File → Import Coordinates

In the panel bar click →

Task → Geometry Optimization

Select the Frequencies box

Model → SCC-DFTB

Parameter directory → DFTB.org/mio-1-1 (or similar)

Properties → Excitations (UV/Vis)

Type of excitations → Singlet

Enter the Number of excitations 5

Calculate the excited state gradient for Excitation number 1

Save the input File → Save as…

Analogously, for the \(T_{\mathrm{1}}\) state:

Open a new instance of AMSinput

Load the previously optimized \(T_{\mathrm{1}}\) structure via File → Import Coordinates

…

Properties → Excitations (UV/Vis)

Type of excitations → Triplet

Enter the Number of excitations 5

Calculate the excited state gradient for Excitation number 1

Save the input File → Save as…

After completing both calculations, the binary result files dftb.rkf containing the normal modes are found in the respective directories <jobname>.results. Create a new directory and copy these rkf files naming them e.g. S1_optFreqDFTB.rkf and T1_optFreqDFTB.rkf. Based on these results the program-tool fcf can be invoked to compute \(S_{eff}\) and \(\omega_{eff}\) values using the following bash script:

$AMSBIN/fcf << eor
STATES T1_optFreqDFTB.rkf S1_optFreqDFTB.rkf
QUANTA 0 0
MODES 144 276
TRANSLATE
ROTATE
eor

Note, that we only couple the quanta 0 of each states as higher vibronic quanta are expected to only have insignificant contributions and because the number of pairs of vibronic wavefunctions to couple quickly becomes overwhelming. We also select the range of modes (Modes 144 276) to include all vibrations above 1000 cm-1. Running this script yields the frequencies and Electron-Phonon couplings for both states as shown in the following example:

=====================
F I R S T   S T A T E
=====================

  Frequency         Displacement  Electron-Phonon coupling
    cm^-1            a0 amu^1/2        (dimensionless)
 ----------          ----------          ----------
 987.027181            0.000001            0.000002
 1041.665759           0.000169            0.000352
 1077.60317           -0.031256            0.06612
 ...
 3075.033683           0.046442            0.165962
 3085.939917           0.014746            0.052789
 3086.713042          -0.000005            0.000017

 ...

We are interested in the \(T_{\mathrm{1}}\rightarrow S_{\mathrm{1}}\) transition and are therefore using the results from the first state. The mode-resolved Huang-Rhys factors \(S_{i}\) are obtained as the square of the above Electron-Phonon coupling parameters. We then sum these squares of all modes with \(\omega_{i} >\) 1000 cm^-1 to get the Huang-Rhys factor \(S_{eff} = \sum\limits_{i}S_{i}\). The frequencies corresponding to the Huang-Rhys factor, \(\omega_{eff}\), can be calculated by using these values: \(\omega_{eff} = \sum_{i}\frac{\omega_{i}S_{i}}{S_{eff}}\). This results in \(S_{eff}\) = 0.792 and \(\omega_{eff}\) = 1500 cm^-1.

These quantities are then used in the above expression for \(\rho_{\mathrm{FCWD}}^{\mathrm{CMT}}\) (with \(\lambda_{\mathrm{M}}\) = \(\lambda\) = 0.074 eV) whereas the summation is stopped at \(n\) = 5. With \(|V_{\mathrm{SOC}}|^{2} =\) 8.2\(\cdot\)10^-10 eV² this results in a reverse intersystem crossing rate that is essentially identical to the result from classical Marcus theory of \(k_{\mathrm{RISC}}\) = 5.24\(\cdot\)10⁴ / s. This supports the finding of Samanta et al. that quantum tunneling effects are negligible in the reverse intersystem crossing process of 4CzIPN and that a classical treatment suffices for its description.