Ab initio theory of the negatively charged boron vacancy qubit in hexagonal boron nitride

Highly correlated orbitals coupled with phonons in two-dimension are identified for paramagnetic and optically active boron vacancy in hexagonal boron nitride by first principles methods which are responsible for recently observed optically detected magnetic resonance signal. Here, we report ab initio analysis of the correlated electronic structure of this center by density matrix renormalization group and Kohn-Sham density functional theory methods. By establishing the nature of the bright and dark states as well as the position of the energy levels, we provide a complete description of the magneto-optical properties and corresponding radiative and non-radiative routes which are responsible for the optical spin polarization and spin dependent luminescence of the defect. Our findings pave the way toward advancing the identification and characterization of room temperature quantum bits in two-dimensional solids.


INTRODUCTION
Hexagonal boron nitrite (hBN) is a laminar van der Waals material with advanced fabrication techniques making it suitable for studying semiconductor physics in two dimensions (2D). In particular, the wide energy gap of hBN may host numerous defects states with internal optical transitions that give rise to color centers. Some of these centers have already shown great potential in quantum technology application 1-7 demonstrated previously in the bulk semiconductors, such as diamond [8][9][10][11] and silicon carbide [12][13][14][15][16] . Color centers in exfoliated hBN, however, may exhibit significant advances over their bulk counterparts owing to the distinct properties of quasi 2D semiconducting flakes.
In recent years numerous color centers in hBN 1-4,17-23 have been reported as new generation high temperature single photon emitters. While the microscopic configuration and electronic structure of these emitters are still not fully understood, it is widely accepted that these color centers can be associated with point defects and point defect complexes. Besides favorable optical properties, point defects often possess spins that can implement isolated quantum bits or qubits. So far previous investigations have focused on the optical properties of the color centers, the spin degree of freedom has been observed only in very recent experiments [5][6][7] .
A key signature of potential point defect qubits is the optically detected magnetic resonance (ODMR) signal which makes optical spin initialization and read out possible through spin dependent decay processes from the excited state to the ground state. Recent experiments have recorded ODMR signal at room temperature for color centers in hBN 6,7 . One of the centers is associated with the negatively charged boron vacancy (VB) based on a previous theoretical study 21 . Although, the charge transition levels [19][20][21]24 and group theory analysis on the electronic structure aided by Kohn-Sham density functional theory (KS DFT) calculations 21 have been reported for this defect but those studies neglected correlation effects and dynamic effects due to phonons and mostly focused on bright states. On the other hand, it has been recently shown from first principles for an exemplary point defect qubit, nitrogen-vacancy center in diamond (see refs [25][26][27] ) that phonon mixing of highly correlated dark states is a key in the optical spin polarization process. In addition, the spin-dependent properties of defect qubits, such as zero-field splitting (ZFS) and hyperfine constants (HF), have not yet been reported in hBN which is a key tool in identification of the point defect qubits and plays a crucial role in the quantum optical control of the qubit.
Here, we carry out a thorough first principles investigation on VB in single sheet hBN. In particular, we compute zero-field splitting, hyperfine splitting, and photo luminescence spectrum, all of which are in good agreement with experiment. Many-body spectrum of VB coupled with phonons is demonstrated allowing us to discuss relevant spin selective decay processes that give rise to optical spin polarization and spin dependent luminescence of the defect. Our first principles results not only reproduce very recent experimental results on VB but also provide a theoretical foundation for developing VB based qubit applications.

