The change in electronic state from one spin multiplicity to another, known as intersystem crossing, occurs in molecules via the relativistic phenomenon of spin-orbit coupling. Current means of estimating intersystem crossing rates rely on the perturbative evaluation of spin-orbit coupling effects. This perturbative approach, valid in lighter atoms where spin-orbit coupling is weaker, is expected to break down for heavier elements where relativistic effects become dominant. Methods which incorporate spin-orbit effects variationally, such as the exact-two-component (X2C) method, will be necessary to treat this strong-coupling regime. We present a novel procedure which produces a diabatic basis of spin-pure electronic states coupled by spin-orbit terms, generated from fully variational relativistic calculations. This method is implemented within X2C using time-dependent density-functional theory and is compared to results from a perturbative relativistic study in the weak spin-orbit coupling regime. Additional calculations on a more strongly spin-orbit-coupled $[UO2Cl4]2\u2212$ complex further illustrate the strengths of this method. This procedure will be valuable in the estimation of intersystem crossing rates within strongly spin-coupled species.

## I. INTRODUCTION

Intersystem crossing (ISC) is the change in quantum state from one multiplicity to another. A formally spin-forbidden process, ISC nevertheless plays an important role in photochemistry.^{1} Due to exchange coupling, there are nearly always triplet excited states lower in energy than the lowest singlet excited state, and these states are by their spin state slow in decaying radiatively to the ground state. In systems where ISC is favorable, these low-lying and long-lasting triplet states can be efficiently populated, with direct consequences for the fluorescent and phosphorescent properties of a material. Accurate calculation of the ISC rates in a molecule, then, can be essential for understanding its photophysical properties.

With respect to spin eigenstates, intersystem crossing between different spin multiplicities occurs due to spin-orbit coupling (SOC), a fundamentally relativistic phenomenon.^{2} In addition to allowing for spin-forbidden transitions, relativistic effects include the contraction of core atomic orbitals in heavier elements, the splitting of degenerate multiplets, and many of the novel characteristics of transition metals, rare earth, and heavy elements.^{3,4} Modern electronic structure methods that include the important effects of spin-dependent quantum mechanics in chemical dynamics tend toward one of the two assumptions: they either treat spin coupling operators perturbatively via the spin-Hamiltonian or include the spin operators variationally via relativistic methods. While perturbative methods have their merits,^{5,6} particularly in the interpretation of ISC, they are limited in their description of the time-dependent nature of the spin-dependent many-electron wavefunction. This is especially true for phenomena driven by *strong* spin couplings, such as in heavy-metal complexes. A proper description of spin dynamics must arise from application of the time-dependent variational principle. However, within variational relativistic theory, interpretation and computation of ISC become challenging because spin is no longer a good quantum number. In this work, we present a novel procedure for calculating SOC matrix elements, which are needed for estimating ISC rates, within the framework of variational relativistic electronic structure theory.

## II. THEORY

### A. Quantum dynamics

To address intersystem crossing rates, it is necessary first to describe how systems may change quantum states. Within the Born-Oppenheimer approximation, the time-dependent electronic wavefunction $\Psi t,R(t)$ is a function of both time *t* and the time-dependent nuclear degrees of freedom **R**(*t*). This wavefunction may be expressed as

The set of electronic basis states *ψ*_{I}(**R**) will parametrically depend on the time-dependent **R**(*t*) but are otherwise assumed to be time-independent, and therefore, the time dependence of Ψ is located within the coefficients *c*_{I}(*t*). The **R**(*t*) dependence of the basis states *ψ*_{I}(**R**) will be assumed for the remainder of the text. After inserting Eq. (1) into the time-dependent electronic Schrödinger equation, the equations of motion for the expansion coefficients are given by

