As is well-known, density functionals generally have difficulties with simultaneously describing local excitation (LE) and CT states with good accuracy. Since we could only afford to do the geometry optimizations and frequency calculations under the DFT and TDDFT levels, a suitable functional that qualitatively reproduces the SDSPT2 excitation energies has to be chosen by comparing the TDDFT vertical absorption energies of a few common functionals with SDSPT2 data. B3LYP and PBE0 are generally common choices for the excited states of metalloporphyrins, and BP86 is often used to optimize their ground-state structures. Pure functionals usually tend to underestimate excitation energies, but empirically, their description of the Q band (an LE state) is better than many hybrid functionals, as will be confirmed by our calculation results. As CT states are involved in the relaxation process of the excited states of copper porphyrin, range-separated hybrid functionals (which provide good descriptions of CT states in general) may prove to be suitable as well. These considerations gave a list of four representative functionals, BP86, B3LYP, PBE0 and ωB97X, that were subjected to benchmark calculations.

Different functionals display distinct behaviors for the excitation energies of CuP compared to the results obtained from SDSPT2, as shown in Figure 3. The two characteristic absorption bands of the porphyrin molecule correspond to the ²S₁ (Q band) and ²S₂ (B band) states, which are the only bright states of most porphyrin complexes in the visible region. They are also the only excited states for which accurate experimental vertical absorption energies are available: in benzene they have been measured as 2.25 and 3.15 eV, respectively (Eastwood and Gouterman, 1969). Moreover, the absorption energy of the ²T₁ state has been measured by fluorescence excitation spectra experiments, but only for certain substituted porphyrins: for example, the ²T₁ absorption energy of CuEtio (Etio = etioporphyrin I) was measured in n-octane as 1.81 eV, while the emission energy from the same state in the same solvent was 1.79 eV (Eastwood and Gouterman, 1969). Assuming that the Stokes shift of the ²T₁ state is independent of the porphyrin substituents, and combined with the experimental emission energy of the ²T₁ state of CuP in the same solvent (1.88 eV) (Eastwood and Gouterman, 1969), we obtain an estimate of the experimental ²T₁ absorption energy of CuP as 1.90 eV. Gratifyingly, the SDSPT2 excitation energies of all three states agree with the experimental values to within 0.3 eV, which is typical of the accuracy of SDSPT2 (Song et al., 2022) and confirms the suitability of SDSPT2 as a benchmark reference for CuP. The B3LYP functional performs better for the two bright states ²S₁ and ²S₂ than the other functionals (Figure 3), with results closer to the SDSPT2 calculations, suggesting its suitability for localized excitations in the porphyrin system. However, it (as well as the pure functional BP86) performs poorly in describing the dark charge transfer (CT) states, significantly underestimating their energies, as expected. In contrast, the range-separated functional ωB97X shows better agreement with the CT states compared to SDSPT2 results, but its description of the ²S₂ state is rather poor, with energies notably higher than the SDSPT2 results. The PBE0 functional represents a compromise between the two classes of functionals and provides more accurate overall descriptions of the LE and CT states, giving results closer to the SDSPT2 calculations. In particular, PBE0 gives the best predictions of all the d-d excitation energies, although a significant overestimation is still seen for the lowest d-d state, ²dd₁. Considering the overall performance in describing different states, we chose to use PBE0 for the remaining part of the present study.

The ²S₁ and ²S₂ states are almost spin-adapted states with minimal spin contamination, even at the U-TDDFT level (Table 1), since they are dominated by singdoublet excitations. As shown in Figure 3, both X-TDDFT and U-TDDFT provide similar descriptions for these two states; note however that functionals with large amounts of HF exchange generally overestimate the excitation energies of these two states, especially ²S₂. At the TDDFT levels, the CT states are dominated by CO-type excitations (from π to 3 d x 2 − y 2 ), which are also spin-adapted. Both U-TDDFT and X-TDDFT show comparable performance in describing the CT states. However, both methods display large errors compared to SDSPT2 for the CT states. Table 1 presents the excitation energies and the corresponding dominant excited state compositions, computed at the ground state structure of CuP. It can be observed that the CT states are predominantly composed of double excitations, which cannot be described by the present TDDFT calculations. Despite this, functionals with large amounts of HF exchange still perform notably better, as is generally expected for CT states. The ²T₁ and ²T₂ states correspond to tripdoublet excitations (from π to π*), and they suffer from significant spin contamination at the U-TDDFT level, since instead of pure doublets, U-TDDFT can only describe these tripdoublet states as a heavy mixture of doublets and quartets, e.g.: Ψ T 1 2 U - TDDFT ≈ − 1 3 Ψ T 1 2 X - TDDFT + 2 3 Ψ T 1 4 , M S = 1 / 2 X - TDDFT , (5)