RESULTS
Electronic structure and magneto-optical properties VB comprises of a vacant boron site and an additional electron captured from the lattice. According to recent ab initio calculations the negative charge state is stable when the Fermi level is located close to the middle of band gap 20 . The highest point group symmetry of VB is D 3h that may reduce to C 2v under the effect of external strain or due to Jahn-Teller distortion. The inplane dangling bonds and the out of plane p z orbitals of the three neighbor nitrogen atoms give rise to two non-degenerate a and two degenerate e single particle defect states, most of which appear in the band gap. The six defect states are occupied altogether by ten electrons in the negative charge state that come from the three fully occupied p z orbitals and the three half-filled dangling bonds of the first neighbor nitrogen atoms and one extra electron captured from the environment. Depending whether the defect state is formed by the in-plane dangling bonds or by the out of plane p z orbitals, we distinguish prime and double prime defect states, respectively. These prime and double prime states play similar role to the parity selection rule for systems with inversion symmetry in the optical transition.
Our first principles DFT calculation provides single particle Kohn-Sham (KS) electronic structure 21 as depicted in Fig. 1a. The electrons of VB occupy the defect states so that the e′ dangling bond state becomes half occupied and give rise to a spin-1 ground state. The spin density of the 3 A 0 2 ground state is depicted in Fig. 1b, which shows similarities to the spin density of known qubits in 3D bulk materials such as the NV center in diamond and the divacancy in SiC 28 . Fine structure of the electron spin sublevels is dominated by the dipolar spin-spin zero-field interaction for which we obtained D s = +3.471 GHz and D b = +3.467 GHz in single sheet and bulk hBN, respectively. These values are in remarkable agreement with the value of D ≈ 3.5 GHz reported by recent experiments for the VB assigned color center in bulk hBN 6 . We note that the experimental results also show an E = 40 MHz splitting, which suggests the presence of strain in the sample or stray electric fields during the ODMR measurements. Nuclear spins of the host hBN lattice give rise to hypefine splitting of the spin sublevels. In Supplementary Tables 1, 2, we provide the calculated hyperfine coupling constants of nitrogen and boron nuclear spins for the most relevant neighboring sites of VB, see Supplementary  Fig. 1. In the calculations, we obtain A z,s = 47.9 MHz and A z,b = 47.2 MHz splitting for the first neighbor 14 N nuclei in single sheet and bulk hBN, which are in excellent agreement with the experimental observation of 47 MHz 6 . Our ab initio ground state spin coupling parameters thus unambiguously identify the experimental color center reported in ref. 6 as VB. Furthermore, our results demonstrate that the ground state coupling parameters are essentially identical in free standing single sheet and bulk hBN.
To obtain the many-body spectrum and the corresponding excited states, we apply the inherently multireference density matrix renormalization group method (DMRG) 29 generalized for arbitrary long-range interactions together with optimization tools based on quantum information theory 30 using Kohn Sham orbitals to span the active space. For further details, see Methods section. The vertical excitation spectrum as obtained on the ground state configuration of VB is provided in Supplementary Table 3 and depicted in Fig. 2a. In the same figure, we provide leading order expression of the excited states electron hole density. We note that the spectrum shows similarities but also difference to previous reports which did not account for high correlation effects 21 . In the triplet manifold, we find two excitation branches. The lower branch contains nearly degenerate 3   Orange and red arrows indicate perpendicular and parallel polarized optical transitions, while blue and green arrows indicate coupling between singlet and triplet states via perpendicular and parallel spin-orbit coupling, respectively. In all cases, leading order Slater determinant expression of the electron hole density is provided. JT and PJT labels Jahn-Teller(JT) and partial Jahn-Teller (PJT) states. Dashed arrows highlight transitions allowed by JT effects. b Simplified electronic structure for understanding spin dependent non-radiative decay processes. Dashed arrows indicate second order, JT, PJT, and strain allowed transitions. ex1 and scd1 label optical excitation and spin conserving decay, while isc1 and isc2 label inter system crossings relevant for the optical spin polarization, respectively. c, d Many-body electron structure as generalized to strained cases by considering c C 3v symmetric and d C 2v symmetric reconstructions.
the mixture of e″ and a 00 2 defect states, respectively. The higher lying branch contains a single state 3 E′ which account for excitation from the a 0 1 state to the e′ state of VB. In the singlet manifold, we find 1 E′ state with the lowest energy, however, the 1 A 0 1 state of similar determinant appears at very high energy. Inbetween these states we find the singlet counterparts of the 3 A 00 2 and 3 E″ excited states, i.e., 1 A 00 2 and 1 E″ excited states. These results demonstrate that static correlation effects play crucial role in forming the excited state spectrum. We note that orbital and structural relaxation effects as well as dynamical screening effects decrease the many-body vertical excitation energies 26 provided in Supplementary Table 3.
Based on group theory considerations in D 3h symmetry of the ground state, we find several allowed optical transitions, yellow and orange arrows in Fig. 2a, and spin-orbit interactions, blue and green arrows in Fig. 2a, that may allow spin dependent nonradiative transitions. In the following, we restrict our considerations to the lowest three excited states, 3 A 00 2 , 3 E″, and 3 E′, in the triplet manifold and to the three lowest excited states, 1 E′, 1 A 00 2 , and 1 E″, in the singlet manifold. First, we investigate optical decay processes that may account for the photoluminescence spectum reported recently for this center 6 . Note that the lowest energy excited states in the triplet manifold, 3 A 00 2 and 3 E″, are optically not allowed 21 and transition from the 3 A 0 2 ground state is only possible to the higher lying 3 E′ states with perpendicular polarization. These selection rules can be relaxed by strain and electric field that may reduce the symmetry to C 3v and C 2v depending on the character of the external field and thus give rise to additional optical excitation pathways, see Fig. 2c, d. Accordingly, 3 E″ and 3 A 0 2 states transform to 3 E and 3 A 2 in C 3v allowing perpendicular polarized photon emission and absorption, while the 3 E″ state splits as 3 A 2 and 3 B 2 in C 2v allowing parallel polarized photon mediated optical transition between the 3 B 1 ground state and the 3 A 2 excited state.
Similarly to the effect of strain, e′ ⊗ E″ Jahn-Teller (JT) effect may give rise to additional transitions. In Fig. 2a, we marked all the JT unstable states in D 3h symmetry. Allowed transitions due the coupling with e′ phonon modes are marked by dashed arrows in Fig. 2a. The dark 3 E″ state combined with an e′ local vibration mode may give rises to polaronic states that can be connected with the ground state through parallel polarized photon emission and absorption.
Structural relaxation and JT effects in the excited states can be investigated by ab initio DFT methods, for details see Methods section. Indeed, our DFT simulations reveal considerable JT distortions that release 0.193, 0.297, and 0.132 eV energy in the 3 E″, 1 E″, and 1 E′ states, see Table 1. We note that the 1 E′ singlet state, dominated by the e′ 2 determinant, is JT unstable only due to the inter mixture of JT unstable e′a 0 1 Slater determinant. Besides e′ phonons, membrane vibrational modes 31 also couple to the electronic structure of VB in the 3 E″ and 1 E″ states that further reduce the symmetry and the energy, see Table 1. While the JT and membrane distortions of the first neighbor shell of VB are in the same order, ≈0.1 Å, the energy gain due to out of plane distortions is an order of magnitude smaller than the JT energy. This indicates that the potential energy landscape is considerably flatter for out of plane structural relaxations. The 1 E′ state partially inherits this behavior by mixing with 1 E″ state when planar reflection symmetry is broken. The resultant out of plane relaxation is ≈0.01 Å for the first neighbor nitrogen atoms, while the energy gain is smaller than 1 meV.
Due to the fact that 3 A 00 2 and 1 A 00 2 states do not experience JT distortion, whereas 3 E″ and 1 E″ states relax substantially due to JT effect, we anticipate that the latter state are lower in energy after geometry reconstruction. The spectrum of VB as obtained on corresponding relaxed excited state configurations is provided in Table 1. Indeed, 3 E″ and 1 E″ are the first and second excited states in the triplet and singlet manifold, respectively. We find allowed optical transitions with parallel polarization and zero phonon luminescence (ZPL) energies 1.710 and 1.318 eV between these states and the ground and the lowest energy singlet state. The calculated 300 K emission side band of the optical transition in the triplet manifold as approximately determined by ground state phonon modes is depicted in Fig. 3a. The spectrum of Huang-Rhys factor at 3.5 fairly agrees with the experimental spectrum. Disagreement between theory and experiment may be related to lower detection efficiency at high wavelength (>1000 nm) in the experiment. We obtain best agreement between the spectra by assuming 1.62 eV energy (765 nm) for the ZPL emission, which is not resolved in the experiment. This value and the first principles triplet ZPL excitation energy agree very well, thus we assign the experimental optical spectrum to the optical transition from the 3 E″ JT system to the 3 A 0 2 ground state. To study the polaronic solution of the triplet excited state in more details, we first calculated the full adiabatic potential energy surface and we obtained barrier energy between different JT distortions of 152 meV. By solving the e′ ⊗ E″ JT Hamiltonian 32 , we obtained the polaronic spectrum as listed in Supplementary  Table 4. The first polaronic e A 00 1 excited state lies only 0.01 μeV above the polaronic e E 00 ground state that indicates effectively static JT distortion of the 3 E″ state. The e A 00 1 state has an allowed optical transition towards the electronic A 0 2 ground state by emitting a photon with parallel polarization. This finding justifies the application of Huang-Rhys theory for calculating the ZPL energy and the phonon sideband of the photo luminesence (PL) spectrum as explained in ref. 27 . The calculated radiative lifetime is  Fig. 3 Photoluminescence spectrum of VB at 300 K as obtained from ab initio simulations for the 3 E″ → 3 A 0 2 transition. To produce the theoretical curve the calculated Huang-Rhys factor of 3.5 is used whereas the position of ZPL at 765 nm (1.62 eV) was aligned to match to the highest intensities of the experimental spectrum which corresponds to about 0.09 eV redshift compared to the KS DFT result (1.71 eV). This difference is within the usual inaccuracy of the KS DFT method. The experimental spectrum is reported in ref. 6 . about 20 μs in the distorted structure that we obtained in single particle approximation. This finding is in accordance with the second order optical transition. We note that the calculated absorption from the electronic ground state is about three orders of magnitude larger for the optically allowed 3 E′ state. Thus, optical driving of the center can be much more effective with photoexcitation towards the second triplet excited state 3 E′ and then the electron will decay to the lower 3 E″ state. Indeed, the 2.33 eV excitation used in the experiment 6 is sufficient to excite the system to the optically allowed 3 E′ state, see Table 1.
Spin selective non-radiative decay The excited state structure of VB offers several alternative decay pathways from triplet optically excited states through singlet excited states to the ground state. These pathways may include spin non-conserving transitions through inter-system-crossings (ISC), where certain many-electron states become nearly degenerate. The interplay of phonons, that drive the system to a configuration of an ISC, and spin-orbit interaction, that flips the spins, may allow triplet-to-singlet and singlet-to-triplet transitions depending on the symmetry of the involved many-electron states. In Fig. 2 the symmetry allowed spin non-conserving transitions are provided in different point group symmetries.
We describe the most relevant decay processes by a six-state (see Supplementary Fig. 2) rate equation model in Supplementary Note 1. The spin selectivity of the decay is predominated by the strong coupling between the triplet and singlet E″ excited states that we find nearly degenerate, have akin atomic configurations, and couple with parallel spin-orbit interaction of approximately 170 GHz strength as obtained from KS DFT calculation. This corresponds to several orders of magnitude larger rate than the calculated radiative rate between the triplets. Due to this strong coupling and the first order allowed optical transition between the lowest energy singlets, we infer that the decay through the redshifted singlet optical emission, with ZPL at ≈1.32 eV (939 nm), dominates over the triplet optical decay for the m S = 0 spin sublevel. The sign of the ODMR contrast predominantly determined by the slow second order decay processes, such as the direct optical and non-radiative 3 E″ → 3 A 0 2 transition with rate r scd1 and the 1 E′ → 3 A 0 2 (m S=0 ) transition with rate r isc1 , see Supplementary Note 1. The negative contrast reported in the experiment 6 is in agreement with the expectation that r isc1 < r scd1 .
Optical pumping induced ground state spin polarization is governed by the 1 E′ → 3 A 0 2 (m S = ±1) transition rate r isc2 and r isc1 . We note that these rates correspond to such decay processes that are allowed only in second order due to out of plane relaxation and partial JT distortion, respectively. This characteristic infers that the result of optical pumping can depend very much on external conditions, such as strain and electric field, that may influence rates r isc1 and r isc2 . We can distinguish two different scenarios. When r isc2 > cr isc1 , where c is a prefactor between 0 and 1 determined by rate r scd1 and rate r isc4 of 3 E″(m S = ±1) → 1 E′, the spin is polarized in the m S = ±1 sublevels, for which we expect luminescence in the triplet manifold with ZPL of ≈1.71 eV energy. For cr isc1 > r isc2 , the spin is polarized in the m S = 0 sublevel of the ground state for which we expect red-shifted luminescence through the singlet states. As the reported room temperature luminescence fits to the triplet transition, the former scenario is more probable and the spin should be polarized in the m S = ±1 sublevels of the ground state accordingly. This finding is supported by the reported out of plane relaxation of the first neighbor atoms at equilibrium and the soft potential for membrane vibrational modes and out of plane distortions.
The above described spin selective decay mechanisms demonstrate distinct behavior of VB from the known solid state qubits. The features provided by VB can be utilized, for example, to boost the ODMR contrast and to control the optical spin polarization mechanism. For achieving enhanced ODMR contrast one may use conventional optical filters to separate triplet emission associated with m S = ±1 states and the red shifted singlet emission associated with m S = 0 state. This way the ODMR contrast may be significantly increased. Furthermore, additional excitation in the singlet manifold can enhance photon count rate for the m S = 0 sublevel to speed up the read out process. As the ground state spin polarization is determined by second order transitions induced by structural distortions, application of strain may allow control over the spin polarization and thus the photoluminescence of the center. For instance, strain that breaks planar reflection symmetry may enhance m S = ±1 polarization, while strain of E′ symmetry may result in m S = 0 polarization and redshifted luminescence.

