The Importance of Vibronic Coupling for Efficient Reverse Intersystem Crossing in Thermally Activated Delayed Fluorescence Molecules

Abstract Factors influencing the rate of reverse intersystem crossing (k rISC) in thermally activated delayed fluorescence (TADF) emitters are critical for improving the efficiency and performance of third‐generation heavy‐metal‐free organic light‐emitting diodes (OLEDs). However, present understanding of the TADF mechanism does not extend far beyond a thermal equilibrium between the lowest singlet and triplet states and consequently research has focused almost exclusively on the energy gap between these two states. Herein, we use a model spin‐vibronic Hamiltonian to reveal the crucial role of non‐Born‐Oppenheimer effects in determining k rISC. We demonstrate that vibronic (nonadiabatic) coupling between the lowest local excitation triplet (3LE) and lowest charge transfer triplet (3CT) opens the possibility for significant second‐order coupling effects and increases k rISC by about four orders of magnitude. Crucially, these simulations reveal the dynamical mechanism for highly efficient TADF and opens design routes that go beyond the Born‐Oppenheimer approximation for the future development of high‐performing systems.


S1 Linear Spin Vibronic Model Hamiltonian
The Hamiltonian used in the present work is based upon a Linear model Vibronic Coupling Hamiltonian described in ref. 1. Here the diabatic potentials, most appropriate for quantum dynamics simulations, are obtained by fixing a point at which the adiabatic and diabatic potentials are equivalent (i.e. the coupling is zero). This is chosen to be the Franck-Condon point. The Hamiltonian is then expanded up to first order as a Taylor series around this point Q 0 , using dimensionless (mass-frequency scaled) normal mode coordinates: The zeroth order term is the ground state harmonic oscillator approximation. The zeroth order coupling matrix contains the adiabatic state energies for each state at Q 0 , i.e. The Franck-Condon excitation energies. The elements of the first order linear coupling matrix (W (1) ) are written: , usually represented by the symbol λ, are present and act between the two triplet states ( 3 CT and 3 LE). An important consideration when setting up a model Hamiltonian for large polyatomic molecules is which nuclear degrees of freedom are required in the model to ensure an accurate description of the dynamics of interest. Here we used symmetry considerations to determine all of the mode for which coupling between the two triplet states could be non-zero. At the equilibrium geometry PTZ-DBTO2 has C s symmetry, and therefore the symmetry rule can be written: where Γ i and Γ j are symmetries of the two states and Γ α is the symmetry of the mode. Therefore off diagonal elements (nonadiabatic coupling elements) will only be non-zero if the product of the two states gives the symmetry of the specific mode. As the 3 CT and 3 LE states are of A" and A' symmetry, respectively only modes of A" symmetry will couple them. Subsequently all A" modes with a frequency <500 cm −1 , i.e. those which can reasonably be thermally activated, (15 in total) were calculated and fitted to determine the magnitude of non-adiabatic coupling. Only modes with significant coupling, > 50cm −1 were included. These are given in Table S1. The potentials are shown in Figure 2 in the main text and Figure S1.
A" 422.60 D-A Torsion Table S1: Mode symmetry and vibration energies (in cm −1 ) for the normal modes of PTZ-DBTO2 found to promote significant vibronic coupling and which are included in the model Hamiltonian.
The potential energy surfaces were calculated using TDDFT(M062X) [2] within the Tamm-Dancoff approximation (TDA) [3] and a def2-TZVP basis set [4] as implemented within the Gaussian quantum chemistry package [5]. While TDDFT is well documented to have difficulties in addressing charge transfer excitations, the large fraction (54%) of non local exact exchange incorporated within the M062X functional remedies, to a large extent, this problem [6]. In addition, the TDA is critical for avoiding problems associated with the triplet instability [7]. The excited state energies for the relevant excited states at the Franck-Condon geometry are shown in Table S2 and are in good agreement with the onset of the experimental absorption spectra [8][9][10]. The SOC matrix element between the 1 CT and 3 LE states were computed at the ground state equilibrium geometry using the perturbative approach developed by Wang and Ziegler [11] as implemented within ADF [12][13][14]. Using the M062X functional and a TZP basis set for all atoms. Scalar relativistic effects were included using the ZORA [15,16] approximation. The SOC matrix element between the two charge transfer states was zero at the Franck-Condon point. The Hyperfine interaction matrix element is known to be small and throughout this work a constant value of 0.2 cm −1 , in constant with typical experimentally recorded values is used [17].

State Symmetry Energy (eV)
Character A" 3.53 HOMO→LUMO where Q i is the nuclear degree of freedom and E rel is the excited state energy, relative to the 3 LE state. This means that E rel 3 LE will be 0, while for E rel 3 CT and E rel 1 CT , the energy represents the gap between this state and the 3 LE state.   Figure 3 in the main text. All energies are in cm −1 unless stated otherwise. In green are highlighted the parameters which are adjusted during each simulation

S2 Quantum Dynamics Simulations
The quantum dynamics addressing the intersystem crossing rate were performed using the Heidelberg Multi Configuration Time Dependent Hartree (MCTDH) package [18,19]. In this approach the wavefunction ansatz is written as a linear combination of Hartree products: where Q 1 ,...,Q f are the nuclear coordinates, A j 1 ...j f (t) are the time-dependent expansion coefficients and ϕ (k) j k are the time dependent basis functions for each k (degree of freedom), known as single particle functions (SPFs). The SPFs used in MCTDH have two advantages: (1) fewer are required as they are variationally determined (2) the functions can be multi-dimensional particles containing more than one degree of freedom thus reducing the effective number of degrees of freedom. The wavepacket simulation describes the evolution of a certain, well defined, initial state. However for a system at finite temperature, important in the context of the simulations addressing the rISC rate, there obviously exists a mixture of different thermally excited states.
To address this we perform simulations within a density operator formalism of MCTDH [20]. Here the single particle functions are replaced with single-particle density operators. Here we adopt a closed quantum system, this is to say that no dissipative operators are included and only the core Hamiltonian described above is used. In this representation the Liouville-von Neumann equation for the system is expressed: For these simulations, the advantage of MCTDH comes into its own. Although the model Hamiltonian used herein is relatively small, for the density operator simulations the dimensionality of the system formally doubles [21] significantly increasing the numerical treatment of the simulations. The full details of both sets of MCTDH simulations are given in Table S4 and ensure convergence for the population kinetics for the entire simulations Table S4: Computational details for the MCTDH simulations within both the wavefunction and density operator formalisms. N i is the number of primitive harmonic oscillator discrete variable representation (DVR) basis functions used to describe each mode. n i are the number of single-particle functions used to describe the wavepacket on each state.
Modes N i n S 1 ,n T 1 ,n The dynamics were performed using the model Hamiltonian described above. The computational details for the quantum dynamics simulations are shown in Table S1 and ensured convergence of the population kinetics shown in the main text.  Figure S2: (a) Population kinetics of the 3 CT state during the rISC dynamics after initial population of the 3 LE state. These simulations are performed within the density operator formalisms with a temperature of 300 K and correspond to the same simulations as shown in Figure 3b in the main text. (b) Zoom in to the first picosecond of (a).