As follows from Eqs 2–4. U-TDDFT thus systematically underestimates the excitation energies of the ²T₁ and ²T₂ states, since the energies of quartets are in general lower than the corresponding tripdoublets, as discussed in the Introduction. In Section 3.3 we will also see that part of the underestimation is due to the failure of U-TDDFT to reproduce the energy degeneracy of Ψ ( T 1 4 , M S = 1 / 2 ) and Ψ ( T 1 4 , M S = 3 / 2 ) . On the other hand, X-TDDFT avoids spin contamination through implicitly incorporating extra double excitations necessary for spin-adapting the tripdoublet states (Eq. 3), and therefore performs systematically better than U-TDDFT for all the functionals studied herein. The improvements of the excitation energies (∼0.05 eV) may seem small, but have profound influences on the magnitude and even the sign of the ²T₁–⁴T₁ gap, and therefore on the ratio of fluorescence and phosphorescence emission, as will be detailed in Section 3.3. Finally, X-TDDFT also improves upon the U-TDDFT excitation energies of the d-d excited states, despite that the latter are almost free from spin contamination.Already from the calculated absorption energies, one can draw some conclusions about the photophysical processes of CuP. The vertical absorption energy differences between the ²T₁ and ²S₁ states, as well as between the ²S₂, ²CT₁ and ²CT₂ states, are very small (0.1–0.2 eV). Therefore, once CuP is excited to the bright ²S₁ state by visible light, the molecule is expected to undergo a cascade of ultrafast IC processes, all the way till the doublet state, ²T₁. While the fate of the ²S₂ state is not immediately obvious from the excitation energies, it is well-known that the S₂ states of porphyrin complexes typically undergo ultrafast IC to the S₁ state on the hundred fs to 1 ps timescale (Bräm et al., 2019), suggesting that the ²S₂ state of CuP will also quickly relax to the ²T₁ state via ²S₁, even without the presence of ²CT₁; the presence of ²CT₁ creates a further intermediate energy level between ²S₁ and ²S₂ and is expected to accelerate the IC further. The availability of an ultrafast IC cascade also means the ISC from these high-lying excited states are probably unimportant, especially considering that copper is a relatively light element. These findings are in qualitative agreement with the experimental observation that the ²S₂ states of substituted copper(II) porphyrins relax to the ²T₁ states in gas phase through a two-step process via the intermediacy of a CT state, with time constants 65 fs and 350–2,000 fs, respectively, depending on the substituents (Ha-Thi et al., 2013); the fact that the CT →²T₁ IC is slower than the ²S₂ → CT IC is consistent with the trend of our calculated energy gaps. In solution, the ²S₁ state of Cu(II) protoporphyrin IX dimethyl ester was known to relax to ²T₁ within 8 ps (Kobayashi et al., 2008), and for CuTPP as well as CuOEP the same relaxation was also found to occur within the picosecond timescale (Magde et al., 1974). Recently, the decay rates of the ²S₁ state were measured as 50 fs and 80 fs for CuTPP and CuOEP, respectively, in cyclohexane (Bräm et al., 2019). The ²S₁ state lifetime of CuP itself was also estimated, although indirectly from the natural width of the 0–0 peak of the Q band, as 30 fs (Noort et al., 1976). Quantitative computation of these IC rates is however beyond the scope of the paper, as the narrow energy gaps and possible involvement of conical intersections probably necessitate nonadiabatic molecular dynamics simulations. Nevertheless, a ²S₂ →²CT₁ →²S₁/²T₂ →²T₁ IC pathway can still be tentatively proposed based on the energy ordering alone. Besides, there are also three d-d excited states with energies slightly lower (within 0.34 eV) than the ²T₁ state. There is thus the possibility that the ²T₁ state may further relax to either one of these three d-d states, which we will return to later. Finally, it is worth noting that the use of the accurate SDSPT2 method, as opposed to TDDFT, is crucial for obtaining a reliable estimate of the qualitative trend of the excited state energies. BP86 predicts that the CT states lie below the ²T states, leading to a qualitatively wrong IC pathway; ωB97X, on the other hand, grossly overestimates the energy of ²S₂ and would underestimate its tendency to undergo IC to the CT states (Figure 3). While B3LYP and PBE0 predict reasonable excited state orderings, they still underestimate the energies of the CT states due to the presence of double excitation contributions, which cannot be correctly described under the adiabatic TDDFT framework.