DISCUSSION
Our ab initio results explain the known magneto-optical properties of VB in hBN. These achievements validate the electronic structure and spin polarization mechanism described in the results section. Having this in hand, one may apply more advanced techniques to harvest functionalities allowed by the electronic structure of VB. So far, only one optical transition is observed in the experiment, however, we predict multitude of them. Application of two or more excitation lasers may be utilized to control and improve optical initialization and read out properties of the defect. Furthermore, we predict extraordinaire capability for VB to control its magneto-optical properties by engineering strain.
We note that our calculations were performed dominantly in a single sheet model, while the experiment were carried out in bulk samples. The good agreement between theory and experiment indicates that VB should function in free standing single sheet hBN as good as in bulk, which is an important requirement for lowdimension applications.
Finally, in our study, we utilized a so far less recognized combination of methods: density matrix renormalization group method on active space spanned by localized defect orbitals defined by hybrid density functional calculations. Successful application of the method on the involved electronic structure of VB defect demonstrate that this method may be a key tool in investigating functional color centers in hBN and in other wide band gap semiconductors.

Density functional theory calculations
We apply two methods, DFT and density matrix renormalization group (DMRG), to describe the electronic structure of VB on different levels of theory. For the DFT calculations plane wave basis set of 450 eV and PAW 33 atomic potentials are used as implemented in VASP 34 as well as plane wave basis set of 750 eV and norm-conserving pseudo potentials are used as implemented in Quantum Espresso (QE) 35 . HSE06 hybrid functional 36 with 0.32 exact exchange fraction 20 is used for hyperfine calculations 37 , excited state calculation in the framework of constrained occupation DFT 38 , and structural relaxation. We use 162 atom super cell of single sheet hBN embedding a single boron vacancy. In perpendicular direction, we use 30 Å supercell size. Our bulk hBN model consists of 972 atoms (9 × 9 × 3 primitive cells) and includes a single boron vacancy. HSE06 functional is used to calculate hyperfine 37 and spin-spin zero-field-splitting parameters. To eliminate spin contamination in the latter case, we apply the correction scheme proposed in ref. 39 . To determine the phonon spectrum of the ground state VB configuration we use PBE functional. Phonon sideband and polaron spectrum are obtained by the machinery described in ref. 25 .
In the hybrid DFT constrained occupation calculations proper excited state configurations are identified by comparing spin density and partial charge density provided by DFT calculations with the true many-body spin and charge densities obtained from DMRG calculations. We note that constrained occupation hybrid DFT calculations often provide non-physical solutions. In particular, this is the case for higher lying excited states. Furthermore, the triplet and singlet E″ states are multi-determinant as V. Ivády et al.
provided by DMRG calculation. We approximate this state by setting the excited hole as a mixture of a 00 2 and e″ states. Note that the spin density of this mixture is a good approximation to the spin density of the many-body excited state.

