Intersystem Crossing and Exciton-Defect Coupling of Spin Defects in Hexagonal Boron Nitride

Despite the recognition of two-dimensional (2D) systems as emerging and scalable host materials of single photon emitters or spin qubits, uncontrolled and undetermined chemical nature of these quantum defects has been a roadblock to further development. Leveraging the design of extrinsic defects can circumvent these persistent issues and provide an ultimate solution. Here we established a complete theoretical framework to accurately and systematically design new quantum defects in wide-bandgap 2D systems. With this approach, essential static and dynamical properties are equally considered for qubit discovery. In particular, many-body interactions such as defect-exciton couplings are vital for describing excited state properties of defects in ultrathin 2D systems. Meanwhile, nonradiative processes such as phonon-assisted decay and intersystem crossing rates require careful evaluation, which compete together with radiative processes. From a thorough screening of defects based on first-principles calculations, we identified the Ti-vacancy complex as a promising defect in hexagonal boron nitride for spin qubits, with a triplet ground state, large zero-field splitting, and a prominent intersystem crossing rate highly desirable for spin-state initialization and qubit operation.


INTRODUCTION
Optically addressable defect-based qubits offer a distinct advantage in their ability to operate with high fidelity under room temperature conditions. 1,2 Despite tremendous progress made in years of research, the existing systems remain inadequate for real-world applications. The identification of stable single photon emitters in 2D materials has opened up a new playground for novel quantum phenomena and quantum technology applications, with improved scalability in device fabrication and a leverage in doping spatial control, qubit entanglement, and qubit tuning. 3,4 In particular, hexagonal boron nitride (h-BN) has demonstrated that it can host stable defect-based single photon emitters (SPEs) [5][6][7][8] and spin triplet defects. 9,10 However, persistent challenges must be resolved before 2D quantum defects can become the most promising quantum information platform. These challenges include the undetermined chemical nature of existing SPEs, 7,11 difficulties in controlled generation of desired spin defects, and scarcity of reliable theoretical methods which can accurately predict critical physical parameters for defects in 2D materials due to their complex many-body interactions.
To circumvent these challenges, design of promising spin defects by high-integrity theoretical methods is urgently needed. Introducing extrinsic defects can be unambiguously produced and controlled, which fundamentally solves the current issues of undetermined chemical nature of existing SPEs in 2D systems. As highlighted by Ref. 12,13, promising spin qubit candidates should satisfy several essential criteria: deep defect levels, stable high spin states, large zero-field splitting, efficient radiative recombination, high intersystem cross-ing rates and long spin coherence and relaxation time. Using these criteria for theoretical screening can effectively identify promising candidates but requires theoretical development of first-principles methods, significantly beyond the static and mean-field level. For example, accurate defect charge transition levels in 2D materials necessitates careful treatment of defect charge corrections for removal of spurious charge interactions [14][15][16] and electron correlations for non-neutral excitation, e.g. from GW approximations [16][17][18] or Koopmans-compliant hybrid functionals. [19][20][21][22] Optical excitation and exciton radiative lifetime must account for defect-exciton interactions, e.g. by solving the Bethe-Salpeter equation, due to large exciton binding energies in 2D systems. [23][24][25] Spin-phonon relaxation time calls for a general theoretical approach to treat complex symmetry and state degeneracy of defective systems, along the line of recent development based on ab-initio density matrix approach. 26 Spin coherence time due to the nuclei spin and electron spin coupling can be accurately predicted for defects in solids by combining first-principles and spin Hamiltonian approaches. 27,28 In the end, nonradiative processes, such as phonon-assisted nonradiative recombination have been recently computed with first-principles electron-phonon couplings for defects in h-BN, 29 and resulted in less competitive rates than corresponding radiative processes. However, the spin-orbit induced intersystem crossing as the key process for pure spin state initialization during qubit operation has not been investigated for spin defects in 2D materials from first-principles. This work has developed a complete theoretical framework which enables the design of spin defects based on the critical physical parameters mentioned above and highlighted in Figure 1a. We employed state-of-the-art first-principles methods, focusing on many-body interaction such as defect-exciton couplings and dynamical processes through radiative and nonradiative recombinations. We developed methodology to compute nonradiative intersystem crossing rates with explicit overlap of phonon wavefunctions beyond current implementations in the Huang-Rhys approximation. 30 We showcase the discovery of transition metal complexes such as Ti and Mo with vacancy to be optically addressable spin triplet defects in h-BN. They are the first extrinsic defects in h-BN to be predicted as stable triplets with large zero-field splitting and spin-selective decay, which will set 2D quantum defects at a competitive stage with NV center in diamond for quantum technology applications. Figure 1. Doping at the divacancy site of h-BN. a Schematic of the screening criteria and workflow developed in this work, where we first search for defects with stable triplet ground state, followed by large zero-field splitting (ZFS), then "bright" optical transitions between defect states required for SPEs or qubit operation by photon, and at the end large intersystem crossing rate (ISC) critical for pure spin state initialization. b Divacancy site in h-BN corresponding to adjacent B and N vacancies (denoted by V B and V N ). c Top-view and d Side-view of a typical doping configuration when placed at the divacancy site, denoted by X VV . Atoms are distinguished by color: grey=N, green=B, purple=Mo, blue=Ti, red=X (a generic dopant).