Here, $\sigma LK(t)=\u27e8\psi L|\u2202\u2202t|\psi K\u27e9$ is the nonadiabatic coupling between states |*ψ*_{L}⟩ and |*ψ*_{K}⟩, and $HLK(t)=\u27e8\psi L|\u0124elR(t)|\psi K\u27e9$ is the electronic Hamiltonian expressed in the same basis. Population transfer between states is driven by the off-diagonal elements of *H* and *σ* in Eq. (2).

To simplify the dynamics, studies frequently perform simulations in one of the two simplifying bases. In the adiabatic basis, the basis states |*ψ*_{I}⟩ are taken from the eigenstates of *Ĥ*_{el}(**R**) at each point **R**(*t*) such that *H*_{LK}(*t*) = *ϵ*_{K}*δ*_{LK} is diagonal, and only the nonadiabatic coupling matrix *σ*_{LK} retains off-diagonal elements. In the other simplifying case, the basis states are rotated to a diabatic basis |*ψ*_{K′}⟩ so as to minimize the nonadiabatic coupling term *σ*_{L′K′} ∼ 0, leaving then only a nondiagonal *H*_{L′K′}. As is difficult to diabatize states over all geometries for any but the smallest systems, an approximate diabatic basis is often chosen according to a variety of methods; for a review, see Ref. 7. Once a basis has been chosen, what remains is to sample the possible nuclear configurations $R$, for example by propagating them along classical trajectories.

### B. Spin-orbit coupling

When studying intersystem crossing, one is interested in the population transfer between electronic states of different spin multiplicities, most often between singlets and triplets.^{1,2} However, within nonrelativistic quantum mechanics, the electronic Hamiltonian does not contain spin-coupling terms capable of driving spin dynamics. Consequently, the integration of spin degrees of freedom over different spin states is zero, causing the nonadiabatic coupling *σ*_{LK} to vanish and preventing the population of one spin state from another.

When the spin-orbit coupling is relatively weak compared to the energy difference between spin states, perturbation theory can be used to introduce the ISC driving force into the nonrelativistic Hamiltonian. Among the most common is the Breit-Pauli spin-orbit operator,^{8–13} derived from the Dirac Hamiltonian, given as

where **s**_{i} is the spin vector for each electron. The first term contains single-electron SOC effects, while the second is a two-electron operator and is commonly neglected,^{12,13} parameterized into the one-electron term,^{14} or approximated using an atomic mean-field approach.^{15–17} After adding such an effective SO operator to the electronic Hamiltonian, states which previously had well-defined spins will couple together and no longer belong to a single multiplicity. In general, the Hamiltonian with the Breit-Pauli operator is unbounded from below, and naive eigenstate calculations will suffer from variational collapse.^{8,9} As a result, approximate eigenstates of the nonrelativistic Hamiltonian are typically calculated first, with the Breit-Pauli operator then applied perturbatively. If the Hamiltonian is again diagonalized in the resultant spin-perturbed basis, a set of adiabatic states naturally result.^{18–20} If the states instead remain spin-pure, they are often considered to be a nearly diabatic basis. This is because, as mentioned above, nonadiabatic coupling between spin-pure electronic states of different multiplicities is automatically zero. In contrast, nonadiabatic coupling between states of the same multiplicity, i.e. singlets, will not necessarily be zero. When singlet-triplet conversions are of particular interest, a first-order estimate of the ISC rate may be obtained by neglecting coupling between states of the same multiplicity.^{21–23} However, a full description of the intersystem crossing process will require the inclusion of all couplings: for example, efficient triplet-triplet internal conversion can increase the rate of ISC by preventing the reverse triplet-singlet conversion.

This approach, equivalent to a perturbation theory that is first order in spin-orbit coupling, is reasonable when treating lighter elements where SOC is weaker. Further down the periodic table, however, when SOC effects become increasingly strong, it is expected that approaches beyond first-order perturbation theory are necessary. In such cases, *variational* inclusion of relativistic effects is desirable and new diabatization techniques are needed for studies of ISC dynamics.