Since all higher lying excited states are predicted to convert to ²T₁ over a short timescale, to study the luminescence of CuP [and probably also other Cu(II) porphyrin complexes bearing alkyl or aryl substituents, given that these substituents do not change excitation energies drastically (Eastwood and Gouterman, 1969)], it should suffice to study the radiative and non-radiative processes starting from the ²T₁ state. Unlike the three ²dd states, whose Franck-Condon emissions are symmetry forbidden, the emission from the ²T₁ state is symmetry allowed, and existing experimental results unequivocally point to a ligand state luminescence, rather than luminescence from a d-d excited state. Therefore, accurately predicting the equilibrium geometry of the ²T₁ state is crucial for subsequent studies.

Some selected bond lengths for the optimized ground state and excited state structures are provided in Table 2. The difference in ground state bond lengths between the UKS and ROKS methods is extremely small ( < 0.0001 Å), as can be seen from their root mean square deviation (RMSD), which can be attributed to the extremely small UKS spin contamination of the ground state of CuP ( ⟨ S 2 ⟩ P B E 0 = 0.7532 ) . The doubly degenerate ²T₁ state, which belongs to the doubly degenerate E _u irreducible representation (irrep) under the D_4h group, undergoes Jahn-Teller distortion to give a D_2h structure, where two of the opposing Cu-N bonds are elongated but the corresponding pyrrole rings remain almost intact, while the other two Cu-N bonds are almost unchanged but the corresponding pyrrole rings exhibit noticeable deformation. The U-TDDFT and X-TDDFT bond lengths of the ²T₁ state show larger deviations than the UKS and ROKS ground state ones, with the largest deviation exceeding 0.001 Å [the C-C (mn) bond], which is also reflected in the RMSD values. However, the structure differences are still small on an absolute scale. This suggests that the coupling of the unpaired Cu(II) d electron and the porphyrin triplet is weak, so that a reasonable tripdoublet state geometry is obtained even if this coupling is described qualitatively incorrectly (as in U-TDDFT). By contrast, our previous benchmark studies on small molecules (where the coupling between unpaired electrons is much larger) revealed that X-TDDFT improves the U-TDDFT bond lengths by 0.01–0.05 Å on average, depending on the functional and the molecule (Li and Liu, 2011; Wang et al., 2020).

As revealed by the above analyses, the relaxation process from high-lying excited states to the ²T₁ state is rapid, and the only energetically accessible relaxation pathways are the radiative (fluorescence) and non-radiative (IC) relaxations from ²T₁ to the ground state ²S₀, the IC from ²T₁ to one of the lower ²dd states, as well as the ISC from ²T₁ to ⁴T₁. The ⁴T₁ state can furthermore convert back to the ²T₁ state through reverse ISC (RISC), or relax to the ground state via radiative (phosphorescence) or non-radiative (ISC) pathways (Figure 4). While the RISC from ⁴T₁ to the ²dd states is also possible, they are expected to be much slower than the spin-allowed ²T₁ →²dd IC processes, and are therefore neglected in the present study.