RESULTS AND DISCUSSION
In the development of spin qubits in 3D systems (e.g. diamond, SiC, and AlN) defects beyond sp dangling bonds from N or C have been explored. In particular, large metal ions plus anion vacancy in AlN and SiC were found to have potential as qubits due to triplet ground states and large zero-field splitting (ZFS). 31 Similar defects may be explored in 2D materials, 32 such as the systems shown in Figure 1b-d. This opens up the possibility of overcoming the current limitations of uncontrolled and undetermined chemical nature of 2D defects, and unsatisfactory spin dependent properties of existing defects. In the following, we will start the computational screening of spin defects with static properties of the ground state (spin state, defect formation energy and ZFS) and the excited state (optical spectra), then we will discuss dynamical properties including radiative and nonradiative (phonon-assisted spin conserving and spin flip) processes, as the flow chart shown in Figure 1a. We will summarize the complete defect discovery procedure and discuss the outlook at the end.

Screening Triplet Spin Defects in h-BN
To identify stable qubits in h-BN, we start from screening neutral dopant-vacancy defects for a triplet ground state based on total energy calculations of different spin states at both semilocal PBE (PerdewBurke-Ernzerhof) and hybrid functional levels. We considered the dopant substitution at a divacancy site in h-BN ( Figure 1b)  We further confirmed the thermodynamic charge stability of these defect candidates via calculations of defect formation energy and charge transition levels. As shown in SI Figure   S1, both Ti VV and Mo VV defects have a stable neutral (q = 0) region for a large range of Fermi level (ε F ), from 2.2 eV to 5.6 eV for Mo VV and from 2.9 eV to 6.1 eV for Ti VV . These neutral states will be stable in intrinsic h-BN systems or with weak p-type or n-type doping.

Zero-Field Splitting and Optically Detected Magnetic Resonance
With a confirmed triplet ground state, we next computed the two defects' zero-field splitting.
A large ZFS is necessary to isolate the m s = ±1 and m s = 0 levels even at zero magnetic field allowing for controllable preparation of the spin qubit. Here we computed the first-order This implies there will be a splitting of the m s = +1 and m s = −1 sublevels even at zero magnetic field. The splitting of these sublevels is typically measured from optically detected magnetic resonance (ODMR) experiments and the simulated resonant frequencies are given in SI Figure S2, which can be used as reference for experimental defect identification. The separation of secondary spin sublevels of the Ti VV and Mo VV defect can be leveraged to potentially initiate well-defined quantum states for implementation of spin-based qubits.

