Thermally activated delayed fluorescence processes for Cu(I) complexes in solid-state: a computational study using quantitative prediction

Abstract

The photophysical properties of four representative Cu(I) complex crystals have been investigated using the combination of an optimally tuned one- and two-dimensional range-separated hybrid functional with the polarizable continuum model, and the thermal vibration correlation function (TVCF) approach. The calculated excited singlet–triplet energy gap, radiative rates and lifetimes match the experimentally available data perfectly. At 300 K, the reverse intersystem crossing (RISC) proceeds at a rate of k^dir._RISC ≈ 10^6–8 s⁻¹, which is 4–5 orders of magnitude larger than the mean phosphorescence rate, k_P ≈ 10^2–3 s⁻¹. At the same time, the ISC rate k^dir._ISC ≈ 10⁹ s⁻¹ is again 2 orders of magnitude larger than the fluorescence rate k_F ≈ 10⁷ s⁻¹. In the case of k^dir._RISC ≫ k_F and k^dir._RISC ≫ k_P, thermally activated delayed fluorescence should occur. Vibronic spin–orbit coupling can remarkably enhance the ISC rates by the vital “promoting” modes, which can provide crucial pathways to decay. This can be helpful for designing novel excellent TADF Cu(I) complex materials.

---

1. Introduction

Recently, the thermally activated delayed fluorescence (TADF) materials displayed by organo-transition metal complexes and organic molecules have attracted great attention because of the remarkable variability in their emission properties for organic light-emitting diodes (OLED).^1–4 Under electrical excitation, the singlet–triplet counter pairs of electrons and holes are equally weighted to recombine and yield excitons, resulting in 25% singlet excitons and 75% triplet excitons in the electroluminescence device according to spin statistics.^5–7 However, the energies of all triplet excitons (75%) are dissipated as heat in the conventional fluorescent materials, which leads to a theoretical upper limit of 25% for the internal quantum efficiency (Fig. 1).^6,7 To obtain high-efficiency OLED materials, recent studies have found that TADF emitters can rely on efficient thermal upconversion from the triplet state T₁ into an emissive singlet state S₁ through reverse intersystem crossing (RISC), which in principle, causes the efficiency of exciton utilization to reach 100%.^8–10

Fig. 1 Diagram of the delayed fluorescence process based on spin-statistics in TADF molecules.

It is well known that efficient TADF must satisfy the vital condition of a small energy gap, ΔE(S₁–T₁), between the T₁ and S₁ states involved in the RISC process. One effective strategy to decrease the ΔE(S₁–T₁) involves using covalently linked electron donor and acceptor units and consequently, such molecules have become the focus of molecular designs adopted for TADF.¹¹ Using this strategy, one obtains the T₁ and S₁ states with strong charge transfer (CT) character from the highest occupied molecular orbital (HOMO) to the lowest unoccupied molecular orbital (LUMO) transitions. These states are characterized by a small spatial overlap between the HOMO and LUMO, which decreases the electron exchange energy and the energy gap, ΔE(S₁–T₁).^11,12 However, in this simple model, one negative effect is inevitably generated: the small spatial 〈HOMO|LUMO〉 overlap for S₁ with CT causes the slow radiative decay rate of the S₁ → S₀ transition and low luminescence efficiency. Therefore, it is an essential challenge to obtain a reasonable balance between a small ΔE(S₁–T₁) and the oscillator strength of the S₁ → S₀ transition in the TADF molecular design. An accurate evaluation of ΔE(S₁–T₁) with a reasonable oscillator strength is necessary to guide the TADF molecular design and expand the range of TADF materials.

Moreover, it is well established that the optical properties of TADF materials can have a considerable impact as a function of their environment. In other words, TADF materials are usually mixed into appropriate host matrixes to remit the concentration quenching and exciton annihilation processes, which restrain the efficiency enhancement of OLEDs. The study reported by Yersin et al.¹³ has shown that the same TADF complex can have a much larger photoluminescence quantum yield in the solid state than that in solution. Therefore, it is also a challenge to be able to accurately describe the solid phase excited-state properties of TADF emitters. However, the theoretical studies of TADF materials, taking into account a solid-state environment, have great limitations so far.

From a theoretical standpoint, a reliable and computationally efficient method for the prediction of the excited state properties for the solid state would thus be highly beneficial since it would allow one to possibly help in determining the origin of the variation in experimental results. The time-dependent density functional theory (TD-DFT) is a useful and reliable tool to compute the excited states of relatively larger systems and is considerably more computationally efficient.¹⁴ Unfortunately, the conventional (semi)local exchange–correlation (XC) functionals may fail completely in predicting the electronic structure in donor–acceptor CT systems.¹⁵ In addition, for molecular crystals, the surrounding environment is substantially different from that of single molecules, ascribed to polarization effects, which is essentially a phenomenon related to nonlocal correlation. These systematic errors are mainly attributed to the inappropriate XC introduction and the potential and density can be incorrect at asymptotically large distances.^16–18 Recently, the range-separated exchange (RS) density functional comprising a suitable fixed amount of exact-exchange (eX) has overcome the incorrect asymptotic behavior in the long-range limit, and provides an improved description of the excited-state properties.^16–21 Abramson et al.²² reported the dielectric constant (ε) in a “screened” RS functional by replacing the 1/R asymptotic behavior with the more general asymptotic 1/(ε × R), which is required when the calculations are performed on the periodic crystals. Very recently, a more computationally efficient method was proposed by Sun et al.²³ for the quantitative characterization of the excited-state properties in molecular crystals. This calculation method is to combine the polarizable continuum model (PCM) and optimally tuned RS functionals, whose advantage comes from the optimal tuning of the range-separation parameter ω for a long-range corrected (LC) functional in the PCM environment. For the simulation of a solid-state environment, the molecules in the crystal phase were “dissolved” in a solvent with the same kind of molecules as the crystal.