Before we discuss the quantitative values of transition rates, we first analyze the relevant electronic states from the viewpoint of point group symmetry. The equilibrium structures of the ²T₁ and ⁴T₁ states are both distorted owing to the Jahn-Teller effect, and possess only D_2h symmetry, compared to the D_4h symmetry of the ground state equilibrium structure. The implications are two-fold: the double degeneracy of the ⁿ T₁(n = 2 or 4) state at the D_4h geometry (where they both belong to the E _u irrep) is lifted to give two adiabatic states, hereafter termed the ⁿ T₁(1) and ⁿ T₁(2) states, respectively, where ⁿ T₁(1) is the state with the lower energy; and the potential energy surface of the ⁿ T₁(1) state has two chemically equivalent D_2h minima, ⁿ T₁(1)(X) and ⁿ T₁(1)(Y), where different pairs of Cu-N bonds are lengthened and shortened (see the schematic depictions in Figure 4). Although ⁿ T₁(1)(X) and ⁿ T₁(1)(Y) are on the same adiabatic potential energy surface, their electronic wavefunctions represent different diabatic states, as they belong to the B_3u and B_2u irreps, respectively. The ⁿ T₁(1)(X) structure is diabatically connected to ⁿ T₁(2)(Y) (i.e., the ⁿ T₁(2) state at the equilibrium structure of ⁿ T₁(1)(Y)) via a D_4h conical intersection, while ⁿ T₁(1)(Y) is diabatically connected to ⁿ T₁(2)(X) via the same conical intersection. Thus, the ⁿ T₁(2)(X) and ⁿ T₁(2)(Y) states are expected to undergo ultrafast IC from the D_4h conical intersection, to give the ⁿ T₁(1)(Y) and ⁿ T₁(1)(X) states as the main products, respectively. The direct transition from ⁿ T₁(2) to states other than ⁿ T₁(1) can therefore be neglected. From the irreps of the electronic states, we conclude that certain ISC transitions are forbidden by spatial symmetry. These include the transitions between ²T₁(1)(X) and ⁴T₁(1)(X), between ²T₁(1)(Y) and ⁴T₁(1)(Y), and between any one of the ⁴T₁(1) structures and ²S₀. All IC and radiative transitions, plus the ISC transitions between ²T₁(1)(X) and ⁴T₁(1)(Y) as well as between ²T₁(1)(Y) and ⁴T₁(1)(X), are symmetry allowed. While symmetry forbidden ISC processes can still gain non-zero rates from the Herzberg-Teller effect, we deem that the rates are not large enough to have any noticeable consequences. On one hand, the two symmetry forbidden ISC pathways between the ²T₁(1) and ⁴T₁(1) states are overshadowed by the two symmetry allowed ones, so that the total ISC rate between ²T₁(1) and ⁴T₁(1) is undoubtedly determined by the latter alone. The ISC from ⁴T₁(1) to ²S₀, on the other hand, has to compete with the IC process from ²T₁(1) to ²S₀ in order to affect the quantum yield or the dominant relaxation pathway of the system noticeably, but the latter process is both spin-allowed and spatial symmetry-allowed, while the former is forbidden in both aspects. We therefore neglect all ISC rates whose Franck-Condon contributions are zero by spatial symmetry.

We then calculated the rate constants for all transitions between ²T₁, ⁴T₁ and ²S₀ whose rates are non-negligible, by the TVCF method. The rates (Figure 4) were calculated at 83 K, the temperature used in the quantum yield studies of Eastwood and Gouterman (1969); the latter studies gave a luminescence quantum yield of 0.09, in a solvent mixture of diethyl ether, isopentane, dimethylformamide and ethanol. The accurate treatment of solvation effects is however complicated and beyond the scope of the paper, so that all transition rates were computed in the gas phase. Our calculated k _ISC from ²T₁ to ⁴T₁ (2.36 × 10⁸s⁻¹) is slightly smaller than the total IC rate from ²T₁ to ²S₀ as well as to the three ²dd states (4.27 × 10⁸s⁻¹), suggesting that treating the ²T₁ and ⁴T₁ states as a rapid equilibrium [as in, e.g., Ake and Gouterman (1969) and Bohandy and Kim (1980)] is not justified at least in the gas phase. The IC from ²T₁ to ²dd₃ is two orders of magnitudes slower than other non-radiative relaxation pathways of ²T₁, and is therefore not considered viable in the remaining discussions. At 83 K, the RISC from ⁴T₁ to ²T₁ is 11% of the forward ISC rate. Both rates are in favorable agreement with the experimental values of Cu(II) protoporphyrin IX dimethyl ester in benzene, k _ISC = 1.6 × 10⁹s⁻¹ and k _RISC = 5.6 × 10⁸s⁻¹, at room temperature (Kobayashi et al., 2008). Our computed ISC and RISC rates give a ²T₁-to-⁴T₁ equilibrium concentration ratio of 1:9.0 when all IC processes are neglected, but our kinetic simulation shows that the steady state concentration ratio is 1:50.7 when the latter are considered, further illustrating that treating the ²T₁-⁴T₁ interconversion as a fast equilibrium can lead to noticeable error. Nevertheless, the fluorescence rate of ²T₁ still exceeds the phosphorescence rate of ⁴T₁ by three orders of magnitude, which more than compensates for the low steady state concentration of ²T₁. Similar conclusions could be derived from the rates reported in Ake and Gouterman (1969) (3.6 × 10³s⁻¹ and 8.3 × 10⁻¹s⁻¹, respectively), calculated from semiempirical exchange and SOC integrals and experimental absorption oscillator strengths, which agree surprisingly well with the rates that we obtained here. Kinetic simulation suggests that 99.2% of the total luminescence at this temperature is contributed by fluorescence, and only 0.8% is due to phosphorescence. This can be compared with the experimental finding by Bohandy and Kim (1980) that the phosphorescence of CuP at 86 K is observable as a minor 0–0 peak besides the 0–0 fluorescence peak, with a fluorescence to phosphorescence ratio of about 5:1 to 10:1 [as estimated from Figure 5 of Bohandy and Kim (1980)]; however note that this study was performed in a triphenylene solid matrix.