Single-Particle Level, Optical Spectra and Radiative Lifetime
Deep defect levels with an allowed isolated and spin-polarized optical transition is another key requirement for potential defects as single photon emitters and spin qubits. Therefore, we will next investigate the electronic structure and optical spectra of the Ti VV and Mo VV spin defects from density functional theory (DFT) and many-body perturbation theory.
We computed the single particle energy levels of Ti VV and Mo VV at three different levels of theory: PBE, PBE0(α) and G 0 W 0 @PBE0(α). Using a hybrid functional as the startingpoint for GW calculations here is advantageous, because hybrid functionals can correct the self-interaction errors for 3d transition metal defects compared to semi-local functional, therefore providing an improved starting point for GW quasiparticle energies. 38,39 For Ti VV defect in h-BN, both wavefunction distribution and ordering of defect states are different between PBE and PBE0(α) as shown in the right panels of Figure 4. However, for sp defects such as N B V N in h-BN, PBE and PBE0(α) have similar wavefunction character and ordering, 29 so G 0 W 0 @PBE may be sufficient.
The single-particle energy levels of Ti VV defect in h-BN are shown in Figure 2. Those of the Mo VV defect can be found in SI Figure S4. The single-particle diagram of Ti VV shows there are several defect states 1a ↑ , 1a ↑ , 2a ↑ and 1a ↓ in the gap of h-BN. Furthermore, the band gap and energy difference between defect states are narrower at PBE, compared with the PBE0(α) and G 0 W 0 @PBE0(α) levels. The band gaps of h-BN between the latter two methods are rather similar, and the separation between occupied and unoccupied defect levels is 0.23 eV larger in G 0 W 0 @PBE0(α) than PBE0(α). We then investigated the optical properties to directly identify bright spin-polarized optical transitions between defect states. We compared the optical excitations with three levels of theory: Random Phase Approximation (RPA) with PBE states, RPA with G 0 W 0 @PBE0(α) states, and solving the Bethe-Salpeter Equation (BSE) with G 0 W 0 @PBE0(α) states, as displayed in Table 1. The difference between RPA and BSE is that the latter includes excitonic (e-h) effects. [40][41][42][43] In Figure 3, we show the BSE absorption spectra for both pristine (dashed lines) and defective h-BN (solid lines). More computational details can be found in method-ology section and SI Figure S3 and Table S3.
We found for the bulk-state transitions above 6 eV, the peak positions are rather similar between pristine and defective systems, consistent with the expectation that these defects are rather localized with limited influence on bulk states. For the defect-state-involved transitions in the left panels of Figure 3 and Table 1, we focus on the first bright transition i.e. 1a ↑ → 2a ↑ for Ti VV and 1a ↑ → 3a ↑ for Mo VV , energetically located within the infrared (IR) region. The first exciton binding energy was computed from the difference of the excitation energies at G 0 W 0 + RPA@PBE0(α) and G 0 W 0 + BSE@PBE0(α) in Table 1. We obtain exciton binding energies of 4.02 eV and 3.97 eV for the Ti VV and Mo VV defects, respectively, much larger than those from pristine monolayer h-BN (2.10 eV) and N B V N defect (2.43 eV).
Interestingly, the stronger exciton binding energy for transition metal defects is consistent with the prominent localization of their exciton wavefunctions (shown in Figure 3 left panels), in sharp contrast to the excitons from pristine h-BN and sp defect (shown in SI Figure S6).
The strong localization of defect-bound exciton also leads to weaker oscillator strength 44 (µ 2 e−h ), as shown in Table 1, and a stronger red shift of defect peaks compared to the bulkstate peaks in the G 0 W 0 + BSE@PBE0(α) spectra relative to the G 0 W 0 + RPA@PBE0(α) spectra (SI Figure S5).
The radiative lifetime for the first defect-defect recombination is computed based on the approach discussed in Ref. 45 (see method section for details). Lifetimes were obtained by using the modulus square of the transition dipole moments (µ 2 e−h ) and excitation energy (E 0 ) at the long wavelength limit from several levels of theory as discussed above, and shown in Table 1. As a comparison, we calculated the excitation energy and radiative lifetime of     An overview of the multiplet structure and excited-state dynamics is given in Figure 4 and briefly summarized here. Beginning from a spin-conserved optical excitation from the triplet ground state 3 0 A to the triplet excited state 3 1 A , the excited state relaxation and recombination can go through several pathways. The excited state can directly return to the ground state via a radiative or nonradiative process (red lines). Or the system may relax to another excited state with lower symmetry through a pseudo Jahn-Teller distortion (PJT; dark blue lines), and ultimately recombine back to the ground state nonradiatively.
A third pathway is to nonradiatively relax to an intermediate singlet excited state through a spin-flip intersystem crossing (ISC), and then again recombine back to the ground state (dashed light-blue lines). In particular, this ISC pathway is critical for the preparation of a pure spin state, similar to NV center in diamond. We will compute the rate for each radiative or nonradiative process, and determine the most competitive pathway under the operation condition.