Within the kinetic-balance condition, we apply a spin separation technique so that the Dirac Hamiltonian can be written as^{8,9}

where *V* is the scalar potential, *c* is the speed of light, and *m* is the electron mass. Spin and orbital angular momenta are coupled through the ** σ** ·

**p**

*V*×

**p**term: the vector

**contains the Pauli spin matrices, and**

*σ***p**is the linear momentum operator. Variational solutions to this Hamiltonian will no longer have a well-defined spin multiplicity, and the eigenstates will have four components. The first term in Eq. (4) is the spin-free portion of the Dirac Hamiltonian, which contains all scalar relativistic effects, while the second term gives rise to spin-couplings.

### C. Diabatization

From the point of view of perturbation theory, the associated spin-diabatic states are those at zeroth order in terms of spin-coupling perturbation. That is, by using only the spin-free part of the Hamiltonian, variationally solving the Dirac equation will generate states that include scalar relativistic effects at the atomic orbital (AO) level but are still spin-pure, producing true singlets, triplets, etc. In the following discussion, {*ψ*_{L}} refer to variational solutions of the Dirac equation including all spin-couplings and ${\psi \u0303K}$ are those obtained using only the spin-free portion of the Hamiltonian.

Unlike the perturbative approaches outlined above, the states {*ψ*_{L}} and ${\psi \u0303K}$ are not simply related by the first-order spin-orbit coupling term, as such couplings are included in the infinite perturbation limit within the variational treatment. In order to calculate the couplings between states of different multiplicities, therefore, we search for a unitary transformation relating {*ψ*_{L}} and ${\psi \u0303K}$. The wavefunction overlap between the two sets of states is

which is non-Hermitian because {*ψ*_{L}} and ${\psi \u0303K}$ are mutually nonorthogonal. To allow for a unitary transformation between the two sets, we perform a Löwdin orthogonalization of the spin-coupled states {*ψ*_{L}} within the basis of spin-free states ${\psi \u0303K}$. That is, we find the matrix *U* that diagonalizes the Hermitian matrix **O**^{†}**O** such that

where ** o** is a diagonal matrix. The transformation

**T**,

^{24}

is subsequently a unitary transformation of the spin-coupled states, projecting them into the spin-free basis. Finally, taking the diagonal spin-coupled Hamiltonian *H*,

and rotating it into the spin-free basis

yields a diabatic Hamiltonian with off-diagonal couplings, whose eigenvalues are the same as those of spin-coupled **H**. In effect, this is a means of spin-purifying spin-orbit-coupled states.

### D. TDDFT implementation

Thus far, this procedure is entirely general for any method that can generate states from a variational relativistic electronic structure theory. However, the means of evaluating the wavefunction overlaps in Eq. (5) will be method-dependent. In this work, we will use the relativistic two-component time-dependent density-functional theory (TDDFT)^{25–27} method to illustrate the application of this approach to obtain the ISC coupling strength; for a recent review, see Ref. 28.

Relativistic electronic structure methods based on the Dirac equation employ a four-component wavefunction ansatz.^{8,9} For most chemically relevant studies, approximate two-component methods are usually of sufficient accuracy. The Douglas-Kroll-Hess (DKH)^{29–33} and exact-two-component (X2C)^{17,25,34–44} methods are popular approaches that decouple the four-component Hamiltonian and yield a reduced-dimension electronic two-component Hamiltonian, though they differ in the approximation used to decouple the large and small components of the relativistic wavefunction. Once an approximate decoupling is achieved, the two-component relativistic electronic Hamiltonian can be used with most electronic structure methods, provided a noncollinear spinor formalism is employed.^{45} It is of particular note that, unlike for effective SO operators like Breit-Pauli, X2C and DKH are bounded from below and do not suffer from variational collapse.^{8,9} Consequently, these two-component methods allow for the variational inclusion of relativistic effects at each stage of the calculation.