The total luminescence quantum yield is predicted by our kinetic simulations to be 5.3 × 10⁻⁶, four orders of magnitude smaller than the experimental quantum yield (0.09) in solution. We believe one possible reason is that the ²T₁–⁴T₁ gap of CuP is larger in solution than in the gas phase. This can already be seen from the experimental ²T₁–⁴T₁ 0–0 gaps of CuP in solid matrices with different polarities: the 0–0 gap was measured in polymethylmethacrylate as 500 cm⁻¹ (Smith and Gouterman, 1968), but 310–320 cm⁻¹ in n-octane (Noort et al., 1976) and 267 cm⁻¹ in triphenylene (Bohandy and Kim, 1980). Therefore, the 0–0 gap in the gas phase is probably smaller than 267 cm⁻¹, and indeed, our X-TDA calculations predict an adiabatic ²T₁–⁴T₁ gap of 92 cm⁻¹ in the gas phase. The larger ²T₁–⁴T₁ gap in solution compared to the gas phase is expected to introduce a Boltzmann factor of exp − ( E s o l − E g a s ) / R T to k _RISC, while changing the other rates negligibly. Setting E _sol = 267 cm⁻¹ and E _gas = 92 cm⁻¹, we obtain a solution phase k _RISC of 9.02 × 10⁵s⁻¹, from which kinetic simulations give a fluorescence-phosphorescence ratio of 8.6:1, in quantitative agreement with experiment (Bohandy and Kim, 1980). Setting E _sol = 500 cm⁻¹ [as appropriate for the polar solvent used in Eastwood and Gouterman (1969)] gives k _RISC = 1.60 × 10⁴s⁻¹, and a total luminescence quantum yield of 4.0 × 10⁻⁵, with 13% contribution from fluorescence and 87% from phosphorescence. The remaining discrepancy (∼×2200) of the experimental and calculated quantum yields can be attributed to the restriction of the molecular vibrations of CuP by the low temperature (and thus viscous) solvent, which is expected to suppress the IC process significantly.

Interestingly, U-TDA completely fails to reproduce the qualitative picture of Figure 4 and predicts a ²T₁–⁴T₁ adiabatic gap of the wrong sign (−276 cm⁻¹), violating Hund’s rule. At first sight, this may seem surprising: since the U-TDA “tripdoublet state” is a mixture of the true tripdoublet state and the quartet state, the U-TDA ²T₁ energy should lie in between the energies of the true ²T₁ state and the ⁴T₁ state, which means that the U-TDA ²T₁–⁴T₁ gap should be smaller than the X-TDA gap but still have the correct sign. However, the U-TDA ²T₁ state is contaminated by the M _S = 1/2 component of the ⁴T₁ state (Eq. 5), while a spin flip-up U-TDA calculation of the ⁴T₁ state gives its M _S = 3/2 component. The two spin components obviously have the same energy in the exact non-relativistic theory and in all rigorous spin-adapted methods, but not in U-TDA, even when the ground state is not spin-contaminated (Li and Liu, 2012; Li and Liu, 2016b). This shows that the restoration of the degeneracy of spin multiplets by the random phase approximation (RPA) correction in X-TDA (Li and Liu, 2011) indeed leads to qualitative improvement of the excitation energies, instead of being merely a solution to a conceptual problem. It also shows that estimating the tripdoublet energy by extrapolating from the energies of the spin-contaminated tripdoublet and the quartet by, e.g., the Yamaguchi method (Soda et al., 2000) does not necessarily give a qualitatively correct estimate of the spin-pure tripdoublet energy. The inverted doublet-quartet gap introduces qualitative defects to the computed photophysics of CuP. Already when the doublet-quartet gap is zero, the Boltzmann factor is expected to raise the k _RISC to 9.24 × 10⁷s⁻¹, reducing the ratio of phosphorescence in the total luminescence to 0.11%. Further raising the quartet to reproduce the U-TDA doublet-quartet gap will reduce the k _ISC to 1.42 × 10⁶s⁻¹, which reduces the ratio of phosphorescence to 0.0007%. These values are obviously in much worse agreement with the experiments (Bohandy and Kim, 1980).