Potential candidates for Cu(I) complexes in the fields of OLEDs have attracted a lot of attention from the academic and commercial communities. A much deeper investigation of the photophysical and chemical properties of the Cu(I) compounds will lead to the development of new materials and material design strategies. In this study, we chose four representative TADF emitters that are Cu(I) complexes, namely, Cu(pop)(NN) (pop = bis(2-(diphenylphosphanyl)-phenyl)ether and NN = bis(pyrazol-1-yl)borohydrate(pz₂BH₂), tetrakis(pyrazol-1-yl) borate (pz₄B), bis(pyrazol-1-yl)-diphenyl-borate (pz₂Bph₂) and dppb = 1,2-bis(diphenylphosphino)-benzene; Fig. 2), as our research objects because their extensive experimental data are available.¹³ Herein, the k_r radiative decay rates are determined via the Einstein relationship and the k_ISC decay rates are quantitatively calculated using the thermal vibration correlation function (TVCF) rate theory in combination with the PCM-tuned LC-BLYP method in the solid-state.^24–29 In addition, vibronic spin–orbit coupling (SOC) has been taken into account from the promoted vibration modes.

Fig. 2 Chemical structure of mononuclear Cu(I) complexes.

2. Computation details

2.1 Optimization of geometries

The optimal tuning method was performed according to the RS functionals. In the expression of RS hybrid functionals, the two-electron repulsion operator 1/R₁₂ was divided into a short-range domain and a long-range domain by means of the Ewald-style partition based on the error function (erf) as follows:³⁰where R₁₂ is the distance between the electrons at the coordinate vectors R₁ and R₂; the range-separation parameter ω (units of bohr⁻¹) determines the ratio between the ranges depending on the value of R₁₂. The first term on the right-hand side in eqn (1) is a short-range interaction described by DFT exchange potential, and the second term is the long-range interaction described by the HF exchange integral. For an asymptotically correct LC functional, α + β = 1, which is enforced for the optimal tuning in this study. The optimal tuning procedure is based on the exact Kohn–Sham (KS) theory, in which the negative HOMO energy for an N-electron system equals the vertical ionization potential;²³ this expression is as follows:

In addition, there have been recent studies^31,32 showing that some amount of short-range HF exchange (i.e. setting α to a nonzero value) can lead to improved electronic properties and excitation energies. We used two different parameterizations: a long-range corrected LC-BLYP functional without any short-range exchange (i.e., α = 0.0, β = 1.0) as well as an LC-BLYP functional containing 20% exchange over the entire range (i.e., α = 0.2, β = 0.8) in conjunction with tuning ω via the non-empirical procedure,³³ which are denoted as .

The optimization of the parameter ω was performed by tuning the optDFTw procedure,³⁴ as described in Fig. 3, and all the ground geometries were optimized using the functional and 6-31+G(d) basis set as shown in Fig. 5. The geometries of the excited states were optimized using the corresponding time-dependent functionals. As for the simulation of the environmental solid polarization effects, the default PCM model using the integral equation formalism variant was employed, where the molecule in the crystal phase was “dissolved” in the solvent with the same kind of crystal molecules. These were performed by adding the “SCRF (PCM, solvent = generic, read)” keyword in the Gaussian 09 package.³⁵ The settings for the PCM parameters have been explicitly described in ref. 23b.

Fig. 3 Functions defined in eqn (2) used for determining of the optimally-tuned range-separation parameter ω (bohr⁻¹) at the level in the solid state. (a) Cu(pop)(pz₂BH₂); (b) Cu(pop)(pz₄B); (c) Cu(pop)(pz₂Bph₂); (d) Cu(dppb)(pz₂Bph₂).

Fig. 4 Diagram of the percentage of exact-exchange (%HF) versus the interelectronic distance (R₁₂) for the functional in the four Cu(I) complexes.

Fig. 5 Calculated mean absolute deviations (MAD) of the vital bond lengths and angles at the equilibrium geometries of the ground state S₀ for the Cu(I) complex molecules in the solid state compared to the experimental values. The MAD values were calculated with respect to the corresponding experimental values, e.g. , and .

2.2 Calculations of the electronic-structure and energy gap ΔE(S₁–T₁)

High level of accuracy of calculations is crucial for the quantitative theoretical prediction of TADF molecules. As described in previous studies, the vertical excitation energies of the S₀ → S₁, E_V(S₁) and S₀ → T₁, E_V(T₁) and the energy gap ΔE(S₁–T₁) = E_V(S₁) − E_V(T₁) were calculated by applying the PCM-tuned and PCM-tuned functionals, combined with the 6-31+G(d) basis set within the Tamm–Dancoff approximation (TDA).³⁶ For comparison, two exchange–correlation functionals, namely, B3LYP (α = 0.2) and M062X (α = 0.54), and two RS hybrid functions, CAM-B3LYP (ω = 0.33, α = 0.19, α + β = 0.65) and ωB97XD (ω = 0.20, α = 0.22, α + β = 1.0), were used.^23b,25 To give a quantitative characterization of the excited states, we also calculated the distance between the centroid of the hole and electron, ΔR_H–L, hole–electron overlap integral S_H–L, and CT excitation energy E_CT corresponding to the S₀ → S₁ transition. These calculations were conducted at the theoretical level of PCM(ε)-tuned using the Multiwfn program.³⁷

2.3 Calculations of fluorescence and phosphorescence rates

The phosphorescence transitions (T₁ → S₀) are strictly forbidden in the regime of non-relativistic treatment. If the SOC is considered, the pure singlet or triplet becomes a mixed state. In this case, the forbidden emission can borrow dipole activity from spin-allowed transitions (S₀ ↔ S_n and T_1,ζ ↔ T_m,ζ; see eqn (3)) through the perturbation SOC interactions, resulting in the non-zero intensity of transitions.³⁸ These calculations were performed using the parallel version of the combined PCM-tuned functional and restricted open-shell configuration interaction with the single excitation (DFT/ROCIS) method by the ORCA 4.0 package.³⁹ DFT/ROCIS calculations with the parameters⁴⁰ c₁ = 0.21, c₂ = 0.49, and c₃ = 0.29 were performed using the converged restricted Kohn–Sham functions at the PCM-tuned level, together with the auxiliary basis set def2/J in order to accelerate the calculations in the framework of the RI approximation. The electric dipole transition, , between T₁ and S₀ becomes allowed and can be expressed as follows:^41,42where the operators μ_α and represent the electric dipole and SOC Hamiltonian, respectively. At ambient temperature, the emission rates can be calculated by eqn (4):^41,42where indicates a vertical emission energy.