In this work, we use a direct atomic-orbital based X2C transformation with the torque-free spin-density approach.^{25,27,42} This transformation is only performed for the one-electron portion of the electronic Hamiltonian. To account partially for two-electron relativistic effects, the SO operator within X2C is scaled by a semiempirical fudge factor,^{14} reminiscent of the use of effective nuclear charges to account for nuclear screening of valence electrons by core electrons.

In the case of linear-response TDDFT, excited-state wavefunctions are not directly produced, but rather transition density matrices, complicating the evaluation of properties between excited states such as the wavefunction overlap in Eq. (5). Pseudowavefunctions, meant to reproduce approximately the linear-response eigenvectors, have been proposed,^{46,47} particularly in the context of evaluating nonadiabatic coupling vectors between excited states. For the latter method, it was shown that pseudowavefunctions produce coupling vectors in close agreement with those obtained from full quadratic response.^{48} If the Tamm-Dancoff approximation (TDA)^{49} is invoked [equivalent to configuration interaction singles (CIS) in the case of time-dependent Hartree-Fock], the auxiliary wavefunctions associated with excited states are expressed as linear combinations of Slater determinants singly excited from the DFT ground state,

where |Φ_{0}⟩ is the ground-state determinant, and $XiaK$ is the amplitude corresponding to the excitation from occupied orbital *i* to virtual orbital *a* for the state *ψ*_{K}. We invoke the TDA for two reasons. First, despite being slightly less accurate overall than full linear-response TDDFT, the TDA in nonrelativistic calculations has been shown to yield more accurate singlet-triplet energy gaps,^{50–52} which are vital for an accurate description of ISC. Second, in the absence of the so-called de-excitation amplitudes present in full linear-response TDDFT, the auxiliary wavefunctions are more intuitively constructed and are orthonormal, which greatly simplifies the procedure outlined here. In the case of TDA-TDDFT, the wavefunction overlap in Eq. (5) is therefore evaluated as

The expectation value $\u27e8\Phi \u03030|a\u0303i\u2020a\u0303a\xe2b\u2020\xe2j|\Phi 0\u27e9$ is the overlap between two Slater determinants expressed in nonorthogonal sets of MOs. Such overlaps are commonly found both within nonorthogonal configuration interaction (NOCI)^{53} and whenever evaluating overlaps between wavefunctions at different time steps or nuclear geometries,^{24} and several methods for evaluating them efficiently have been discussed.^{54–56} The specific method employed in this work is outlined in the Appendix.

It should be noted that a previous study^{57} evaluated the approximate overlap between excited states calculated using TDDFT at the scalar relativistic level and those with spin-orbit effects included at the zero-order regular approximation (ZORA)^{58} level of theory. This was done in order to estimate the composition of spin-orbit-coupled states in terms of spin-pure states, and in principle, that technique would also be capable of estimating SOC matrix elements in a similar fashion. However, as the previous work was principally interested in a qualitative comparison, it neglected rotations between occupied and virtual orbitals (in effect assuming that the two ground states were identical), an assumption justified only when SOC is very weak. In addition, X2C treats relativistic effects substantially more accurately than does ZORA, and as such, the method introduced here should yield more accurate couplings.

## III. COMPUTATIONAL DETAILS

All calculations were performed using a locally modified copy of the developmental version of GAUSSIAN16.^{59} The geometry of the [Ru(bpy)_{3}]^{2+} complex was optimized using the PBE functional,^{60} with the Ahlrichs double-*ζ* basis set with polarization on light atoms,^{61} and the LanL2DZ effective core potential and associated double-*ζ* basis set^{62} on Ru. X2C calculations used the same functional and basis set for the light atoms, with the relativistic all-electron Sapporo triple-*ζ* basis^{63} on Ru. The geometry of the $[UO2Cl4]2\u2212$ complex was optimized using the CAM-B3LYP functional,^{64} with the 6-31g(d) basis set on light atoms and the LanL2DZ effective core potential and associated double-*ζ* basis set on U. X2C calculations then used the cc-pVTZ^{65,66} basis set for light atoms and the relativistic SARC-DKH2 basis set^{67} for U. All nonrelativistic density functional theory calculations were performed on a pruned (75,302) grid (grid = UltraFine in GAUSSIAN16), while the convergence of relativistic calculations was improved by switching to a pruned (175,974) grid for light atoms and an unpruned (250,974) grid for heavy atoms (grid = SuperFine).