Direct Radiative and Nonradiative Recombination
First, we will consider the direct recombination process for the Ti VV defect. Figure 5 shows the configuration diagram of the ground state triplet configuration 3 0 A of Ti VV alongside the first triplet excited state 3 1 A . The zero-phonon line (ZPL) of the Ti VV defect for direct recombination is calculated to be 0.49 eV based on the total energy difference of ground and excited states by constrained occupation DFT (CDFT) at PBE shown in Figure 5. ZPL can be also computed by a two-step process: vertical excitation energy from the ground state subtracting the relaxation energy to equilibrium geometry at the excited state (∆E in Figure 5). Here with the vertical excitation energy computed by BSE as 0.556 eV in Table 1, subtracting relaxation energy 26 meV, we obtain 0.53 eV for ZPL, in reasonable agreement with the CDFT total energy calculation.
As discussed earlier and shown in Table 1, the radiative lifetime of direct recombination at Ti VV defect is rather long compared to sp defects due to the strongly localized nature of its exciton wavefunctions. On the other hand, the Huang-Rhys factor (S f ) is calculated to be 0.91 for this recombination process, which implies extremely small electron-phonon coupling and potentially an even slower nonradiative process. Following the formalism presented in

Nonradiative Recombination through the PJT Triplet Excited State
Although the direct recombination process is slow for Ti VV defect, we found that the triplet excited state of Ti VV can undergo a PJT, reducing its symmetry from C S to C 1 (analogous to distortions for the N B V N defect 51,52 ). This C 1 geometry configuration has a 3 1 A multiplet state. It can recombine with ground state 3 0 A nonradiatively along the configuration diagram in Figure 5, with a reduced ZPL 0.482 eV compared to direct recombination. Notably, the barrier at the potential energy surface crossing between 3 1 A and 3 0 A is low and the Huang-Rhys factor is significantly larger S f = 14.95. We again computed the corresponding nonradiative recombination rate following Eq. 11. Consistent with the value of S f , the phonon term X if was several orders of magnitude larger than the one of direct recombination. Therefore, the resulting nonradiative recombination is extremely fast with a lifetime of 10 −2 ps at 10 K (dark blue dashed lines in Figure 4). Figure 6. The pseudo Jahn-Teller effect in the triplet excited state of Ti VV defect in h-BN. a Configurational energy diagram showing potential distortion to a lower energy, lower symmetry C 1 geometry from the C S geometry. b Distortion due to the pseudo Jahn-Teller effect with the changes in Ti-B/Ti-N bond lengths provided inÅ (green=B, grey=N, blue=Ti).
As such, it is possible that upon optical excitation, the recombination via PJT may play a crucial role. We approximated the rate of PJT between two triplet excited states through a classical rate equation, Γ C N R = νe −βEact , 53,54 as we found the transition between 3 1 A and 3 1 A is purely structural without electron transfer. Here, ν is approximated by the effective phonon frequency hν = 13.30 meV. The classical activation barrier E act is computed by linear extrapolation between the C S and C 1 configurations, as shown in Figure 6 with the Jahn-Teller barrier δ JT = 0.006 eV. We then obtain the C S to C 1 transition with a lifetime of 10 2 ps at 10 K. Combining with the extremely fast recombination between 3 1 A and 3 0 A , the indirect nonradiative process through PJT is significantly faster than direct recombination