2.4 Calculations of ISC rates

ISC rates of the T₁ ↔ S₁ mainly depend on the energy gap (ΔE(S₁–T₁)) as well as on the strength of the SOC between T₁ and S₁ states. However, in TADF Cu(I) complexes, their SOC interaction is strictly forbidden due to the same CT electronic structures. Therefore, the vibronic SOC contribution to the ISC and RISC rates should be considered,⁴³ and these calculations were conducted based on the TVCF theory using the MOMAP suite of programs.^24–29 The potential energy surface of S₁ was considered by using Q_S1 = SQ_T1 + D, relative to that of T₁. The vector D is the displacement between S₁ and T₁ geometries and S is the Duschinsky rotation matrix representing the mixing of normal modes in the S₁ and T₁ states; Q represents the nuclear normal mode coordinates. In the Franck–Condon approximation, the ISC rate constant at temperature T is then given by eqn (5).^26,44here, is the canonical partition function; v and u are the vibrational quantum numbers of T₁ and S₁ states. Further details of the formula can be found in ref. 26.

3. Results and discussion

3.1 Parameter ω optimization and the geometric structures

The optimally tuned derived ω values are given in Fig. 3 for the four Cu(I) molecular systems in the solid states; i.e., using the PCM model with the respective screening effects of the dielectric constant of the molecular crystal during the tuning process. Under the environment of the larger dielectric constants, the optimal ω values were greatly reduced to the range of 0.0491–0.0559 bohr⁻¹ for the solid phase system, which is compared with the default ω = 0.47 bohr⁻¹ for LC-BLYP. As discussed in detail elsewhere, since ω reflects the global delocalization degree, there is an inverse relationship between the tuned ω value and the spatial extension of the delocalization. The smaller ω values demonstrate that the electron density in the simulated crystal environment is more delocalized than that for the isolated molecule in the gas phase.

In addition, from the perspective of the equation of RS functionals as shown in eqn (1), a smaller ω value corresponds to a larger interelectronic distance, R₁₂, where the description of exchange switches from the short-range DFT-type to the long-range exact-exchange (HF-type) or, in other words, an effective electron delocalization length; the previously optimal ω values were found to decrease with increasing system size and conjugation length. However, in some instances, this correlation was not monotonous and varied strongly for systems with different electronic structures. Here, the Cu(I) complexes have different sizes, being smallest for Cu(pop)(pz₂BH₂) (ω = 0.0559 bohr⁻¹) and largest for Cu(pop)(pz₂Bph₂) (ω = 0.0495 bohr⁻¹) and Cu(dppb)(pz₂Bph₂) (ω = 0.0491 bohr⁻¹), which has a strong dependence of ω on the size of ligands around the central Cu(I) atom.

To illustrate the above behavior, the dependence of the percentage (%HF) of eX in LC-BLYP and versus the R₁₂ for four Cu(I) complexes is displayed in Fig. 4. We take the tuned ω value of the functional for the Cu(dppb)(pz₂Bph₂) system as an example. At R₁₂ = 2.5 a.u. for the solid state, the HF exchange of the optimally tuned functional is 14.6% (ω = 0.0491 bohr⁻¹), whereas the default LC-BLYP (ω = 0.47 bohr⁻¹) gives more than 89.9% HF exchange. This indicates that the precise description of the solid phase requires the functionals to include less “localized” HF exchange and more “delocalized” DFT-BLYP exchange for the Cu(I) complexes. The eX percentage functional of the other three complexes is similar to that of Cu(dppb)(pz₂Bph₂).

The ground state (S₀) geometries of the solid state molecules were optimized using the optimally tuned functional and conventional functional B3LYP. These structural parameters are illustrated in Fig. 5 and Table S1 in the ESI.† It is easily noted that the optimal tuning of functionals is very important for the geometric optimization of S₀. For the B3LYP, deviations in the optimized results are much larger. For example, the mean absolute deviations (MADs) in variation between the lower and the upper deviations are 0.0563 Å and 0.0851 Å, respectively, for the bond lengths in the four Cu(I) complexes and between 1.41° and 2.53° for the bond angles, as compared with the crystal data. For the PCM-tuned (ω = 0.0559 bohr⁻¹) functional, the MADs were merely 0.0224 Å for the bond lengths and 0.51° for the bond angles in the Cu(pop)(pz₂BH₂) complexes; these estimated geometric parameters are in good agreement with experimental crystal data. In this study, we ultimately chose the PCM-tuned (TD-) level, unless otherwise stated, to optimize the geometry structures (including the S₀, S₁ and T₁ states) of the Cu(I) molecular crystals.

3.2 Energy gap ΔE(S₁–T₁) and transition properties

A concise summary of all four Cu(I) molecule excited energies obtained by TD-DFT (B3LYP, M062X, ωB97XD, CAM-B3LYP, ) is given in Table 1. Overall, the two RS functionals provide significantly better predictions for the excited state energies and properties than the conventional hybrid functionals. Compared to the experimental ΔE(S₁–T₁) values, the energy gaps of Cu(I) molecules are most accurately calculated by , including some short-range exchange, and least accurately predicted by B3LYP (this can be attributed to the large electron self-interaction error of the B3LYP functional, which spuriously favors the CT character of the excitations). We also found that significantly improves the accuracy of these excitations in comparison to the results. The calculated values of ΔE(S₁–T₁), 1129, 1048, 887, and 565 cm⁻¹ for the four Cu(I) molecules, namely, Cu(pop)(pz₂BH₂), Cu(pop)(pz₄B), Cu(pop)(pz₂Bph₂), and Cu(dppb)(pz₂Bph₂), are in very good agreement with the measured values of 1300, 1000, 800, and 370 cm⁻¹, respectively in the powder.¹⁴ In the following investigations, was selected to predict the excited state properties, corresponding to optimal ω values listed in Table 2.