## IV. RESULTS AND DISCUSSION

To test the validity of the approach outlined above for generating spin-orbit-coupling matrix elements, we compare our method against a recent DFT study^{68} of intersystem crossing in a ruthenium trisbipyridine ([Ru(bpy)_{3}]^{2+}) complex (Fig. 1), shown experimentally to have very rapid ISC.^{69,70} In the previous work, SOC matrix elements were evaluated by perturbatively applying the spin-orbit operator within the zero-order regular approximation (ZORA)^{58} to spin-pure states that included scalar relativistic effects at the ZORA level of theory.

The energy levels of [Ru(bpy)_{3}]^{2+}, calculated using X2C, are plotted in Fig. 2. Due to the high symmetry of the molecule, there are several pairs of degenerate states in the spin-pure calculation. Additionally, within each set of triplets, the states are degenerate to within a millielectronvolt. In contrast, upon the inclusion of spin-orbit effects, there is a dramatic change in the energetics. Total energies shift by nearly a tenth of an electronvolt, and the degeneracy of some triplet states is broken by up to 30 meV. In this system, the relative energetic positions of the states make it easy to identify their predominant spin state; in general, however, once relativistic effects are included, it is no longer necessarily meaningful to speak of states as singlets or triplets, but rather they are linear superpositions of several multiplicities. Finally, Fig. 2 illustrates the results of using the diabatization scheme presented in this work. By rotating the spin-coupled states into the spin-pure basis, the degeneracy of triplet states is almost completely recovered: T_{4} has a splitting of 6 meV, while no other set of triplet states is split by more than 3 meV. Importantly, though, their relative energies are distinct from those of the initial spin-pure calculation; the states retain information about spin-orbit coupling effects, and as shown below, the resultant off-diagonal Hamiltonian elements contain the necessary information to evaluate intersystem crossing rates.

The SOC matrix elements for the six lowest-energy singlets and nine lowest-energy triplets are compared in Table I: for degenerate pairs of states such as S_{2} and S_{3}, the root-mean-square of both matrix elements is reported. Despite the SOC matrix elements spanning three orders of magnitude, we see close agreement between the two approaches. The mean absolute deviation for each singlet (or pair of singlets) is only a few tens of wavenumbers; over all 24 matrix elements, only one has a deviation greater than 80 cm^{−1}. Given that ZORA uses different functional forms for both scalar and relativistic effects than X2C and also that SOC is evaluated perturbatively rather than variationally, it is not obvious that the two methods should yield similar values, yet the predictions for SOC in [Ru(bpy)_{3}]^{2+} are in very close agreement. It appears likely, then, that this procedure is valid in the calculation of SOC matrix elements between singlets and triplets.

. | Weakest coupling . | Strongest coupling . | . | ||||
---|---|---|---|---|---|---|---|

Singlet . | Triplet . | X2C . | ZORA . | Triplet . | X2C . | ZORA . | MAD . |

S_{0} | T_{1} | 1 | 1 | T_{2}/T_{3} | 367 | 380 | 5 |

S_{1} | T_{1} | 0 | 0 | T_{5}/T_{6} | 626 | 626 | 15 |

S_{2}/S_{3} | T_{1} | 22 | 29 | T_{7}/T_{8} | 632 | 622 | 23 |

S_{4}/S_{5} | T_{4} | 17 | 23 | T_{1} | 572 | 583 | 16 |