Spin-Orbit Coupling and Nonradiative Intersystem Crossing Rate
Lastly, we considered the possibility of an ISC between the triplet excited state To compute the ISC rate, we considered a new approach which is a derivative of the nonradiative recombination formalism presented in Eq. 11:  Table 2 and light blue lines in Figure 4.
The results of all the nonradiative pathways are summarized in Table 2 and are displayed in Figure 4 along with the radiative pathway. In short, the spin conserved optical excita-

CONCLUSIONS
In summary, we proposed a general theoretical framework for identifying and designing optically-addressable spin defects for the future development of quantum emitter and quantum qubit systems. We start from searching for defects with triplet ground state by DFT total energy calculations which allow for rapid identification of possible candidates. Here we found that the Ti VV and Mo VV defects in h-BN have a neutral triplet ground state. We then computed zero-field splitting of secondary spin quantum sublevels and found they are sizable for both defects, larger than that of NV center in diamond, enabling possible control of these levels for qubit operation. Next the electronic structure and optical spectra of each defect were computed from many-body perturbation theory. Our results show optically addressable transitions for these defects, with much larger exciton binding energy than intrinsic sp defects, due to strongly localized exction wavefunctions. Finally, we analyzed all possible radiative and nonradiative dynamical processes with first-principles rate calculations. In particular, we identified a dominant spin-selective decay pathway via intersystem crossing at the Ti VV defect in h-BN, an order of magnitude faster than the competitive pseudo Jahn-Teller induced nonradiative recombination, and shows a key advantage for initial pure spin state preparation and qubit operation.
This work emphasizes that the theoretical discovery of new spin defects requires careful treatment of many-body interactions and various radiative and nonradiative dynamical processes such as intersystem crossing. We demonstrate high potential of extrinsic spin defects in 2D host materials as qubits for quantum information science. Future work will involve further examination of spin coherence time and its dominant decoherence mechanism, as well as other spectroscopic fingerprints from first-principles calculations to facilitate experimental validation of these defects.

Thermodynamic Charge Transition Levels and Defect Formation Energy
The defect formation energy (F E q ) was computed for the Ti VV and Mo VV defects following: where E q is the total energy of the defect system with charge q, E pst is the total energy of the

Zero-Field Splitting
The first-order ZFS due to spin-spin interactions was computed for the dipole-dipole interactions of the electron spin: Here, µ 0 is the magnetic permeability of vacuum, g e is the electron gyromagnetic ratio, is the Planck's constant, s 1 , s 2 is the spin of first and second electron, respectively, and r is the displacement vector between these two electron. The spatial and spin dependence can be separated by introducing the effective total spin S = i s i . This yields a Hamiltonian of the form H ss = S TD S, which introduces the traceless zero-field splitting tensorD. It is common to consider the axial and rhombic ZFS parameters D and E which can be acquired from theD tensor: Following the formalism of Rayson et al., 33 the ZFS tensorD can be computed with periodic boundary conditions as: Here the summation on pairs of i, j runs over all occupied spin-up and spin-down states, with χ ij taking the value +1 for parallel spin and −1 for anti-parallel spin, and Ψ ij (r 1 , r 2 ) is a two-particle Slater determinant constructed from the Kohn-Sham wavefunctions of the ith and jth states. This procedure was implemented as a post-processing code interfaced with Quantum ESPRESSO. To verify our implementation is accurate, we computed the ZFS of the NV center in diamond which has a well-established result. Using ONCV pseudopotentials, we obtained a ZFS of 3.0 GHz for NV center, in perfect agreement with previous reported results. 31