Table 1 Calculated vertical excited energies (eV), E_v(S₁) and E_v(T₁), and the corresponding energy difference ΔE(S₁–T₁) at the S₀ minimum using the different functional methods, and compared to the experimental values^a

Methods	Cu(pop)(pz2BH2)			Cu(pop)(pz4B)			Cu(pop)(pz2Bph2)			Cu(dppb)(pz2Bph2)
	Ev(S1)	Ev(T1)	ΔE(S1–T1)	Ev(S1)	Ev(T1)	ΔE(S1–T1)	Ev(S1)	Ev(T1)	ΔE(S1–T1)	Ev(S1)	Ev(T1)	ΔE(S1–T1)
a The experimental values of ΔE(S1–T1): 1300, 1000, 800, and 370 cm^−1 for Cu(pop)(pz2BH2), Cu(pop)(pz4B), Cu(pop)(pz2Bph2), and Cu(dppb)(pz2Bph2), respectively; vertical excited energies between about 3.35 and 3.99 eV.
B3LYP	3.34	3.20	0.14(1129)	3.34	3.22	0.12(968)	3.27	3.16	0.11(887)	3.03	2.96	0.07(565)
M062X	4.42	4.10	0.32(2581)	4.43	4.03	0.40(3226)	4.39	4.09	0.30(2420)	3.83	3.54	0.29(2339)
ωB97XD	4.29	3.89	0.40(3226)	4.30	3.89	0.41(3307)	4.26	3.88	0.38(3065)	3.83	3.51	0.32(2580)
CAM-B3LYP	4.28	3.81	0.47(3791)	4.29	3.80	0.49(3952)	4.25	3.80	0.45(3629)	3.81	3.50	0.31(2500)
	2.51	2.45	0.06(484)	2.52	2.47	0.05(403)	2.42	2.35	0.07(565)	2.28	2.23	0.05(403)
	3.34	3.20	0.14(1129)	3.35	3.22	0.13(1048)	3.27	3.16	0.11(887)	3.03	2.96	0.07(565)

---

Table 2 Optimal range-separated parameters ω (bohr⁻¹), dielectric constant ε, calculated ionization potential (IP), electron affinity (EA), energetic gap (E_g), CT excitation energy (E_CT) of the molecules in the solid phase; energy is in eV

Species	PCM(ε)-tuned
	ε	ω	IP(-εH)	EA(-εL)	Eg	ECT^a
a The CT excitation energy: ECT = IP − EA − 1/R; εH = HOMO energy, εL = LUMO energy.
Cu(pop)(pz2BH2)	3.19	0.0169	5.01	0.59	4.42	0.38
Cu(pop)(pz4B)	3.23	0.0172	5.10	0.68	4.42	0.30
Cu(pop)(pz2Bph2)	3.49	0.0144	4.93	0.65	4.28	0.14
Cu(dppb)(pz2Bph2)	3.66	0.0131	4.84	0.78	4.06	−0.63

---

For the studied complexes, their first S₁ and T₁ excited states mainly come from the electronic transitions from HOMO to LUMO, and they possess obvious CT character and have the same configurations, ¹CT and ³CT (see Fig. 6). Based on the quantum theory, the S₁ and T₁ states are separated by twice the electron exchange energy J, as illustrated in eqn (6) and (7) with ϕ_L and ϕ_H corresponding to the electron (LUMO) and hole (HOMO) orbitals involved in the transitions, respectively. The J value is determined by the spatial separation and overlap integral of ϕ_L and ϕ_H, the higher overlap of HOMO and LUMO and smaller spatial separation lead to higher J and ΔE(S₁–T₁).¹²

ΔE(S₁–T₁) = E(S₁) − E(T₁) = 2J (6)

Fig. 6 Computed hole (HOMO)-electron (LUMO) indexes of the mononuclear Cu(I) complexes for selected S₁ states at the PCM(ε)-tuned theoretical level. The distance between the centroid of the hole and electron, ΔR_H–L is in Å, the hole–electron overlap integral S_H–L is in a.u.; ε_H, HOMO energy; ε_L, LUMO energy; E_g, energetic gap; E_CT, CT excitation energy; all energy units are in eV.

Thus, to gain a quantitative understanding, parameters such as the extent of hole–electron overlap (S_H/L) and mean separation distance (ΔR_H–L) of HOMO and LUMO associated with the S₁ transitions were calculated using the Multiwfn program,³⁷ and the calculated results are given in Fig. 6. As expected, a correlation between the ΔE(S₁–T₁) and S_H/L was observed for the solid state Cu(I) complexes, except for Cu(dppb)(pz₂Bph₂) and those with a smaller overlap S_H/L and therefore a smaller ΔE(S₁–T₁). The S_H/L overlap gradually decreases, 0.1241 > 0.1189 > 0.1119 for Cu(pop)(pz₂BH₂), Cu(pop)(pz₄B), and Cu(pop)(pz₂Bph₂), respectively, and the corresponding ΔE(S₁–T₁) values also have a decreasing trend with 1129 > 1048 > 887 cm⁻¹. For the S₁ state of Cu(dppb)(pz₂Bph₂) in the solid state (ε = 3.66), the calculated S_H/L and ΔR_H–L values are 0.1557 and 3.06 Å, respectively, which indicate that it should have a larger ΔE(S₁–T₁); however, contrary to the theoretical analysis, it had a smaller ΔE(S₁–T₁) of 565 cm⁻¹. One key factor is that this does not consider the localisation and spatial confinement of the HOMO–LUMO. This degree of spatial confinement is very crucial, which will alter the coulombic and exchange energies. Indeed, since the S₁ state of Cu(dppb)(pz₂Bph₂) is more delocalised than the three other S₁ excited states, it will exhibit the smallest ΔE(S₁–T₁). This is not addressed by the relation ΔE(S₁–T₁) ∝ S_H/L.