. | Weakest coupling . | Strongest coupling . | . | ||||
---|---|---|---|---|---|---|---|

Singlet . | Triplet . | X2C . | ZORA . | Triplet . | X2C . | ZORA . | MAD . |

S_{0} | T_{1} | 1 | 1 | T_{2}/T_{3} | 367 | 380 | 5 |

S_{1} | T_{1} | 0 | 0 | T_{5}/T_{6} | 626 | 626 | 15 |

S_{2}/S_{3} | T_{1} | 22 | 29 | T_{7}/T_{8} | 632 | 622 | 23 |

S_{4}/S_{5} | T_{4} | 17 | 23 | T_{1} | 572 | 583 | 16 |

To illustrate the effects of this method when applied to a more strongly spin-orbit-coupled system where perturbational methods might be expected to break down, we extend this analysis to the uranyl tetrachloride anion, $[UO2Cl4]2\u2212$. The actinide center leads to very strong SOC effects in this complex, and previous theoretical studies of this molecule exist for comparison.^{71,72} Excitation energies from spin-orbit-free calculations are given in Table II. When compared to results from a multireference CASPT2 calculation,^{71} very close agreement is seen, with most states agreeing to within a few tens of millielectronvolts. More striking are the results with spin-orbit coupling included, shown in Table III. Perturbative application of SOC in the CASPT2 calculation drops the energy of the lowest-lying excited states by less than 0.2 eV, while TDDFT with variational SOC whether two-component or four-component^{72} presents a massive change of nearly 0.5 eV. Clearly, variational treatment of SOC has a much larger effect in this complex than in the ruthenium complex above.

. | Excitation energy (eV) . | |
---|---|---|

State . | X2C-TDDFT . | CASPT2 . |

a ^{3}B_{1g} | 2.81 | 2.77 |

a ^{3}B_{2g} | 2.90 | 2.87 |

a ^{3}E_{g} | 3.13 | 3.28 |

a ^{1}E_{g} | 3.58 | 3.71 |

a ^{1}B_{1g} | 3.82 | 3.80 |

a ^{1}B_{2g} | 3.89 | 3.90 |

. | Excitation energy (eV) . | |
---|---|---|

State . | X2C-TDDFT . | CASPT2 . |

a ^{3}B_{1g} | 2.81 | 2.77 |

a ^{3}B_{2g} | 2.90 | 2.87 |

a ^{3}E_{g} | 3.13 | 3.28 |

a ^{1}E_{g} | 3.58 | 3.71 |

a ^{1}B_{1g} | 3.82 | 3.80 |

a ^{1}B_{2g} | 3.89 | 3.90 |

. | Excitation energy (eV) . | ||
---|---|---|---|

State . | X2C-TDDFT . | 4c-TDDFT
. | CASPT2 . |

a E_{g} | 2.36 | 2.38 | 2.61 |

a B_{2g} | 2.41 | 2.48 | 2.64 |

a B_{1g} | 2.36 | 2.41 | 2.74 |

b E_{g} | 2.54 | 2.62 | 2.83 |

b B_{1g} | 2.67 | 2.70 | 2.98 |

b B_{2g} | 2.67 | 2.71 | 3.02 |

c E_{g} | 3.16 | 3.19 | 3.41 |

a A_{2g} | 3.32 | 3.38 | 3.70 |

a A_{1g} | 3.33 | 3.38 | 3.70 |

d E_{g} | 3.59 | 3.69 | 3.96 |

c B_{1g} | 3.79 | 3.79 | 3.97 |

c B_{2g} | 3.82 | 3.82 | 4.03 |

. | Excitation energy (eV) . | ||
---|---|---|---|

State . | X2C-TDDFT . | 4c-TDDFT
. | CASPT2 . |

a E_{g} | 2.36 | 2.38 | 2.61 |

a B_{2g} | 2.41 | 2.48 | 2.64 |

a B_{1g} | 2.36 | 2.41 | 2.74 |