Radiative Recombination
In order to quantitatively study radiative processes, we computed the radiative rate Γ R from Fermi's Golden Rule and considered the excitonic effects by solving BSE: 45 Here, the radiative recombination rate is computed between the ground state G and the two-particle excited state S(Q ex ), 1 q L ,λ and 0 denote the presence and absence of a photon, H R is the electron-photon coupling (electromagnetic) Hamiltonian, E(Q ex ) is the exciton energy, and c is the speed of light. The summation indices in Eq. 8 run over all possible wavevector (q L ) and polarization (λ) of the photon. Following the approach described in Ref. 45, the radiative rate (inverse of radiative lifetime τ R ) in SI unit at zero temperature can be computed for isolated defect-defect transitions as: where e is the charge of an electron, 0 is vacuum permittivity, E 0 is the exciton energy at Q ex = 0, n D is the reflective index of the host material and µ 2 e−h is the modulus square of exciton dipole moment with length 2 unit. Note that Eq. 9 considers defect-defect transitions in the dilute limit; therefore the lifetime formula for zero-dimensional systems embedded in a host material is used 8,54,64 (also considering n D is unity in isolated 2D systems at the long-wavelength limit). We did not consider the radiative lifetime of Ti VV defect at a finite temperature because the first and second excitation energy separation is much larger than kT . Therefore a thermal average of the first and higher excited states is not necessary and the first excited state radiative lifetime is nearly the same at 10 K as zero temperature.

Phonon-Assisted Nonradiative Recombination
In this work, we compute the phonon-assisted nonradiative recombination rate via a Fermi's golden rule approach: Here, Γ N R is the nonradiative recombination rate between electron state i in phonon state n and electron state f in phonon state m, p in is the thermal probability distribution of the initial state |in , H e−ph is the electron-phonon coupling Hamiltonian, g is the degeneracy factor and E in is the energy of vibronic state |in . Within the static coupling and onedimensional (1D) effective phonon approximations, the nonradiative recombination can be reduced to: Here, the static coupling approximation naturally separates the nonradiative recombination rate into phonon and electronic terms, X if and W if , respectively. The 1D phonon approximation introduces a generalized coordinate Q, with effective frequency ω i and ω f . The phonon overlap in Eq. 12 can be computed using the quantum harmonic oscillator wavefunctions with Q − Q a from the configuration diagram ( Figure 5). Meanwhile the electronic overlap in Eq. 13 is computed by finite difference using the Kohn-Sham orbitals from DFT at the Γ point. The nonradiative lifetime τ N R is given by taking the inverse of rate Γ N R .
Supercell convergence of phonon-assisted nonradiative lifetime is shown in SI Table S4. We validated the 1D effective phonon approximation by comparing the Huang-Rhys factor with the full phonon calculations in SI Table S5.

Spin-Orbit Coupling Constant
Spin-orbit coupling (SOC) can entangle triplet and singlet states yielding the possibility for a spin-flip transition. The SOC operator is given to zero-order by: 65 where c is the speed of light, m e is the mass of an electron, p and S are the momentum and spin of electron i and V is the nuclear potential energy. The spin-orbit interaction can be rewritten in terms of the angular momentum L and the SOC strength λ as, 65 H so = i λ ⊥ (L x,i S x,i + L y,i S y,i ) + λ z L z,i S z,i .
where λ ⊥ and λ z denote the non-axial and axial SOC strength, respectively. The SOC strength was computed for the Ti VV defect in h-BN using the ORCA code by TD-DFT. 66,67 More computational details can be found in the SI.

DATA AVAILABILITY
The data that support the findings of this study and the code for the first-principles methods proposed in this study are available from the corresponding author (Yuan Ping) upon reasonable request. wavefunctions; 1D vs full-phonon HR factor comparison; details of spin-orbit coupling constants.