An understanding of the degree of spatial confinement can be gained from the energetic gap (E_g) between the HOMO and the LUMO and the CT excitation energy (E_CT), which becomes larger with increasing localisation. According to the description, we would consequently expect that the larger E_g and E_CT would lead to the larger ΔE(S₁–T₁) as confinement increases for the Cu(I) complexes and this trend is reflected in Table 2. Here, a clear trend is found and those with a larger E_g and a larger ΔE(S₁–T₁) are observed.

3.3 Phosphorescence and fluorescence

For excited T₁ states, if their phosphorescence and non-radiative decay are slow and the energy gap ΔE(S₁–T₁) is small enough, then after vibration thermalisation, a reverse of ISC back to the S₁ can occur followed by delayed fluorescence; therefore, a detailed study of the T₁ states is very important. We know that the degeneracy of the three spin sublevels (T_1,I, T_1,II, and T_1,III) can be lifted due to the presence of zero-field splitting (ZFS) due to the SOC interaction, which can lead to quite different radiative and non-radiative properties. Yersin et al.^13,45 have reported the experimental ZFS values in the range of 1 to 10 cm⁻¹ for Cu(I) complex crystals. The results of the calculated ZFS parameters show that the D values are in the range of 0.1 to 4.7 cm⁻¹ at the PCM-tuned levels, consistent with the experimental values (see Table 3).

Table 3 Calculated spin–orbit coupling values (cm⁻¹) between T₁ and S₁ for the four Cu(I) complexes and the ZFS tensors D and E with units in cm⁻¹

Species	SOCx^a	SOCy	SOCz	SOCC^b	D	E
a b
Cu(pop)(pz2BH2)	−72.74	−58.56	11.07	54.29	4.7	0.1
Cu(pop)(pz4B)	11.52	10.45	13.23	11.78	1.4	0.2
Cu(pop)(pz2Bph2)	−1.18	−1.16	−7.11	4.20	2.6	0.4
Cu(dppb)(pz2Bph2)	−0.50	−2.19	0.12	1.29	0.1	0.0

---

Taking into account ZFS, k_P (phosphorescence) and k_F (fluorescence) rates or (τ) lifetimes, have been obtained by the PCM-tuned calculation at the T₁ geometry, including SOC interactions, which are listed in Table 4. At room temperature, due to the usually small ZFS for the four Cu(I) complexes, an average emission decay time τ_av of the three sub-states can be calculated by the three individual decay times according to τ_av = 3(τ_I⁻¹ + τ_II⁻¹+ τ_III⁻¹)⁻¹ (where τ_I, τ_II, and τ_III represent the emission decay times of the T_I, T_II, and T_III substates, respectively).⁴³ The phosphorescence of the four Cu(I) complexes started from T₁ and all the average radiative decay rate constants k_P are relatively slow (less than 10³ s⁻¹), particularly for Cu(dppb)(pz₂Bph₂), for which the k_P was calculated to be 9.01 × 10² s⁻¹ because the decreased rigidity of the dppb ligand leads to a smaller k_P in comparison with the three Cu(I) molecules, including the pop ligand. The calculated means for the three phosphorescence rates are k_P,av = 7.69 × 10³ (τ_av = 130 μs), 1.53 × 10³ (τ_av = 654 μs), 2.02 × 10² (τ_av = 495 μs), and 9.01 × 10² s⁻¹ (τ_av = 1109 μs) for Cu(pop)(pz₂BH₂), Cu(pop)(pz₄B), Cu(pop)(pz₂Bph₂) and Cu(dppb)(pz₂Bph₂), respectively, and the corresponding fluorescence rates k_F (see Table 4) are in good agreement with the experimentally observed values.

Table 4 Computed vertical transition energies ΔE, oscillator strength f, radiative rates k_r, and lifetimes τ of the spin sublevels at the T₁ minimum and the corresponding experimental values in parentheses

State	ΔE (cm^−1)	f	kr (s^−1)	τ (μs)
a ref. 13a.b ref. 13b.
Cu(pop)(pz2BH2)
T1,I	22343.9	1.68 × 10^−6	5.59 × 10^2	1789
T1,II	22344.0	9.80 × 10^−8	3.26 × 10^1	30670
T1,III	22348.7	6.76 × 10^−5	2.25 × 10^4	44
Average			kP = 7.69 × 10^3	τav (T1) = 130 (610)^a
S1	22993.3	5.47 × 10^−2	kF = 2.89 × 10^7	τ (S1) = 0.035 (0.03)^a
Cu(pop)(pz4B)
T1,I	22728.3	1.45 × 10^−6	5.28 × 10^2	1894
T1,II	22728.5	8.28 × 10^−7	2.85 × 10^2	3509
T1,III	22729.9	1.09 × 10^−5	3.75 × 10^3	267
Average			kP = 1.53 × 10^3	τav (T1) = 654 (450)^a
S1	23200.0	1.06 × 10^−2	kF = 3.80 × 10^6	τ (S1) = 0.263 (0.2)^a
Cu(pop)(pz2Bph2)
T1,I	20643.3	4.07 × 10^−6	1.16 × 10^3	862
T1,II	20643.7	1.58 × 10^−5	4.49 × 10^3	223
T1,III	20646.3	1.40 × 10^−6	3.98 × 10^2	2513
Average			kP = 2.02 × 10^3	τav (T1) = 495 (480)^a
S1	21380.5	3.44 × 10^−2	kF = 1.05 × 10^7	τ (S1) = 0.095 (0.12)^a
Cu(dppb)(pz2Bph2)
T1,I	22576.1	2.43 × 10^−7	8.26 × 10^1	12106
T1,II	22576.1	1.73 × 10^−6	5.88 × 10^2	1701
T1,III	22576.2	5.96 × 10^−6	2.03 × 10^3	493
Average			kP = 9.01 × 10^2	τav (T1) = 1109(1200)^b
S1	22948.66	2.28 × 10^−2	kF = 8.06 × 10^6	τ (S1) = 0.124 (0.18)^b