DMRG calculations
The studied sample, consisting of more than a hundred atoms, is described by several hundreds of Kohn-Sham orbitals making prohibitively expensive to be tackled directly by the DMRG approach. Therefore, a careful selection of active space to be treated on post-DFT level is important describing the remaining single particle states, i.e., which are not included in the active space, in frozen core approximation. Based on previous studies on defect systems 26 , it is expected that the low lying energy spectrum is predominantly determined by the highly localized defect orbitals. Therefore, the discussed results are based on DMRG calculations performed on active space of 27 spin restricted Kohn-Sham DFT orbitals which are selected according to the localization of the orbitals around the defect, i.e., occupied orbitals kept in the active space whose degree of localization on the first neighbor three nitrogen atoms reaches at least 0.10. In order to preserve the point group symmetry of the orbitals, we studied canonical orbital set without further localization.
To verify the applicability of the active space selection based on localization, we also performed DMRG calculations defining the active orbitals based on an energy window around the localized defect levels. The structure of the many-body excited state spectrum yielded on such active space is essentially identical to the presented one observing upward shifts in energy in the range of 0.01-0.15 eV. Obtaining slightly lower energies with the applied variational computational procedure, e.g., by −0.06 eV for the ground state, indicates that the orbital selection based on localization is preferable to the one based on energy. We apply the DMRG method on spin-restricted Kohn-Sham orbitals computed by QE using PBE functional. In these SCF calculations we use fixed geometries obtained on HSE06 level of theory by VASP. One-particle and two-particle integrals of the Hamiltonian matrix of the QE Kohn-Sham orbitals are generated by our in-house code. DMRG calculations are performed using the Budapest DMRG code 40 . Considering the large number of excited states to be computed, in order to enhance convergence, distinct DMRG calculations are performed fixing the total spin of the target states to 0 and 1. The numerical accuracy of the calculations was controlled by using elements of quantum information theory and by the dynamic block state selection approach 41 keeping up to thousands of block states for the priory set quantum information loss threshold value χ = 10 −5 . The density matrix of the system block in the DMRG truncation procedure is formed of the equally weighted linear combination of all target states.

DATA AVAILABILITY
The data that support the findings of this study are available from the authors on reasonable request, see author contributions for specific data sets.

CODE AVAILABILITY
The DFT supercell codes that were used in this study are standard codes: VASP (version 5.4) is a commercial code (see www.vasp.at) whereas Quantum Espresso (version 6.1) is an open source code (see www.quantum-espresso.org/). The electronphonon code, the QC-DMRG-Budapest code and its interface to the Quantum Espresso code that were used in this study are available upon reasonable request to the corresponding author.