Finally, we briefly comment on the luminescence lifetimes. The luminescence of CuP is known to decay non-exponentially (Smith and Gouterman, 1968), so its luminescence lifetime can only be approximately determined. The luminescence lifetime of CuP has been determined as 400 μs (Smith and Gouterman, 1968) at 80 K in polymethylmethacrylate, and a biexponential decay with lifetimes 155 and 750 μs was reported (Eastwood and Gouterman, 1969) at 78 K in methylphthalylethylglycolate. The same references also reported that the luminescence lifetimes of CuOEP and CuTPP are also within the 50–800 μs range. However, in a room temperature toluene solution the luminescence lifetimes of CuOEP and CuTPP were reported to be 115 and 30 ns, respectively (Liu et al., 1995), and a few nanoseconds in the gas phase (Ha-Thi et al., 2013). If we define the luminescence lifetime as the time needed for 1–1/e ≈ 63.2% of the luminescence to be emitted, then kinetic simulations from our X-TDDFT rate constants give a gas-phase luminescence lifetime of 5.7 ns at 83 K, which is much shorter than the low-temperature condensed phase results, but in reasonable agreement with the room-temperature solution phase and especially the gas-phase experimental results. However, using the ISC and RISC rate constants consistent with the U-TDA doublet-quartet gap, one obtains a lifetime of 2.4 ns, differing noticeably from the X-TDDFT result. Thus, our results suggest that X-TDDFT/X-TDA gives non-negligible corrections upon the luminescence lifetime of U-TDDFT/U-TDA, and also confirm that the discrepancy of the experimental and calculated quantum yields is probably due to suppression of the IC of ²T₁ by the low temperature solvent.

As mentioned in the Introduction, the simple orbital energy difference model based on a restricted open-shell determinant (Figure 1) predicts that the excitation energy of the lowest tripdoublet of any doublet molecule is at least the sum of the excitation energies of the first two excited states (as long as the ROKS ground state satisfies the aufbau rule). It therefore comes as a surprise that the lowest tripdoublet state of CuP (²T₁) is only barely higher than the lowest excited state (²dd₁) by 0.34 eV (Table 1), even though the ROKS ground state of CuP is indeed an aufbau state (Figure 5). This suggests a failure of the ROKS orbital energy difference model.

To understand why the ROKS orbital energies fail qualitatively for describing the excited state ordering of CuP, despite that the X-TDDFT method (which uses the ROKS determinant as the reference state) still gives reasonable excitation energies as compared to SDSPT2, we note that the α and β Fock matrices of an ROKS calculation are in general not diagonal under the canonical molecular orbital (CMO) basis. Only the unified coupling operator R, assembled from blocks of the CMO Fock matrices, R = 1 2 F CC α + F CC β F CO β 1 2 F CV α + F CV β F OC β 1 2 F OO α + F OO β F OV α 1 2 F VC α + F VC β F VO α 1 2 F VV α + F VV β , (6)

Is diagonal (Hirao and Nakatsuji, 1973). Note that herein we have used the Guest-Saunders parameterization (Guest and Saunders, 1974) of the diagonal blocks of R, which is the default choice of the BDF program, although our qualitative conclusions are unaffected by choosing other parameterizations. However, the leading term of the X-TDDFT calculation is not simply given by the eigenvalue differences of R, Δ i a σ , j b τ ′ = δ σ τ δ i j δ a b R a a − R i i , (7) but rather from the α and β Fock matrices themselves via Δ i a σ , j b τ = δ σ τ δ i j F a b σ − δ a b F j i σ . (8)