b E_{g} | 2.54 | 2.62 | 2.83 |

b B_{1g} | 2.67 | 2.70 | 2.98 |

b B_{2g} | 2.67 | 2.71 | 3.02 |

c E_{g} | 3.16 | 3.19 | 3.41 |

a A_{2g} | 3.32 | 3.38 | 3.70 |

a A_{1g} | 3.33 | 3.38 | 3.70 |

d E_{g} | 3.59 | 3.69 | 3.96 |

c B_{1g} | 3.79 | 3.79 | 3.97 |

c B_{2g} | 3.82 | 3.82 | 4.03 |

The strengths of this method may be more clearly seen in the energy-level diagram of the $[UO2Cl4]2\u2212$ anion, shown in Fig. 3. Spin-orbit coupling leads to an enormous energy splitting of nearly half an eV, and the identification of states as singlets or triplets is clearly impossible. However, despite very powerful spin-orbit splitting, this method still recovers near-degeneracy of triplet states, with correspondingly large SOC matrix elements preserving the splitting of eigenvalues. The question remains of why exact degeneracy is not recovered. One possibility is that more excited states must be included in the sample space in order to adequately express spin-pure states in the spin-orbit-coupled basis, which is explored in Table IV. When only the first 12 states (primarily arising from the four lowest-energy triplets) are included, degeneracy is quite poor, with each given set of triplets possessing a standard deviation in energy of over 50 meV. However, including even the four next lowest-energy states (corresponding to the four lowest-energy singlets) dramatically improves the triplet degeneracies, and by including 35 states, degeneracy is recovered to the millielectronvolt level. Full degeneracy is still not completely recovered, demonstrating the incomplete nature of this transformation, but sufficient degeneracy is recovered that we believe this method will be applicable even to strongly coupled systems.

## V. CONCLUSION

In this work, we have taken an important step toward the estimation of ISC rates using variational relativistic theory. While some perturbative treatments of spin-orbit coupling using scalar-relativistic methods have been proposed, for some of which^{18,73,74} analytical nuclear gradients have been implemented, the true strength of relativistic methods is their ability to include spin-orbit coupling effects fully variationally and at the molecular orbital level. However, variational treatment of ISC has not previously been implemented because (1) in the absence of a perturbative SO operator, it is nontrivial to evaluate couplings in a “spin-diabatic” basis and (2) the analytical excited-state nuclear gradients needed to evaluate nonadiabatic coupling^{75} are not yet available for two-component methods. The procedure outlined here, in contrast, takes states with spin-orbit effects included variationally and rotates them into a spin-pure basis. This effective spin-purification produces states of known multiplicities, making their identification easier, while also coupling them together in a manner that allows for intersystem crossing.

It should be noted that, in a finite atomic-centered basis and including only a subset of all possible excitations, the spin-pure states do not form a complete basis for the spin-coupled states, and care must be taken to ensure that not too much information is discarded in the change in the basis. In principle, the completeness of the two bases could be improved by increasing both the number of atomic basis functions and the number of excited states calculated, as shown in the uranyl tetrachloride example above. However, if the calculation is restricted to excited states with a finite number of excitations (e.g., singles, singles and doubles, etc.), the two sets of states may still not fully span one another, as they are comprised of different molecular orbitals. In practice, however, we have found that the ground state determinants overlap by at least 0.85 (0.99 in the case of the Ru complex above), allowing for the use of this method.