---

3.4 ISC and RISC

As discussed previously, owing to the small energy gap ΔE(S₁–T₁) between the S₁ and T₁, TADF should, in principle, be possible. However, to actually take place effectively, S₁ has to be repopulated, which indicates that the RISC rate should be larger than the rates of radiative and nonradiative decay of the T₁ state to the S₀ state. In fact, the S₁ and T₁ states have the same electronic configurations and almost equivalent weights, and their direct SOC interaction is almost forbidden. In this case, vibronic spin–orbit coupling effects should be considered.

Intersystem crossing rate constants were calculated by considering the direct spin–orbit coupling, (k^dir._ISC/RISC) and the vibronic spin–orbit coupling, (k^vib._ISC/RISC) between the T₁ and S₁ states, as shown in Table 5. As can be seen from Tables 4 and 5, for the four Cu(I) complexes, RISC rates are approximately k^dir._RISC ≈ 10^6–8 s⁻¹ at 300 K, which is 4–5 orders of magnitude larger than the mean phosphorescence rate, k_P ≈ 10^2–3 s⁻¹. This indicates that the S₁ state can be redistributed from the T₁ state by the RISC process. At the same time, the ISC rate k^dir._ISC ≈ 10⁹ s⁻¹ is again 2 orders of magnitude larger than the fluorescence rate k_F ≈ 10⁷ s⁻¹. Based on the kinetic analysis of a three-level system proposed by Kirchhoff et al.,⁴⁶ in the case of k^dir._ISC > k_F + k_IC0 and k^dir._RISC > k_P + k_ISC0, the S₁ and T₁ state populations rapidly equilibrate before decaying radiatively at room temperature, and TADF should be observable. Here, the rate constants of the internal conversion of S₁ to S₀ (k_IC0) and of ISC from T₁ to S₀ were not calculated owing to the larger energy difference between S₁ (or T₁) and S₀, assuming that these processes are much slower than the radiative rates and will be omitted in the following discussion.

Table 5 Calculated intersystem crossing rate constants (s⁻¹) considering the direct spin–orbit coupling, (k^dir._ISC/RISC) and the vibronic spin–orbit coupling, (k^vib._ISC/RISC) between the T₁ and S₁ states, respectively

T = 300 K	k^dir.ISC	k^dir.RISC	k^vib.ISC	k^vib.RISC
Cu(pop)(pz2BH2)	4.72 × 10^9	6.97 × 10^6	1.39 × 10^11	1.97 × 10^8
Cu(pop)(pz4B)	5.68 × 10^9	1.23 × 10^7	1.74 × 10^10	3.79 × 10^7
Cu(pop)(pz2Bph2)	6.14 × 10^9	2.93 × 10^7	7.94 × 10^11	3.78 × 10^9
Cu(dppb)(pz2Bph2)	1.57 × 10^9	1.04 × 10^8	4.85 × 10^11	1.08 × 10^10

---

The first equilibrium condition, k^dir._ISC ≫ k_F is easily satisfied at all temperatures, while the second condition, namely k^dir._RISC ≫ k_P, because of the strong dependence on temperature, is not. A spin-forbidden T₁ ↔ S₁ transition requires the effect of SOC that provides a major mechanism for RISC. SOC induces a different spin-multiplicity mixing that allows the wave function to break spin symmetry, but direct SOC between S₁ and T₁ states tends to zero; thus, a kind of geometrical alteration must exist, which either leads to the increase in the spin–orbit interaction or causes direct T₁–S₁ intersection. For the T₁ state, among many vibration modes, only a few vibration modes play a significant role in altering the geometric structure, inducing stronger SOC interaction. In the present study, these favorite modes are called “promoting” modes, which are determined by Huang-Rhys factors and related reorganization energies of the normal modes, i (see Table S2 in ESI†). Their relationship can be expressed as follows:⁴⁷

where S_i and ω_i represent the Huang-Rhys factor and the vibrational frequency for the ith normal mode, respectively; D_i is the coordinate displacement along the mode i. The extent of geometry relaxation between T₁ and S₁ states can be measured by the value of the Huang-Rhys factor. Here, taking the Cu(dppb)(pz₂Bph₂) complex as an example, the calculated S_i and λ_i versus the normal mode i on the corresponding S₁ and T₁ potential surfaces are plotted in Fig. 7; those of the remaining complexes are given in Fig. S2–S4 in the ESI.† We noted that some low-frequency vibration modes have larger S_i, and one of these vibration frequencies has the largest Huang-Rhys factor and reorganization energy, i.e., it is that vital “promoting” mode, which corresponds to the distortion motion of the phenylene ring of the dppb ligand (see Fig. 7d). This implies that the distortion rotation of the phenylene ring can provide an important channel for conversion between T₁ and S₁.

Fig. 7 Calculated Huang-Rhys factors and reorganization energies of each normal mode on the corresponding S₁ (a and b, red) and T₁ (c and d, blue) potential surfaces for Cu(dppb)(pz₂Bph₂). The biggest displacement vector is depicted.