Here i, a represent occupied CMOs, j, b virtual CMOs, and σ, τ spin indices. For the diagonal matrix element of an arbitrary single excitation, Eqs 7, 8 differ by the following term: Δ i a σ , i a σ − Δ i a σ , i a σ ′ = 1 2 F a a σ − F a a σ ′ − F i i σ − F i i σ ′ , (9) where σ′ is the opposite spin of σ. For a general hybrid functional, the Fock matrix element differences in Eq. 9 are given by (where p is an arbitrary CMO, c _x is the proportion of HF exchange, and v ^xc is the XC potential) F p p β − F p p α = c x p t | p t + v p p β x c − v p p α x c , (10)

Assuming, for the sake of simplicity, that there is only one open-shell orbital t in the reference state. Assuming that the XC potential behaves similarly as the exact exchange potential, the difference Eq. 10 is positive, and should usually be the largest when p = t, while being small when p is spatially far from t. The corollary is that the orbital energy difference approximation Eq. 7 should agree well with the X-TDDFT leading term Eq. 8 for CV excitations (where the difference is proportional to the small exchange integral (pt|pt)), but underestimate the excitation energies of CO and OV excitations by a correction proportional to the large (tt|tt) integral.

The underestimation of CO and OV excitation energies by ROKS orbital energy differences opens up the possibility of engineering a system to break the ω _ia = ω _it + ω _ta constraint inherent in the ROKS orbital energy difference model, and make the tripdoublet state the lowest excited state or only slightly higher than the lowest excited state. Possible approaches include:

1. Increase the difference Eq. 9 for the CO and OV states, while keeping it small for the lowest CV state, so that all CO and OV states are pushed above the lowest CV state. This is most easily done by making the open-shell orbital t very compact, which naturally leads to a larger F _ttβ − F _ttα (due to a larger (tt|tt)) but a smaller F _ppβ − F _ppα , p ∈ {i, a} (due to a small absolute overlap between the p and t orbitals).

Now, it becomes evident that CuP fits the above design principles very well. The unpaired electron in the ground state of CuP is on the Cu 3 d x 2 − y 2 orbital (Figure 6), which is spatially localized. Moreover, the Cu 3 d x 2 − y 2 orbital occupies a different part of the molecule than the ligand π and π* orbitals, which results in a small absolute overlap between the orbitals and helps to reduce the effect of Eq. 9 on the CV excitation energies. To quantitatively assess the effect of Eq. 9 on the CO and OV excitation energies, we note that the X-TDDFT leading term Eq. 8 is nothing but the UKS orbital energy difference, if the shape differences of the UKS and ROKS orbitals are neglected. Therefore, we have plotted the UKS orbital energies of CuP in Figure 5 as well. Intriguingly, the α Cu 3 d x 2 − y 2 orbital now lies below the porphyrin π(a_1u ) and π(a_2u ) orbitals, while the β Cu 3 d x 2 − y 2 orbital lies above the porphyrin π*(e _g ) orbitals. Therefore, the differences of UKS orbital energies predict that the lowest excited states of CuP are the CV states obtained from exciting an electron from π(a_1u ) and π(a_2u ) to π*(e _g ). This is not only consistent with our U-TD-PBE0 excitation energies, but also the X-TD-PBE0 results (save for the ²S₂ state, which is higher than ²CT₁ computed at the respective levels of theory), despite that the latter method is spin-adapted. Note also that although the (tt|tt) integral leads to a huge splitting between the α and β Cu 3 d x 2 − y 2 orbitals, the splitting is only barely enough for the UKS orbital energy differences to predict a tripdoublet first excited state: if the β Cu 3 d x 2 − y 2 orbital were just 0.2 eV lower, one would predict that the CO-type CT excitation π(a_2u ) → Cu 3 d x 2 − y 2 is lower than the lowest tripdoublet π(a_2u ) → π*(e _g ). Alternatively, one may say that the HOMO-LUMO gap of the porphyrin ligand is barely narrow enough to fit within the energy window between the α and β Cu 3 d x 2 − y 2 orbitals, which clearly illustrates the importance of using a narrow-gap ligand for designing systems with a low-lying tripdoublet excited state.