While SOC matrix elements are necessary in the treatment of intersystem crossing, more work is needed before variational relativistic theory can be used to estimate ISC rates. Needed next is the inclusion of nuclear dynamics effects. At present, in the absence of excited-state nuclear gradients, dynamics will need to be performed using a lower level of theory, a possibility explored in an upcoming publication. The development of such nuclear gradients would be invaluable to several extensions of relativistic theory: geometry optimizations, nuclear dynamics of spin-coupled electronic states, and of course the explicit evaluation of nonadiabatic coupling. Even without such gradients, however, the machinery for calculating wavefunction overlaps in the manner introduced here could be used to probe coupling in the adiabatic basis by the numerical evaluation of $\u27e8\psi K|\u2202\u2202t|\psi L\u27e9$,^{76,77} or by implementing the local diabatization method,^{24,78} where states at one time step are rotated into the basis of states at the preceding time step. Also to be determined is whether spin-pure states must be evaluated at each time step of a trajectory or whether a single representative set (or at least a smaller number of sets) of spin-pure states could provide a sufficiently accurate diabatic basis. By combining the evaluation of spin-orbit coupling using variational relativistic theory with a method for incorporating nuclear dynamics, we will be able to directly model population transfer between states of different multiplicities after photoexcitation, allowing us to accurately and efficiently estimate intersystem crossing rates.

## SUPPLEMENTARY MATERIAL

The full set of spin-orbit coupling matrix elements for the ruthenium complex calculated using X2C-TDDFT and their comparison to those from ZORA^{68} are presented in the supplementary material.

## ACKNOWLEDGMENTS

The development of the two-component electronic structure method is funded by the U.S. Department of Energy (Grant No. DE-SC0006863). The development of nonperturbative spectroscopic methods is supported by the National Science Foundation (Grant No. CHE-1856210). Computational intersystem crossing is partially supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences, through Argonne National Laboratory under Contract No. DE-AC02-06CH11357. Computations were facilitated through the use of advanced computational, storage, and networking infrastructure provided by the Hyak supercomputer system at the University of Washington, funded by the Student Technology Fee and the National Science Foundation (Grant No. MRI-1624430).

### APPENDIX: WAVEFUNCTION OVERLAP

For two sets of nonorthogonal molecular orbitals (MOs) defined by the matrices $C\u0303$ and **C**, constructed from the same atomic orbital (AO) basis, the overlap between the MOs $|\varphi \u0303p$ and |*ϕ*_{q}⟩ is given by

where **S** is the AO overlap matrix. The wavefunction overlap between two *n*-electron Slater determinants $|\Phi \u0303$ and |Φ⟩ constructed from these MOs is given by the determinant

Naively implemented, the overlaps in Eq. (11) would require calculating the overlap $\u27e8\Phi \u0303ia|\Phi jb\u27e9$ for all *n*^{2}$v2$ possibilities $i,j,a,b$, where $|\Phi ia$ is the Slater determinant formed by promoting an electron from occupied orbital *i* to virtual orbital *a*. However, a substantial reduction in complexity can be realized by noting that rotations among occupied orbitals do not change the energy of the DFT ground state. One is free, then, to rotate the occupied orbitals $C\u0303O$ and **C**_{O}, respectively, by the unitary transformations **U** and **V**, found by performing a singular value decomposition (SVD) of the ground-state overlap

where **s** is a diagonal matrix. The overlap between the respective ground state determinants $\u27e8\Phi \u03030|\Phi 0\u27e9$ is now (to within a phase) equal to the product of the singular values ∏*s*_{i}. The new occupied orbitals, given as $C\u0303O\u2032=C\u0303OU$ and $CO\u2032$ = **C**_{O}**V**, in the two sets are then biorthonormal: that is, $\u27e8\varphi i\u0303|\varphi j\u27e9=\delta ij$. By rotating **M** into the new basis,

as well as the TDA-TDDFT excitation coefficients,

the overlap $\u27e8\Phi \u0303ia|\Phi jb\u27e9$ can now much more simply be evaluated as

reducing the complexity of the overall overlap calculation from $O(n5v2)$ to $O(n2v2)$.

## REFERENCES

_{2,3}-edge x-ray absorption spectroscopy with real-time exact two-component relativistic time-dependent density functional theory

_{3}]

^{2+}

_{2}Cl

_{4}]

^{2−}