The values of component-averaged 〈T₁|H_SOC|S₁〉 SOCs calculated along the displacement of the vital “promoting” modes with respect to the coordinates at the equilibrium geometry are illustrated in Fig. 8. For Cu(pop)(pz₂BH₂), Cu(pop)(pz₄B), Cu(pop)(pz₂Bph₂) and Cu(dppb)(pz₂Bph₂), the promoting modes ν₇₈, ν₁₆, ν₂₇, and ν₁ with frequencies of 699.6, 63.7, 153.6 and 28.9 cm⁻¹, corresponding to Huang-Rhys factors S_i = 2.07, 1.07, 4.27, and 2.42 (see Table S2 in ESI†), respectively, lead to the remarkable increase in the 〈T₁|H_SOC|S₁〉 SOCs in the ranges of 82.5–105.2, 48.5–55.3, 16.2–65.3, and 4.0–4.9 cm⁻¹, respectively. Recalling that we calculated the direct SOC to be several orders of magnitude lower, the vibronic SOC may thus be suspected to be the dominant coupling interaction here. The results show that vibronic SOC can remarkably enhance the ISC rates between T₁ and S₁. The largest ISC rates, k^vib._ISC, including vibronic SOC, were actually obtained for the non-radiative transition T₁ ↔ S₁ and were of the order of ≈10^10–11 s⁻¹, which is about 2 orders of magnitude larger than k^dir._ISC. Moreover, the RISC rates from T₁ to S₁ are affected by the same intensity. Therefore, considering the case of the vibronic SOC, TADF should be obviously observable.

Fig. 8 Change in the SOC (〈T₁|H_SOC|S₁〉) upon deformation of the vital promoting vibrational normal modes in the triplet state surfaces for the RISC processes. Derivatives of SOCs with respect to the corresponding (dimensionless) coordinates at the equilibrium geometry, |∂〈T₁|H_SOC|S₁〉|∂q_i, and the vibrational normal frequencies, ν_i, are inserted in the box. (a) Cu(pop)(pz₂BH₂); (b) Cu(pop)(pz₄B); (c) Cu(pop)(pz₂Bph₂); (d) Cu(dppb)(pz₂Bph₂).

4. Conclusions

The accuracy of the calculation level is crucial for the quantitative theoretical prediction of the photophysical properties of TADF molecules. In this study, the geometries and photophysical properties of a series of TADF Cu(I) crystal complexes have been quantitatively investigated using a methodology proposed by Sun et al. This method uses the combination of an optimally tuned one- and two-dimensional range-separated hybrid functional with the polarizable continuum model, in which the solid-state screening effects are achieved via the introduction of the solid-state dielectric constant, ε. When using an optimal tuning ω value, we find that the LC technique gives an outstanding description of the charge-transfer excited state properties. The calculated results have demonstrated that the two-dimensional tuning of the functional provides excellent agreement with experimental data. In contrast, we find that reoptimizing the fraction of the Hartree Fock exchange in conventional hybrid functionals still yields an inconsistent description of excitation energies. For geometries, the one-dimensional tuning of the range-separated hybrid functional performs slightly better than both the B3LYP and functionals.

Photophysical radiation and non-radiation processes by combining an optimal tuning ω value of LC functional methods and TVCF rate theory have been comprehensively investigated for four Cu(I) complexes. At 300 K, the RISC rate was k^dir._RISC ≈ 10^6–8 s⁻¹, which is 4–5 orders of magnitude larger than the mean phosphorescence rate, k_P ≈ 10^2–3 s⁻¹. This indicates that the S₁ state can be redistributed from the T₁ state by the RISC process. At the same time, the ISC rate k^dir._ISC ≈ 10⁹ s⁻¹ was again 2 orders of magnitude larger than the fluorescence rate k_F ≈ 10⁷ s⁻¹. In the case of k^dir._ISC ≫ k_F and k^dir._RISC ≫ k_P, the TADF phenomenon should be observable. Hence, the electronic-vibrational coupling activated by temperature is responsible for the occurrence of delayed fluorescence in Cu(I) complexes. It is worth noting that vibronic spin–orbit coupling can remarkably enhance the ISC rates between T₁ and S₁ by the vital “promoting” modes with the larger Huang-Rhys factor. These “promoting” modes correspond to the distortion motion of the phenylene ring of the pop or dppb ligands, which can provide an important pathway to decay between T₁ and S₁, and would be helpful for designing excellent novel TADF Cu(I) complex materials.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The study was supported by the National Natural Science Foundation of China (Grant No. 21263022, 21663025, 21663024). We gratefully thank Dr Haitao Sun (East China Normal Universty) for providing the help in using the PCM-tuned LC-BLYP* approach.

References