Besides, we note that our SDSPT2 vertical absorption energy results suggest that three d-d excited states, ²dd₁, ²dd₂ and ²dd₃, lie barely below the ²T₁ state, although the actual energy ordering of ²dd₃ and ²T₁ cannot be said with full certainty. This means that CuP undergoes anti-Kasha fluorescence, which (as mentioned in the Introduction) is difficult to achieve especially for open-shell molecules. The ²dd states are very close in energy to the ²T₁ state, strongly favoring an ultrafast IC process, while the ²dd states themselves relax radiationlessly to the ground state on the hundred picosecond timescale, mainly through the cascade ²dd₂ →²dd₁ →²S₀ (Figure 4). In reality, however, our calculations suggest a ²T₁ →²dd₁ IC rate of only 2.00 × 10⁸s⁻¹ (with other ²T₁ →²dd IC processes even slower), which although competitive with the ²T₁ →²S₀ IC rate, is very small for an IC process with an adiabatic gap of only 0.42 eV (the estimated SDSPT2 adiabatic excitation energies of the ²dd₁, ²dd₂, ²dd₃ and ²T₁ being 1.66, 1.90, 2.06 and 2.08 eV, respectively). The reason is that the ²T₁ and ²dd states differ by a double excitation, while the nuclear derivative operator in the NACMEs is a one-electron operator, making the IC forbidden to first order. Moreover, to transform the ²T₁ state to any of the ²dd states, one has to perform an LMCT excitation as well as an MLCT excitation (from α π*(e _g ) to α Cu 3 d x 2 − y 2 , and from β Cu 3d _xy /3d _xz /3d _yz /3 d z 2 to β π(a_2u )), both of which involve occupied/virtual orbital pairs with poor overlap. Therefore, the drastically different natures of the ligand-centered ²T₁ state and the metal-centered ²dd states effectively prevented the quenching of the luminescence by the ²dd states, despite that the ²T₁ and ²dd states are extremely close in energy.

To conclude this section, we briefly note that even making the first doublet excited state a tripdoublet state still does not guarantee the realization of tripdoublet fluorescence. Two remaining potential obstacles are (1) the IC of the tripdoublet state to the ground state and (2) the ISC of the tripdoublet to the lowest quartet state (which is almost always lower than the lowest tripdoublet state owing to Hund’s rule). Both can be inhibited by making the molecule rigid, which is indeed satisfied by the porphyrin ligand in CuP. Alternatively, if the ISC from the first quartet state to the ground state is slow (as is the case of CuP, thanks to the spatial symmetry selection rules), and the gap between the first doublet and the first quartet is comparable to the thermal energy kT at the current temperature, then the quartet state can undergo RISC to regenerate the tripdoublet state, which can then fluoresce. This is well-known as the thermally activated delayed fluorescence (TADF) mechanism (Parker and Hatchard, 1961; Endo et al., 2011; Yang et al., 2017), although existing TADF molecules typically fluoresce from singlets and use a triplet “reservoir state” to achieve delayed fluorescence. In order for the TADF mechanism to outcompete the phosphorescence from the first quartet state, both the phosphorescence rate and the doublet-quartet gap have to be small. While the low phosphorescence rate of CuP can be explained by the fact that copper is a relatively light element, the small ²T₁–⁴T₁ gap of CuP can be attributed to the distributions of the frontier orbitals of CuP. Recall that the X-TDDFT gap between a tripdoublet excitation Eq. 3 and the associated quartet excitation Eq. 4 is exactly given by the X-RPA gap (Li and Liu, 2011), which is equal to 3 2 ( i t | i t ) + ( t a | t a ) . However, both of the two integrals are small for the ²T₁ and ⁴T₁ states of CuP, since the orbitals i and a reside on the ligand while t is localized near the metal atom (Figure 6). Such a clean spatial separation of the metal and ligand CMOs (despite the close proximity of the metal and the ligand) can further be attributed to the fact that the Cu 3 d x 2 − y 2 orbital has a different irrep than those of the ligand π and π* orbitals, preventing the delocalization of the open-shell orbital to the π system of the porphyrin ligand; while the Cu 3 d x 2 − y 2 orbital can still delocalize through the σ bonds of the ligand, the delocalization is of limited extent due to the rather local electronic structures of typical σ bonds (Figure 6). Incidentally, the only other class of tripdoublet-fluorescing metalloporphyrins that we are aware of, i.e., vanadium(IV) oxo porphyrin complexes (Ake and Gouterman, 1969; Gouterman et al., 1970), are characterized by a single unpaired electron in the 3d _xy orbital, whose mixing with the ligand π and π* orbitals is also hindered by symmetry mismatches. Whether this can be extended to a general strategy of designing molecules that fluoresce from tripdoublet states (or more generally, molecules that possess small doublet-quartet gaps) will be explored in the future. Finally, we briefly note that the design of doublet molecules with TADF and/or phosphorescence is also an interesting subject and deserves attention in its own right.