1. H. Xu, R. F. Chen, Q. Sun, W. Lai, Q. Q. Su, W. Huang and X. G. Liu, Chem. Soc. Rev., 2014, 43, 3259–3302.
2. Z. Y. Yang, Z. Mao, Z. L. Xie, Y. Zhang, S. W. Liu, J. Zhao, J. R. Xu, Z. G. Chi and M. P. Aldred, Chem. Soc. Rev., 2017, 46, 915–1016.
3. Y. Tao, K. Yuan, T. Chen, P. Xu, H. H. Li, R. F. Chen, C. Zheng, L. Zhang and W. Huang, Adv. Mater., 2014, 47, 7931–7956.
4. H. Yersin, A. F. Rausch, R. Czerwieniec, T. Hofbeck and T. Fischer, Coord. Chem. Rev., 2011, 225, 2622–2652.
5. T. Ogiwara, Y. Wakikawa and T. Ikoma, J. Phys. Chem. A, 2015, 119, 3415–3418.
6. M. A. Baldo, D. O'brien, Y. You, A. Shoustikov, S. Sibley, M. Thompson and S. Forrest, Nature, 1998, 395, 151–154.
7. F. B. Dias, K. B. Bourdakos, V. Jankus, K. C. Moss, K. T. Kamtekar, V. Bhalla, J. Santos, M. R. Bryce and A. P. Monkman, Adv. Mater., 2013, 25, 3707–3714.
8. C. Adachi, M. A. Baldo, M. E. Thompson and S. R. Forrest, J. Appl. Phys., 2001, 90, 5048–5505.
9. B. Zhao, T. Zhang, B. Chu, W. Li, Z. Su, Y. Luo, R. Li, X. Yan, F. Jin and Y. Gao, Org. Electron., 2015, 17, 15–21.
10. H. Uoyama, K. Goushi, K. Shizu, H. Nomura and C. Adachi, Nature, 2012, 492, 234–238.
11. M. Y. Wong and E. Zysman-Colman, Adv. Mater., 2017, 29, 1605444.
12. T. Chen, L. Zheng, J. Yuan, Z. F. An, R. F. Chen, Y. Tao, H. H. Li, X. J. Xie and W. Huang, Sci. Rep., 2015, 5, 10923.
13. (a) R. Czerwieniec, J. Yu and H. Yersin, Inorg. Chem., 2011, 50, 8293–8301; (b) R. Czerwieniec and H. Yersin, Inorg. Chem., 2015, 54, 4322–4327.
14. J. Föller, M. Kleinschmidt and C. M. Marian, Inorg. Chem., 2016, 55, 7506–7516.
15. H. T. Sun, C. Zhong and J.-L. Brédas, J. Chem. Theory Comput., 2015, 11, 3851–3858.
16. J. Autschbach and M. Srebro, Acc. Chem. Res., 2014, 47, 2592–25602.
17. D. J. Tozer, J. Chem. Phys., 2003, 119, 12697–12699.
18. T. Körzdörfer and J. L. Brédas, Acc. Chem. Res., 2014, 47, 3284–3291.
19. S. P. Huang, Q. S. Zhang, Y. Shiota, T. Nakagawa, K. Kuwabara, K. Yoshizawa and C. Adach, J. Chem.Theory. Comput., 2013, 9, 3872–3877.
20. U. Salzner and A. Aydin, J. Chem. Theory Comput., 2011, 7, 2568–2583.
21. T. Stein, L. Kronik and R. Baer, J. Am. Chem. Soc., 2009, 131, 2818–2820.
22. S. Refaely-Abramson, S. Sharifzadeh, M. Jain, R. Baer, J. B. Neaton and L. Kronik, Phys. Rev. B, 2013, 88, 081204.
23. (a) H. T. Sun, S. Ryno, C. Zhong, M. K. Ravva, Z. R. Sun, T. Körzdörfer and J.-L. Brédas, J. Chem. Theory Comput., 2016, 12, 2906–2916; (b) Z. Hu, B. Zhou, Z. R. Sun and H. T. Sun, J. Comput. Chem., 2017, 38, 569–575.
24. Q. Peng, Y. P. Yi, Z. G. Shuai and J. S. Shao, J. Chem.Phys., 2007, 126, 114302.
25. Y. L. Niu, Q. Peng, C. M. Deng, X. Gao and Z. G. Shuai, J. Phys. Chem. A, 2010, 114, 7817–7831.
26. Q. Peng, Y. L. Niu, Q. Shi, X. Gao and Z. G. Shuai, J. Chem. Theory Comput., 2013, 9, 1132–1143.
27. Y. L. Niu, W. Q. Li, Q. Peng, H. Geng, Y. P. Yi, L. J. Wang, G. J. Nan, D. Wang and Z. G. Shuai, Mol. Phys., 2018, 116, 1078–1090.
28. L. L. Lv, K. Yuan and Y. C. Wang, Org. Electron., 2017, 51, 207–219.
29. L. L. Lv, K. Yuan and Y. C. Wang, Phys. Chem. Chem. Phys., 2018, 20, 6548–6561.
30. T. Yanai, D. P. Tew and N. C. Handy, Chem. Phys. Lett., 2004, 393, 51–57.
31. A. E. Raeber and B. M. Wong, J. Chem. Theory Comput., 2015, 11, 2199–2209.
32. M. A. Rohrdanz, K. M. Martins and J. M. Herber, J. Chem. Phys., 2009, 130, 054112.
33. D. A. Egger, S. Weismann, S. Refaely-Abramson, S. Sharifzadeh, M. Dauth, R. Baer, S. Kummel, J. B. Neaton, E. Zojer and L. Kronik, J. Chem. Theory Comput., 2014, 10, 1934–1952.
34. T. Lu, optDFTw and scanDFTw program v1.0, https://sobereva.com/346, (2017-3-8).
35. M. J. Frisch, G. W. Trucks and H. B. Schlegel, et al., Gaussian 09, Revision-D.01, Gaussian Inc, Wallingford, CT, 2009.
36. S. Hirata and M. Head-Gordon, Chem. Phys. Lett., 1999, 314, 291–299.
37. T. Lu and F. Chen, Multiwfn: a multifunctional wavefunction analyzer, J. Comput. Chem., 2012, 33, 580–592.
38. K. Sugisaki, K. Toyota, K. Sato, D. Shiomi and T. Takui, J. Phys. Chem. A, 2016, 120, 9857–9866.
39. F. Neese, The ORCA program system, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 73–78.
40. M. Roemelt and F. Neese, J. Phys. Chem. A, 2013, 117, 3069–3083.
41. B. Minaev, G. Baryshnikov and H. Agren, Phys. Chem. Chem. Phys., 2016, 16, 1719–1758.
42. G. Baryshnikov, B. Minaev and H. Ågren, Chem. Rev., 2017, 117, 6500–6537.
43. T. J. Penfold, E. Gindensperger, C. Daniel and C. M. Marian, Chem. Rev., 2018 DOI:10.1021/acs.chemrev.7b00617.
44. Q. Peng, D. Fan, R. H. Duan, Y. P. Yi, Y. L. Niu, D. Wong and Z. G. Shuai, J. Phys. Chem. C, 2017, 121, 13448–13456.
45. R. Czerwieniec, M. J. Leitl, H. H. H. Homeier and H. Yersin, Coord. Chem. Rev., 2016, 325, 2–28.
46. J. R. Kirchhoff, R. E. Gamache Jr, M. W. Blaskie, A. A. Del Paggio, R. K. Lengel and D. R. McMillin, Inorg. Chem., 1983, 22, 2380–2384.
47. P. K. Samanta, D. Kim, V. Coropceanu and J.-L. Brédas, J. Am. Chem. Soc., 2017, 139, 4042–4051.

---

Footnote

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

---