Skip to content
BY 4.0 license Open Access Published by De Gruyter February 1, 2023

Recent advances in the ab initio theory of solid-state defect qubits

  • Ádám Gali ORCID logo EMAIL logo
From the journal Nanophotonics

Abstract

Solid-state defects acting as single photon sources and quantum bits are leading contenders in quantum technologies. Despite great efforts, not all the properties and behaviours of the presently known solid-state defect quantum bits are understood. Furthermore, various quantum technologies require novel solutions, thus new solid-state defect quantum bits should be explored to this end. These issues call to develop ab initio methods which accurately yield the key parameters of solid-state defect quantum bits and vastly accelerate the identification of novel ones for a target quantum technology application. In this review, we describe recent developments in the field including the calculation of excited states with quantum mechanical forces, treatment of spatially extended wavefunctions in supercell models, methods for temperature-dependent Herzberg–Teller fluorescence spectrum and photo-ionisation thresholds, accurate calculation of magneto-optical parameters of defects consisting of heavy atoms, as well as spin-phonon interaction responsible for temperature dependence of the longitudonal spin relaxation T1 time and magneto-optical parameters, and finally the calculation of spin dephasing and spin-echo times. We highlight breakthroughs including the description of effective-mass like excited states of deep defects and understanding the leading microscopic effect in the spin-relaxation of isolated nitrogen-vacancy centre in diamond.

1 Introduction

Quantum information is physical [1]. Solid-state defect spins are a conceivable platform to realize the elementary unit of quantum information, i.e., quantum bits or qubits [2]. Two prototypical representatives are the phosphorus donor (P-donor) in silicon and the nitrogen-vacancy (NV) centre in diamond. From electronic structure point of view, these two defects reside the opposite sides of the spectrum: the P-donor can be described by the so-called effective mass state with hydrogen-like Rydberg series of excitation energies split from the conduction band of silicon, which are weakly localized wavefunctions, whereas the carbon dangling bonds in NV centre create strongly localized orbitals with deep levels in the fundamental band gap of diamond. Kane proposed to apply P-donor spins as a qubit [3]. However, the read-out of the single spin in a controlled fashion had been great challenge for relatively long time that could be realized in a single-electron transistor device operating at hundreds of millikelvin temperatures [4]. The read-out of the single spin of diamond NV centre has been realised optically, i.e., optically detected magnetic resonance (ODMR), which was the first optically read single defect spin in a solid [5]. In this case, the readout and initialization of the electron spin of diamond NV centre could be readily carried out by optical means at room temperature. Recently, the single electron spin electrical read-out via photo-ionised electrons and holes has been realised for diamond NV centre which is a hybrid scheme: photo-excitation is required for creating spin-dependent photocurrent from a single NV defect and optical initialisation of the spin, and then the spin-dependent photocurrent is observed to read the electron spin state [6, 7]. The coherent manipulation of these electron spins were realised electron spin resonance techniques [8, 9]. The coherent control and readout of single defect spins define the underlying defects as quantum defects, and the quantum defect with its host material can be called as a quantum-coherent material.

Since the discovery and realisation of these quantum defects an intense research has begun to seek alternative solid-state defect spins both in the experimental and theoretical fronts, which might have favourable properties for certain quantum technology applications [10, 11]. Recently, data mining techniques with machine learning algorithms have been spread at the theoretical fronts. The data mining can be approached either towards the host materials [1214] or to extend this with creating defect structures and calculating their key qubit parameters [15, 16]. In these approaches, there are assumptions for selection of materials and defects that might be too restrictive and might lead to overlooking important candidates. For instance, it was assumed that the host materials should have wide band gap with low density of nuclear spins, at least, for the defect qubits alike diamond NV centre [10, 12]. However, certain quantum technology applications such as quantum communication does not require room temperature operation, and small-band-gap silicon has become a promising platform to realize spin-to-photon interface with quantum memory [11, 17], [18], [19], [20], [21], [22], [23]. One of the most promising platforms to host qubits in two-dimensional (2D) materials is hexagonal boron nitride (hBN) which has 100% nuclear spin abundance but the coherence times of defect spins can be well extended with using good control of these nuclear spins [2427]. Certainly, the selection criteria can be changed based on these recent findings. Nevertheless, quantum-coherent materials can only be interpreted together with their defect qubits, thus selection of host materials should be followed by finding defects for which the electron spin can be initialised and read-out with sufficiently long coherence times. Calculation of the coherence times for any hypothetical defects, e.g., Ref. [13], could lead to misleading results in certain cases where the spin density distribution and so the strongest hyperfine interaction between the electron spin and nearest nuclear spins characteristic for the actual defect would strongly affect its electron spin’s coherence time [24]. The automatically generated defects are often selected from the thermodynamically most stable ones [15, 16]. It has been found that this selection may omit very important complex defects realizing qubits, e.g., G-centre in silicon [28], which is one of most promising qubit candidates in silicon. Furthermore, it has been shown for diamond NV-centre [29] that the strongly coupled electron-phonon states are inevitable for understanding the optical spin-polarisation of its electron spin. These polaronic states are also usually ignored in these databases. These examples clearly call for improving the ab initio magneto-optical spectroscopy methods, in order to increase the credibility and prediction power of these databases.

Ab initio methods have significantly contributed to understanding and control of diamond NV qubits and exploration of alternative quantum defects which was summarized in a recent review paper [30]. In that review paper, the diamond NV qubit is thoroughly described including the electronic structure and polaronic solutions and spectra within Jahn–Teller theories, and the desiderata of the target of computations and the developed computational methods are summarized in detail which are not perpetuated here. We assume that the readers are aware of the basic description of diamond NV-like qubits and the previously implemented methods to compute their magneto-optical properties which are the basis of further developments.

For instance, new findings have been reported for the very basic property of this qubit such as the spin-lattice relaxation (T1) time [31] which apparently demonstrates that even the most studied diamond NV-centre has been not fully understood to date. It is an immediate quest to further advance ab initio methods, in order to accurately calculate the excited states together with quantum mechanical forces, the electron-phonon coupling and the basic magnetic parameters such as the zero-field-splitting between the electron spin’s levels (ZFS), or understanding the temperature dependence of these parameters and the spin-lattice relaxation (T1) time, and the coherence time (T2) also as a function of spin bath around the target defect qubit. This paper provides a comprehensive review about recent years ab initio developments on solid-state defect spin qubits along these directions.

We illustrate the advances of computational methods on various defects in solids: we shall discuss (i) deep defects such as NV centre, silicon-vacancy centre and nickel-vacancy centre in diamond, divacancy and vanadium centres in silicon carbide, and boron-vacancy in two-dimensional hexagonal boron nitride; (ii) the shallow excited states of deep defects such as the neutral silicon-vacancy centre in diamond and various interstitial defects in silicon; (iii) and shallow donors in silicon.

2 Computational methods

The ab initio investigation of solid-state defect qubits alike diamond NV-centre or silicon P-donor requires full and accurate description of the host material and the embedded isolated defect. The photostability of the quantum defects depends on the ionisation threshold energies, therefore it is critical that the crystalline bands and the resonant or localized defect states are computed at equal footing. Green-function methods are principally ideal to represent the topology of the problem, i.e., the embedded defect in a perfect solid. However, practical implementations of Green-function methods suffer from the consistent calculation of quantum mechanical forces which is required to calculate the ionic coordinates of the given point defect. As the geometry of the defect is highly decisive in their magneto-optical properties, the supercell method is most often employed to model the quantum defects which methodology readily offers the calculation of quantum mechanical forces based on the Hellman-Feynman theorem [30]. In this review article, we focus on those method developments and implementations which work within this formalism because the supercell formalism guarantees the simultaneously accurate calculation of the ionisation threshold and intra-defect optical transition energies. The ground state properties are typically calculated by means of Kohn-Sham density functional theory (KS DFT) [32], which can be a starting point for the calculation of excited states. A natural choice for the basis set for supercell formalism is the plane-wave basis which is combined by pseudopotential or projector augmentation wave (PAW) methods (see Ref. [33] and references therein). The computation methods of the ground state thermodynamic properties of point defects in solids within this formalism was already described in detail in the literature (e.g., Ref. [34]) that can be applied to the specific quantum defects [35]. In short, the quantum defect’s local properties and associated parameters are computed within plane wave supercell KS DFT methods, and implementation based on these methods or relied on the parameters obtained by these methods will be discussed in the review paper.

3 Method developments and results

In this section, we collect recent developments on ab initio calculation of properties of quantum defects. We start with the treatment of the excited states and geometries, and then we continue with the discussion of the ionisation threshold energies with the photo-ionisation spectrum, the supercell treatment of shallow defect states, the zero-field-splitting parameters, hyperfine parameters and gyromagnetic tensor, and the spin-lattice relaxation T1 times as well as the spin coherence T2 times. Although, it is unconventional to start the discussion with the excited states prior the description of finite size effects of supercells we decided to do so because we shall consider the finite size effects of shallow excited states in the section of supercell modelling of defects which requires description about the calculation of excited states. We shall also show that accurate description of the excited states is needed to calculate some ground state related magnetic parameters such as the zero-field splitting parameter.

3.1 Excited states

3.1.1 Following the topology of the problem: DFT + CI multiscale methods

The accurate calculation of the excited states is still under extensive research. The popular many-body perturbation method, GW + BSE (see Ref. [36] and references therein), fails for the highly correlated singlet states of the diamond NV centre [37, 38]. The highly correlated states may be recognized of which many-body electronic wavefunction (Ψ) can be genuinely described by two or more Slater-determinants with significant weights. This may be quantified with the function of one-particle density matrix ρ̃. This can be defined as ρ̃ij=Ψ|aiaj|Ψ with the creation and annihilation operators of single electrons i and j, respectively. The degree of correlated electronic state (Λ) is then

(1)Λ=Trρ̃ρ̃2=Trρ̃Trρ̃2,

as if and only if the Ψ can be described as a single Fock-state when ρ̃=ρ̃2. Those cases are in particular pathological in which double-excitation Slater-determinants appear with relatively strong contribution (i.e., aiajakal|ΨGS with the ground state many-body wavefunction ΨGS) which is the case for the so-called 1A1 state of the diamond NV centre (see Figure 1) which plays an important role in the optical spin-polarisation and read-out processes (see Ref. [30] and references therein).

Figure 1: 
Defect levels diagramme. (a) Single-particle scheme of the electronic structure of diamond NV centre. The spin-polarisation between spin-up and spin-down electrons in the Kohn–Sham density functional theory results in different levels for spin-up and spin-down orbitals which are not depicted here but rather their average. VB and CB label valence band and conduction band, respectively. The fundamental band gap of diamond is 5.4 eV. (b) Many-body level structure of diamond NV centre. The blue and red shaded areas show the ionisation bands excited from the 3A2 and 3E states, respectively. In the bottom row the non-physical picture is shown which often appears in many scientific papers.
Figure 1:

Defect levels diagramme. (a) Single-particle scheme of the electronic structure of diamond NV centre. The spin-polarisation between spin-up and spin-down electrons in the Kohn–Sham density functional theory results in different levels for spin-up and spin-down orbitals which are not depicted here but rather their average. VB and CB label valence band and conduction band, respectively. The fundamental band gap of diamond is 5.4 eV. (b) Many-body level structure of diamond NV centre. The blue and red shaded areas show the ionisation bands excited from the 3A2 and 3E states, respectively. In the bottom row the non-physical picture is shown which often appears in many scientific papers.

In KS DFT theory, these many-body states cannot be captured by the known exchange-correlation potentials (Vxc); however, it provides a very good description for the host materials and simple defect states, i.e., for which Λ = 0. In general, the bands of the host crystal can be well described by KS DFT method whereas the strong Coulomb interaction between localized defect orbitals confined to a small place, e.g., around a vacancy, represents a problem for KS DFT method. The topology of the problem implies a method in which the Coulomb interaction between the localized defect orbitals is directly calculated, called configurational interaction (CI), which are in contact with the bath of itinerant electrons that can be treated by KS DFT. This can be considered as a multiscale problem where the interaction between electrons are calculated with different precisions in a single system. As we will see below the challenges are to define the set of orbitals, often called “active space”, for which the precise method should be applied (here, CI) and to find an interface between the levels of different approaches, here DFT and CI, which can produce self-consistent results. Interfacing different levels of theory in a single system is a common problem for all the multiscale methods.

By following the above mentioned topology of the problem, the many-body Hamiltonian may be described as

(2)Ĥ=ijAtijeffaiaj+12ijklAvijkleffaiajalak,

where veff is the partially screened Coulomb-interaction W0R and the ijkl KS wavefunctions are within the active space A. Partial screening is computed from all the KS wavefunctions except for the set of KS wavefunctions within A. The definition of veff includes contributions to the Hartree and exchange correlation energies that are also included in the KS DFT calculations for the whole solid. Exactly this is the point where an interface between DFT and CI methods should be developed which can be considered as a typical problem for all the related electronic structure embedding schemes called double counting (dc) error. Therefore, the KS DFT Hamiltonian ĤijKS, occurring in the first term in Eq. (2), requires double counting correction tijdc, i.e.,

(3)tijeff=ĤijKStijdc.

The methodology was first implemented by Bockstedte and co-workers within VASP code [39, 40] which was called DFT + CI-cRPA method [41]. Here, CI refers to the configurational interaction between the electrons in the active space whereas constrained random phase approximation was applied to compute W0R [41]. In that implementation, Heyd–Scuzeria–Ernzerhof (HSE06) [42] functional was applied which includes a quarter of Fock-exchange (α = 25%) in the KS potential for the KS orbitals. To calculate tijdc, a heuristic approach was applied, i.e., the quarter of Fock-exchange was used for the electrons in the active space. This approach resulted in a very good agreement between the computed spectra and the spectra derived from experiments for diamond NV centre and silicon carbide divacancy centres [41]. Later this treatment was also implemented to the Quantum Espresso software package [43] and the CI part of the Hamiltonian was also interfaced with a code running on quantum computers [44]. The quality of present quantum computers precludes to achieve accurate results [44], thus the results achieved by classical computers will be summarized below. The DFT + CI multiscale method was then rephrased to quantum embedding or quantum defect embedding method but it is essentially the same with the originally invented CI-cRPA method with the caveat that the latter implementation used a full Hartree–Fock dc correction (HFDC), i.e., the exchange term was not scaled by the fraction of the applied hybrid DFT functional [41, 44]. Ma and co-workers also went beyond RPA (b-RPA) by taking the exchange-correlation effects into account in the calculation of W0R, and they applied it to diamond NV centre [44]. It was found that b-RPA screening substantially modifies only the 1A1 level referenced to the ground state energy with respect to the results within RPA screening approach (see Table 1). The origin of this effect is not explained and well understood. Pfäffle and co-workers implemented a semiempirical DFT + CI method [45]: they basically applied HFDC to the DFT or GW quasi-particle (GWQP) levels in Eq. (3). The screening was taken from the bulk diamond with an analytical formula developed for semiconductors which provided the accuracy of the RPA method [46]. Resta found that the space-dependent dielectric screening only varies within the range of nearest-neighbour distance of ions in semiconductors, if the origin is chosen to an ion in the crystal, with a well-defined function and then it approaches the dielectric constant of the material, ɛ(0) [46]. Taking the host material’s dielectric screening function for a defect might be too inaccurate, in particular, for defects which contain vacancy as no ion would “modulate” the dielectric constant of the host material in the region of vacancy. Therefore, Pfäffle and co-workers introduced a semi-empirical formula to modify the dielectric screening function near the vacancy [45] of diamond NV defect which they called “masked” solution. It is unclear how the parameter in that formula was picked up. With the “masked” solution for the diamond NV centre they obtained triplet optical transition close to the DFT + CI-cRPA solution [41] whereas the energy gap between the singlets closely follows that of b-RPA screening with HFDC [44]. We note that this method was further applied to the neutral and positively charged NV defect in diamond. The optical spectrum is known for the neutral NV defect in diamond [47, 48] for which the ZPL energy results in 2.156 eV between 2E ground state and 2A2 excited state. Although, the order of the many-body levels is consistent with the experimental data by the semi-empirical DFT + CI method but the calculated energy gap for 2E2A2 transition is at 1.65 eV which is disappointingly low. It is yet unclear what is the origin behind the discrepancy.

Table 1:

The calculated electronic structure of diamond NV centre with using different types of DFT + CI embedding methods. The ground state is 3A2 and all the levels are referenced to it with fixed ionic coordinates of the electronic ground state. The type of DFT functional is given: the DDH functional provides α = 17.8% Fock-exchange for diamond [44]. In HSE06 functional [42] α = 25% was applied. The applied approximation for calculating W0R is either cRPA or goes beyond (b-RPA). Finally, the dc correction is either Hartree–Fock type (HFDC), or scaled with α in the DFT functional (αFDC). The exact dc (EDC) correction is applied together with quasiparticle (QP) levels within self-consistent GW method (see text). The experimental data are only estimates from PL spectrum and rate modelling of the observed decay rates (see Ref. [30] and references therein).

NV HSE06 [41] HSE06 [49] DDH [44] DDH [44] GWQP [50] GWQP [50] Exp.
cRPA cRPA cRPA b-RPA cRPA cRPA
αFDC HFDC HFDC HFDC HFDC EDC
1 E 0.49 0.49 0.476 0.561 0.375 0.463 ∼0.40
1 A 1 1.41 1.39 1.376 1.759 1.150 1.270 ∼1.55
3 E 2.02 2.06 1.921 2.001 1.324 2.152 ∼2.15

In any embedding methods it is a crucial question how to select the active space orbitals. The in-gap defect states are per se localized, therefore, they are naturally involved in the active space. However, defect states resonant with the valence band or conduction band may exist which are coupled to the Bloch-states of the host material, thus they do not appear as sharp resonances but are rather broadened. In practice, the coupling strength and the width of broadening depend on the size of supercell and k-point sampling of the Brillouin-zone for the given supercell. For instance, it was found that whereas the lowest energy spectrum of diamond NV centre may converge relatively well for the minimal active space taken from the in-gap defect orbitals but this treatment leads to a serious inaccuracy for the isovalent silicon carbide divacancy defect for which explicit involvement of the resonant states in the valence band is necessary to achieve converged results [41]. This result clearly demonstrated that simplified Hubbard-U model with using only the in-gap defect orbitals (e.g., Ref. [38]) is generally insufficient for accurate description of the many-body electronic structure of quantum defects. In practice, an energy region of about 3 eV around the Fermi-level was used to pick-up the states for the active space with 512 and 576-atom supercell models in the original study [41]. In later studies, the choice of the active space was further investigated [49, 50]. Muechler and co-workers also implemented the DFT + CI-cRPA method to the VASP code but they first constructed Wannier wavefunctions [51] with preserving the position of the in-gap defect levels which could maximally localize the defect wavefunctions including those that are resonant with the host bands, and then the Wannier orbitals were applied in Eq. (2). With this treatment the results upon the number of orbitals in the active space can rapidly converge [49]. Another possibility is to pick-up the states based on localization of the KS orbitals that can be simply quantified as LVψiKS=VΩψiKS2dV, where V is a chosen volume including the defect, smaller than the supercell volume Ω. It was found that L V = 5% is needed to obtain convergent results for the notorious 3E state which corresponds to about 40 (30) states in the active space for 216-atom (512-atom) supercell model [50]. The 3E state is notorious in terms of L V because a1 hole orbital is involved in this state of which level lies close to the valence band which leads to exchange-correlation coupling to resonant states in the valence band. This is not the case for the 1E and 1A1 levels with very little contributions from the a1 hole orbital, so they converge fast in this regard.

Another critical question in the embedding methods is the treatment of double counting correction which has been briefly considered above. In Hartree–Fock methods, the double counting correction can be readily derived and it is applied routinely in quantum chemistry. However, DFT is applied to calculate the electronic structure of the ground state in the supercell modeling of defects, and there is no theoretical rationale to apply HFDC in this case. By applying hybrid DFT functionals, Fock-exchange is employed to KS orbitals. By assuming a tiny contribution of the semilocal exchange-correlation in dc correction, one may apply αFDC meaning that the corresponding fraction of the Fock-exchange for the active space orbitals is employed. This heuristic treatment may follow the idea of the DFT + U treatment in which U is an orbital-dependent on-site correction, and that is scaled by α when hybrid DFT is applied [52]. However, hybrid DFT is still a functional of the electron density and not the many-body wavefunction, thus this heuristic treatment strictly cannot be justified.

Recently, Sheng and co-workers derived the dc correction for KS DFT with using the Green-function approach [50] which might be motivated by the success of interfacing KS DFT and dynamical mean field theory (DMFT) with the same approach [53]. They assumed that the non-diagonal terms of the self-energy coupling the active space and the environment is negligible, and static approximation is used for the G0 and W0. This results in

(4)tdc=Vxc+W0RρAiG0RW0,

where ρ A is the density matrix of the electrons within the active space, and G0 and W0 are computed at the quasiparticle energies (ϵQP) of these electrons [50]. The exact dc correction is called here EDC in which the quasiparticle (QP) levels are solved self-consistently with calculating ∑xc self-energy within GW method as ϵiQP=ϵiKS+ψiKS|xcϵiQPVxc|ψiKS. This makes the ϵQP results almost independent from the applied DFT functional within 0.1 eV, it was either the semilocal Perdew–Burke–Ernzerhof (PBE) [54] or the dielectric-dependent hybrid (DDH) [55] functional. The results are summarized for diamond NV centre in Table 1 where they applied both HFDC and EDC corrections. They reported extremely low energy for 3E state with HFDC which was not explained in detail [50]. The EDC results are quite comparable to the genuine CI-cRPA result with the heuristic αFDC correction, nevertheless, the good agreement might be specific to diamond NV centre.

Sheng and co-workers applied the EDC method also to the neutral silicon-vacancy (SiV) defect and their sister group-IV vacancy defects in diamond [50]. The ground state of the defect is 3A2g which is a similar wavefunction to that of diamond NV centre. However, the excited state of the neutral SiV defect is highly complex with leaving a hole in each double degenerate orbital which results in three triplet states (3A2u, 3E u , 3A1u). These triplet states are highly correlated and they are also coupled by phonons via product Jahn–Teller interaction [56, 57]. According to that study, the product Jahn–Teller interaction leads to a strong ionic relaxation upon excitation (0.3 eV), so the vertical excitation energy of the optically allowed 3E u state should lie at around 1.6 eV above the ground state level. The EDC calculation yielded disappointedly high energy for 3E u at 2.161 eV when L V = 5% was applied [50]. We note that the calculated excitation energies for the triplet excited states do not converge monotonously with decreasing the value of L V and their values are close to the desired energies with setting L V = 20%. It is unclear yet what is the origin of this behaviour.

Despite the remaining issues with the DFT + CI multiscale methods, this approach is very promising to calculate the highly correlated defect states within supercell formalism. The very important next step is to derive the quantum mechanical forces which is then can be used to reoptimise the geometry of the ground state and the excited states. This is required to calculate the characteristic zero-phonon-line energies and Debye–Waller factors of the defects which are key parameters of quantum emitters. The big challenge here is that the original derivation of the EDC method assumed that the off-diagonal term of the self-energy coupling the active space and the environment can be neglected but once the forces are considered one has to take into account the change of the character of the underlying KS wavefunction self-consistently with the self-energy through the coupling between the active space and the environment. In general, it is a known problem in quantum chemistry CI methods that the final results may depend on the initial single particle wavefunctions if restricted active space is applied and the underlying single particle wavefunctions are fixed at the electronic ground state manifold. The DFT + CI multiscale methods face to the same problem, principally.

3.1.2 Density matrix renormalization group methods: an alternative multiscale method

Density matrix renormalization group (DMRG) was originally developed to describe one-dimensional quantum models in solid-state physics with local interactions. The underlying mathematical framework, however, is not restricted to models studied in condensed matter physics or applications to molecular clusters but among many others, it can be also used to study nuclear shell models, particles in confined potential, or problems in the relativistic domain. The success of these developments relies on the efficient factorization of interactions and the optimization of the DMRG network topologies based on concepts of quantum information theory leading to tremendous reduction in computational costs (see Ref. [58] and references therein). In particular, the factorization for the many-body wavefunction Ψ with L spatial orbitals reads as

(5)|Ψ=(n)C(n)i=1Lainiaini|0

with

(6)C(n)=i=1LAi(nini),

where now the spin state (σ) is explicitly written by arrows and n = (n1↑n1↓n2↑n2↓ … nLnL) where n ∈ {0, 1}. The components of the state specific C tensor increase exponentially with system size L scaling as 22L which becomes untractable for few hundreds of electrons. However, the dimension of the matrix product states A i can be optimized in DMRG approach, AiDMRG, truncated to a fixed manageable bond dimension, M, that is, dimAiDMRG[M,M], i.e., M × M tensor. Increasing M, the precision of the approximation is well-controlled approaching variationally the exact solution. In the DMRG protocol, the matrix product state matrices are locally optimized and truncated by minimizing the discarded entanglement between the left and right neighbouring blocks of the matrix product state chain, obtained from the reduced density matrix of the block. The algorithm iterates through the matrix product state chain in a sequential order back and forth until reaching convergence.

In short, application of DMRG method on electron-nuclei systems can be considered as a special wavefunction method that can be used to accurately calculate the static correlation between electrons. Barcza and co-workers extended this method to interface with DFT calculations of quantum defects [58]. In this post-DFT method, the Coulomb integrals of the KS orbitals are directly taken from the supercell DFT calculation in the electronic ground state manifold which is post-processed by DMRG algorithms. Despite the advantage of DMRG method, hundreds of KS orbitals cannot be directly treated owing to the computational costs. Therefore, an optimal selection of orbitals with tractable size is needed which is responsible for the strong static correlations. This was carried by the complete active space (CAS) self-consistent field method which is well-known in quantum chemistry (see Ref. [58] and references therein). The CAS method classifies the set of orbitals to three categories; that is, the so-called core and virtual orbitals are frozen to the mean field level and filled with two and zero electrons, respectively. The third class comprises of the so-called active space orbitals which are populated with the rest of electrons minimizing the energy. Accordingly, the virtual orbitals does not play any role in the corresponding CAS Hamiltonian, whereas the core electrons affect the electrons of the active space through the Coulomb interactions, that is, the Hamiltonian of the active space reads

(7)ĤCAS=Enuc+Ecore+ijtijCASaiaj+12ijklVijklaiajakal.

The one-electron integrals of the CAS space, tijCAS, describe not only the kinetic energy of the electrons in the active space and their attraction to nuclei but also their interaction with the core electrons. Describing the electrons in the active space with the DMRG method, which treats exactly the electron exchange, the one-electron interactions are written as

(8)tijCAS=tij+12c(2ViccjVicjc)

to treat the Coulombic effects of the frozen electrons on the active space orbitals. In other words, this is a multiscale method where the interface between the active space and the environment can be well managed within Hartree–Fock level. Here, the summation runs only on the indices of the core orbitals (c). Finally, the additional energy contribution of the core electrons is summed up in term Ecore, that is

(9)Ecore=2ctcc+cc2VccccVcccc.

In practice, the active space is restricted to the most important orbitals featuring strong correlation as was discussed for DFT + CI methods. Even though the method has limitations to provide correct description of dynamic correlations using the relatively small active space, it captures static correlations with high accuracy providing valuable insights into the low-lying energy spectrum and the essential structure and symmetry properties of the corresponding electronic eigenstates. Note also that, contrary to DFT + CI methods, the CAS Hamiltonian in Eq. (7) does not include the KS energies explicitly but only the KS orbitals by construction. Also, the absolute energies of the states computed from the CAS Hamiltonian are not trivially comparable with counterparts obtained on the DFT level of theory due to the different description of the exchange and correlation effects.

The methodology was first applied to the negatively charged boron-vacancy defect in hexagonal boron nitride (hBN) [59]. Previous ab initio calculations predicted this defect in hBN as a qubit [60] of which ODMR signal was confirmed in experiment [61] and has become a leading quantum defect in two-dimensional (2D) materials. This defect consists of three nitrogen dangling bonds which introduce four closeby levels in the fundamental band gap. Since the vacancy lobes are closely spaced with strong Coulomb interaction they form highly correlated states. The order of the many-body levels and the nature of the states were determined by DMRG method [62]. In the calculations about 50 KS orbitals were used in the active space selected by localization criterion. It was shown that the flake model with about 80 B and N atoms with hydrogen termination of the edge of the flake provides the same electronic structure as the periodic supercell model at DMRG level [58]. This makes it realistic to calculate the electronic structure of defects within the flake model by means of traditional quantum chemistry wavefunction methods [63, 64]. We note that the validity of the flake (molecular cluster) model could be specific to (planar) defects in hBN and may be not applied in general for defects in other 2D materials.

DMRG method was also applied to magnesium-vacancy (MgV) centre in diamond [65]. MgV centre was created by Mg implantation to diamond which has a unique photostability [66] which makes it an interesting quantum emitter. The PL signal was associated with the negatively charged MgV defect which has a similar electronic structure to that of SiV centre [65] with a caveat that the resonant a2u level lies very close to the top of the valence band and the two degenerate e u and e g levels lie in the fundamental band gap. These states are localized on the vacancy lobes which create highly correlated many-body states. In particular, the neutral MgV has overwhelmly complicated electronic structure which needs wavefunction method such as DMRG.

Beside the optimization of the selection of the active space orbitals, a critical issue to take the dynamical correlation effects into account. The present Hartree–Fock treatment of the core orbitals [58] could be insufficient for many materials to achieve accurate electronic structure. In order to achieve an extremely accurate total energy of crystals, a coupled-cluster-single-double with perturbative triple [CCSD(T)] wavefunction approach is required [67]. It is likely that accurate low-energy excitation spectrum around the Fermi-level may be achieved at much lower complexity of wavefunction approaches than CCSD(T). One possible route of development is the implementation of more complex levels of wavefunction approaches than Hartree–Fock core and full CI active space model in the DMRG multiscale method, in order to converge towards highly accurate low-energy excitation spectrum. Another important issue, similarly to the DFT + CI methods, is to compute quantum mechanical forces acting on ions. The concept of the force calculation does exist for DMRG method that was already implemented for quantum chemistry codes (see Ref. [68] and references therein).

3.1.3 Spin-flipping time-dependent density functional theory and BSE methods

Time-dependent DFT (TDDFT) based on KS-DFT in the kernel [69, 70] can be principally applied to calculate the low energy excitation spectrum of quantum defects [71, 72]. In order to achieve accurate results, the proper choice of the DFT functional is essential [71]. Unlike the present implementations of GW + BSE and post-DFT multiscale methods, TDDFT framework and implementations exist to calculate the quantum mechanical forces acting on ions in the electronic excited state. In a pioneer work it was shown (see Supplementary Materials in Ref. [73]) that the observed Stokes-shift of NV centre can be well reproduced by TDDFT calculation when the NV centre is embedded in the core of 1.4-nm sized nanodiamond. The optimized geometry by TDDFT method well reproduced the optimized geometry by the ΔSCF method [74] for diamond NV centre. This issue has been recently investigated for the supercell model of diamond NV centre and silicon carbide divacancy centres [75]. It was found that the optimized geometries, the adiabatic potential energy surfaces (APES) of the 3E state and the zero-phonon-line (ZPL) energies are very close to each other as obtained ΔSCF method and TDDFT method based on DDH functional (see Figure 4 in Ref. [75]). The reason behind this good agreement is that Λ ≈ 0 for 3E so it is a very good approximation to describe 3E state as promotion a single electron from a1 in-gap defect level to the e in-gap defect level in the spin minority channel which is exactly constructed by ΔSCF method.

In the recent years, spin-flipping TDDFT (sf-TDDFT) theory has been developed and applied to molecules [76, 77] which is based on the Casida equations [70] but the original equations are modified in order to calculate the states and energies associated with Δm s = ±1 spin transition. In principally, the supercell implementation of the sf-TDDFT can be applied to diamond NV centre to obtain the singlet 1A1 and 1E states with geometry optimization in these states that has not yet been achieved by other means. In the DFT + CI-cRPA calculations [41], it was already recognized for diamond NV centre that the 1E state contains Slater-determinants associated with symmetry breaking solutions, thus 1E state is dynamically distorted from the C3v high symmetry geometry unlike the 1A1 state which stays in C3v symmetry. Insights from group theory, ΔSCF and DFT + CI-cRPA calculations with electron-phonon Hamiltonian models made it possible to construct the absorption and emission spectra between the 1A1 and 1E states for NV centre [29]. Nevertheless, the construction of the absorption spectrum was not accurate as principally it could not use the true APES of 1E and 1A1 states. As a consequence, the sharp resonance in the absorption band at 170 meV above the diamond phonon bands was missing in the constructed absorption spectrum.

In a seminal work by Jin and co-workers [78], the singlet states of diamond NV centre were directly calculated by sf-TDDT including quantum mechanical forces. The ab initio APES could be calculated both for the 1E and 1A1. Although, the sf-TDDFT excitation does not involve the double excitation electronic configurations in 1E and 1A1 states but their contributions may influence their energy levels – higher than that by a DFT + CI method [50, 78] – but not significantly their optimized geometries. The previously developed model Hamiltonians based on ΔSCF calculations [29] have been confirmed by sf-TDDFT calculations [78]: the 1E state is strongly anharmonic whereas the 1A1 state shows an almost perfect parabolic APES but the effective phonon frequencies are higher than that for 3A2 ground state. On the other hand, the 1A1 state shows a slight anharmonicity due to its phonon coupling to the 1E state by the symmetry breaking e phonons. Taking this correction into account, the calculated absorption spectrum shows a perfect agreement with the observed absorption spectrum including the sharp vibration resonance at 170 meV [78]. It was found that the e phonons dominantly contribute to the absorption spectrum, in stark contrast to the optical spectra between the triplet states of the defect. The luminescence spectrum between singlets was not calculated in this study [78] which is a Herzberg–Teller optical transition [29], and the highly accurate calculation of the shape of the phonon sideband would need to solve the multi-mode Jahn–Teller problem [79].

Spin-flip BSE (sf-BSE) method can be basically also applied to calculate Λ ≫ 0 states. Here the original idea is that |ΨGS⟩ of the system could be a close shell singlet but the low-energy singlet excited state might be a aiajakal|ΨGS type. This may be addressed by a high-spin reference state, e.g., aiaj|ΨGS type of S = 1 shelving state, by flipping the spin state by m s = −1 to reach |ΨGS⟩ (negative excitation energy) or by raising with m s = +1 to reach the excited state of aiajakal|ΨGS type (positive excitation energy). The sf-BSE spectrum can be calculated by ignoring the exchange terms in the BSE kernel. This methodology was first applied to atoms and molecules [80]. Parallel to this effort, Barker and Strubbe applied this method to diamond NV centre [81] with using PBE DFT functional as implemented in the Octopus code [82, 83]. By choosing S = 1 reference state raises the issue of spin-polarization within KS-DFT approach which naturally results in spin contamination error and an “artificial” gap between the occupied an unoccupied GWQP levels. Barker and Strubbe rather took the KS DFT energy levels instead of the GWQP levels in the BSE kernel as a workaround because PBE KS levels do not show this problem owing to the semilocal nature of the PBE functional. They picked up the 3A2 state as a reference state which is the ground state of diamond NV centre. The computed spectrum basically agreed with the previous GW + BSE result [37], i.e., the 1A1 level lies too low in the spectrum [81]. In order to calculate the 1A1 level correctly with sf-BSE method, the 3E reference state may be chosen which can access both single and double excitation Slater-determinants with respect to 3A2 ground state. Even though sf-BSE may work out for the singlet states accurately, the calculation of quantum mechanical forces should be implemented within GW + BSE and GW + sf-BSE methods to calculate the ZPL and Debye–Waller factor of the optical excitation spectrum of the defect qubits.

3.1.4 Temperature broadening of optical excitation spectrum with BSE method

The many-body perturbation theory of electron-phonon coupled optical transition with non-equlibrium Green-functions has been developed in Ref. [84] which was further derived in Ref. [85]. They considered the adiabatic limit for the dipole matrix elements, while they retained dynamical effects only in transition energies [85]. In terms of the phonon-dependent optical dipole transition moments, it goes beyond the Huang–Rhys theory. The central equation is

(10)IBSE(ω,T)λ,ν2|Πλ|2xν2fλ<δ(ωEλων)nB(ων,T)2ων+δ(ωEλ+ων)1+nB(ων,T)2ων,

where Π λ are the exciton (λ) dipole matrix elements within BSE theory, and fλ< are their occupations with E λ energy; x ν is the normal coordinate of phonon μ with the energy of ω ν . Here, fλ< is non-vanishing only if the excitons are composed by transitions between bands occupied by excited electrons and holes, and the two Dirac δ in the square bracket correspond to the cases where an exciton recombines with the creation [δ(ωE λ + ω ν )] or annihilation [δ(ωE λ ω ν )] of a phonon; n B is the Bose–Einstein occupation function for phonon ν at temperature T. The no-phonon optical transition is given by

(11)I0BSE=λ|Πλ|2fλ<δ(ωEλ).

In the usual implementation of BSE, it does not take into account the polaron shift. In other words, the no-phonon or zero-phonon line energy is calculated at the fix coordinate of the ground state from which the GW + BSE calculation starts. The theory was first applied to bulk hexagonal boron nitride which has indirect band gap, so I0BSE=0 (see Ref. [85]). Libbi and co-workers implemented the theory and applied it to the negatively charged boron-vacancy in hexagonal boron nitride [86]; the importance of defect was already mentioned in this review paper.

The observed fluorescence comes from a first-order forbidden transition between the 3E′′ excited state and 3A2 ground state [62] which becomes only allowed by participation of phonons. It was suggested based on ΔSCF HSE06 calculations that symmetry distorting Jahn–Teller distortions could lead to optical transitions. While it was noted based on the comparison DFT HSE06 and DMRG methods that HSE06 ΔSCF method has a limitation in describing the APES of 3E′′ excited state, still the Huang–Rhys fluorescence spectrum was calculated with relatively good agreement with the overall fluorescence energy but the calculated phonon sideband seemed to be too wide when compared to experimental data (see Ref. [62] and references therein). The in-plane Jahn–Teller distortion resulted in about 193 meV energy gain compared to that of the high D3h configuration. Further symmetry reduction and energy gain of 13 meV were also observed by HSE06 ΔSCF calculations by out-of-plane phonon modes [62]. The out-of-plane phonon modes are transformed as A2 and E′′.

Libbi and co-workers rather applied G0W0 method on PBE DFT calculations [86]. Then they calculated the non-equlibrium BSE optical spectrum based on Eq. (10). Because of the forbidden nature of the optical transition, the equilibrium BSE results in exactly zero optical dipole moment for the lowest energy transition in agreement with the previous HSE06 DFT result [62]. However, it becomes visible by applying non-equilibrium BSE method (I0BSE) at room temperature. The non-equilibrium BSE spectrum appears via fλ< in I0BSE through the thermalisation of electrons and holes. The pseudo-equilibrium occupations equal to

(12)fnk=1eϵnkμelkBT+1,f̄nk=1eϵnkμhokBT+1

for electrons and holes, respectively. Here ϵ n k is the quasiparticle energy of the state {nk} and μel (μho) the chemical potential for electrons (holes). The chemical potential was set for the electrons in such a way that a whole electron is promoted to the excited state manifold as usual in ΔSCF procedure which guarantees a neutral excitation. The chemical potential of the holes is tuned in such a way that the number of holes coincides with that of the electrons excited to the empty states. The non-equilibrium (NEQ) occupations induce a renormalization of the quasi-particle energy levels. It is calculated as

(13)ϵNEQ=ϵEQG0W0+ϵNEQCOHSEXϵEQCOHSEX,

where the first term is the quasi-particle energy determined at the G0W0 level of theory using the equilibrium occupations, while the second and the third terms represent the quasi-particle energy at the COHSEX level of theory calculated using the non-equilibrium and equilibrium occupations, respectively.

The calculated IBSE(ω, T) spectrum indicates that the phonons associated with the in-plane Jahn–Teller distortion has a very small optical dipole transition moments but the out-of-plane phonons significantly amplify the optical dipole transition moments [86]. As a consequence, the shape of the PL spectrum is governed by the out-of-plane phonons and not the in-plane phonons. We note that the combination of the 3E′′ electronic state with A2 and E′′ phonons results in Ẽ polaronic state which has allowed optical transition towards A2 ground state for each polarization of the emitted light. The overall width and shape of the calculated fluorescence spectrum agreed well with the experimental spectrum (see Ref. [86] and references therein) but the calculated spectrum is shifted to lower energy (of about 150 meV). The origin of the discrepancy was not explained in Ref. [86] whether it comes from the G0W0 approach or the neglect of polaron shift. Nevertheless, the methodology demonstrated its strength in analysing the Herzberg–Teller type optical transitions which is also dominant in the fluorescence spectrum of the singlet states in diamond NV centre [29].

3.2 Photo-ionisation thresholds: many-body and temperature effects

Photo-ionisation threshold energies and cross-sections are key properties for such quantum defects for which the qubits are initialised and read out optically. Photo-ionisation may promote an electron from the filled in-gap state to the conduction band (positive ionisation) or an electron from the valence band to the empty in-gap state (negative ionisation). For a given defect in semiconductors and insulators, both process may occur depending on the illumination wavelength. Similar to the neutral photo-excitation processes, phonons may assist the photo-ionisation processes with yielding temperature-dependent photo-ionisation thresholds. Ab initio simulations using the usual Born-Oppenheimer approximation with separating the problem of motion of ions and electrons will result in the low-temperature photo-ionisation thresholds that may be not accurate at elevated temperatures.

Photo-ionisation may occur by simply single photon absorption which basially goes the same way as neutral photo-excitation just the initial |Ψ i ⟩ or the final state |Ψ f ⟩ is the host band in the photo-ionisation process, i.e., Ψi|er̂|Ψf which is a one-body operator. Another possibility is the Auger-process for such defects that have multiple levels occupied by electrons such as the diamond NV centre (see Figure 1). In semiconductors, the Auger-process is considered to be important for the cases when the number of carriers is high in the band edges of the host material. However, in wide band gap materials, the Auger-process can well compete and be even dominating one when compared to other non-radiative processes such as the electron-phonon coupling because the large energy spacing between the defect level and the band edges leads to too slow multi-phonon processes. The critical matrix element is ϕiϕj|v̂eff|ϕkϕl where v̂eff is the screened Coulomb interaction appeared also in Eq. (2) where the many-body Ψ wavefunctions are now expressed by the appropriate single-particle wavefunctions ϕs (see Ref. [30] and references therein). The Auger-process is described by this two-body operator.

Another important consideration about describing the photo-ionisation is to distinguish the observable many-body picture and the auxiliary single-particle picture or band structure (see Figure 1). The band structure diagramme is an effective single particle picture including the localised defect states inside the band gap. In other words, band structure virtually plots the single-particle levels of which single-particle wavefunctions build up the many-body Slater-determinant solution. In the case of shallow donor state such as phosphorus donor (P-donor) in Si, the defect level is occupied by a single electron in the band gap, and that electron has a relatively weak exchange-correlation interaction towards the valence band states. Thus, the single electron orbital represents well the many-body total energy with respect to the total energy of the system when the electron is promoted to a higher energy effective mass state. As a consequence, the occupied donor defect level with respect to the conduction band minimum (CBM) can be considered as the photo-ionisation energy (when the relaxation of the ions upon photo-ionisation is neglected or minor). However, if the donor level is occupied by two electrons, e.g., sulfur donor in Si, then the two electrons have considerable Coloumb interaction with each other. As a consequence, once an electron is removed from the doubly occupied donor level then we break the Coulomb repulsion between the two electrons and their contribution to the total energy of the system before ionisation which is missing for the single electron left in the donor level after ionisation. Therefore, the difference betwen the doubly occupied donor level and singly occupied donor level does not correspond to the many-body total energy difference of the neutral and singly ionised sulfur defects. For deep defects with multiple electrons localised on the defect this effect is severe: it results in a strong Coulomb and exchange interactions among the localised electrons. Therefore, the many-body 3A2 level, of which state is ΨGS Slater-determinant for the diamond NV centre, should not be drawn to the band structure of perfect diamond as 3A2 state already contains all the diamond bands. The false level diagramme with mixing the single particle and the many-body pictures will lead to a false impression about the ionisation energies. For instance, if the 3A2 level is drawn at the position of the ionisation energy (about 2.65 eV at room temperature [87, 88]) w.r.t. CBM then the 3E level is usually drawn above the 3A2 level with the ZPL energy at 1.945 eV. According to this false consideration, the ionisation energy of 3E level w.r.t. CBM would be that of the 3A2 minus the ZPL energy which yields about 0.7 eV. The problem with this consideration is that it completely neglects the exchange-correlation effects between the electrons or other words the many-body effects as will be explained below.

In the example of diamond NV(−) → NV(0) ionisation process, if we start from the ground state 3A2 electronic configuration of NV(−) (Figure 2(b)) then a single electron from the e level is promoted to the CBM which results in the 2E ground state of NV(0). Starting the ionisation from 3E state of NV(−), we first have to consider the 3E state which can be well described as a hole is left in the a1 level and the double degenerate e level is occupied by three electrons (see left panel in Figure 2(b)). Two ionisation processes are viable at this point: (i) direct process with an electron promoted from the e level to the conduction band which may leave the defect [89]; (ii) Auger-process occurs after absorption of the second photon where the promoted electron in the conduction band recombines with the a1 hole, and the energy gain of this process is used to simultaneously promote another electron from the e level to a high energy state in the conduction band which finalizes the ionisation process (see Refs. [90, 91] and Figure 2(c)). According to the Slater–Condon rules [92, 93], only one spin orbital can change in the many-body wavefunction upon direct ionisation process described by one-body operator, therefore 3E of NV(−) arrives at the metastable 4A2 of NV(0) plus an electron in the conduction band by photo-excitation. Since the shelving 4A2 level lies above the ground state 2E level of NV(0) the ionisation threshold energy of this process starting from 3E of NV(−) will be higher by about 0.48 eV according to HSE06 calculations [89] then the ionisation threshold energy starting from the 3A2 ground state level minus the ZPL energy of the 3E3A2 optical transition, which finally yields about 1.2 eV ionisation threshold energy [89]. In the alternative Auger-process, the two-body operator nature of the process makes it possible to arrive at the 2E ground state of NV(0) from 3E NV(−). If the electron is excited to high energy above the CBM then phonons can very quickly (within picoseconds) cool it down to the CBM which is might be bound by the weakly attractive potential of NV(0) towards the CBM electrons which makes this process viable with about 1 nanosecond rate [90]. In the two-photon absorption process of NV(−), the typical excitation energy is about 4.66 eV (green light at 532-nm wavelength) which subsequently goes through the 3E level with phonon cooling upon absorption of the first photon and probably a phonon cooling process after absorbing the second photon as explained above. From energetics point of view, it is feasible to arrive at 2E ground state of NV(0) from the 3A2 ground state of NV(−). A recent experimental study concludes that the two-photon absorption based charge conversion of NV(−) → NV(0) can be well explained by a dominating Auger-process. In the rate modelling, a single-photon optical transition should also occur from the shelving 1E state to a higher lying state [91]. A strong and broad transition around 2.58 eV from 1E to 1E′ has been observed in numerical simulations providing a possible candidate for such a mechanism [41]. 1E′ can decay to the triplet excited states via an inter-system crossing [30, 94].

Figure 2: 
Photo-excitation and ionisation of diamond NV centre. (a) Calculated cross section as a function of photon energy [89]. Solid blue: photo-ionisation from the excited state 3E, σph; dark red: stimulated emission, σst; orange: intradefect absorption, σintra; dashed blue: photo-ionisation from the singlet state 1E. Photo-ionisation thresholds from 3E and 1E are indicated (estimated error bar 0.1 eV), together with the experimental values of the ZPL energy for NV(−) and NV(0). (b) Kohn–Sham (single particle) states of the NV(−). The dashed arrows show the excitation processes from the 3A2 ground state (left panel) and 3E excited state (right panel), where both the stimulated emission and photo-ionisation processes are depicted for the latter. 4A2 shelving state of NV(0) is left after completing the ionisation whereas stimulated emission brings the system back to 3A2 ground state of NV(−). (c) Photo-ionisation via Auger-process after the second photon was absorped. 2E ground state of NV(0) is left after completing the ionisation.
Figure 2:

Photo-excitation and ionisation of diamond NV centre. (a) Calculated cross section as a function of photon energy [89]. Solid blue: photo-ionisation from the excited state 3E, σph; dark red: stimulated emission, σst; orange: intradefect absorption, σintra; dashed blue: photo-ionisation from the singlet state 1E. Photo-ionisation thresholds from 3E and 1E are indicated (estimated error bar 0.1 eV), together with the experimental values of the ZPL energy for NV(−) and NV(0). (b) Kohn–Sham (single particle) states of the NV(−). The dashed arrows show the excitation processes from the 3A2 ground state (left panel) and 3E excited state (right panel), where both the stimulated emission and photo-ionisation processes are depicted for the latter. 4A2 shelving state of NV(0) is left after completing the ionisation whereas stimulated emission brings the system back to 3A2 ground state of NV(−). (c) Photo-ionisation via Auger-process after the second photon was absorped. 2E ground state of NV(0) is left after completing the ionisation.

Nevertheless, such spin-to-charge-conversion (SCC) protocols exist for diamond NV centre which apply low-power excitation to avoid two-photon absorption processes but dual wavelengths of excitation in which the green illumination is used for the optical cycle between 3A2 and 3E states whereas longer-wavelength-than-ZPL illumination is applied to ionise it from the 3E level towards the conduction band [95, 96]. In this case, direct ionisation may occur from the 3E level. In the seminal work of Razinkovas and co-workers [89], the absolute ionisation cross section from the 3E state together with the induced emission was calculated for NV centre as a function of the excitation energy (Figure 2(a)). The calculated ratio of the photo-ionisation cross section and the cross section for stimulated emission is greater than 2 for the energy around 1.2 eV and 1.93 eV photo-excitation energies, which agree well with the applied photo-ionisation energies in the SCC experiments [95, 96].

We note that understanding the reionisation process, NV(0) → NV(−), is critical for stabilising the diamond NV(−) qubit state. Typically, the reionisation automatically occurs by the applied green illumination to drive NV(−) by two-photon absorption process which is also subsequential process going through the 2A2 excited state: the first photon absorption brings 2E to 2A2 and the second photon is then absorbed in the 2A2 excited state. The 2A2 state consists of a hole left in the a1 level and two electrons in the e level. The photo-excitation of the 2A2 excited state of NV(0) may also occur by either direct process (promoting an electron from the valence band to the empty a1 defect level in the gap) or Auger-process (occupying the in-gap a1 hole by an electron from the in-gap e level and then promoting an electron from the valence band to the empty e defect level in the gap). Both processes leave a hole in the valence band. In the direct process, the system arrives at the ground state of NV(−) because of the alluded Slater–Condon principle [91]. On the other hand, the Auger-process enables to arrive at the 1A1 and 1E singlet states of NV(−) too, beside the 3A2 ground state. The energy cost of these processes varies with the final state. The calculated adiabatic acceptor charge transition level of the NV defect is at about 2.75 eV from the conduction band edge [97, 98], whereas the calculated energy gap between the 3A2 ground state and 1A1 state is at about 1.6 eV (see Ref. [30] and references therein). The total energy cost to convert the NV(0) ground state to the 1A1 NV(−) excited state is then about 4.3 eV which coincides with twice the ZPL energy of NV(0) [91]. This means that a special excited state of 1A1 of NV(−) binding a hole resonant with the valence band maximum (VBM) develops. This hole is Coulombically bound which is a special bound exciton state or Rydberg state which has been observed for the SiV defect [11] and has been recently implied and modelled for the 3A2 plus bound hole system for the NV(−) defect [99, 100]. The bound hole is weakly localised with following the effective mass theory. By even taking into account the possible relaxation energy of the ions caused by the change in the electronic states, we may claim that a green laser excitation can reach the 1A1 plus bound hole state of NV(−) by two-photon excitation of NV(0). Scattering rather to 1E and 3A2 states of NV(−) via Auger-process leaves a hole deep in the valence band at around 1.2 eV and 1.6 eV from the VBM, respectively. According to the calculations [91], a resonant a1 state broadened by the diamond bands occurs in this energy region which is originated from the dangling bond orbitals of the carbon and nitrogen atoms near the vacant site. The resonant state is weakly localized unlike the usual diamond bands that are completely delocalized. This should lead to a larger direct and Auger-ionization rates of NV(0) than those of NV(−) because no such a high-energy resonant state sharing the same spin state with that of the ground state exists in the conduction band, critical in the photo-ionisation of NV(−). Interestingly, the Auger-process should lead to a preferential occupation of the m s = 0 state via spin-selective intersystem crossings between the 1E state and the 3A2 spin states; however, direct ionisation would result in 1/3:2/3 relative population of the m s = 0:m s = ±1 states of 3A2 [91]. The ab initio calculation of the rates of the direct and Auger-processes requires an accurate computation of the excited states of NV(0) which has not yet been solved as we briefly discussed it in Section 3.1.

The afore-mentioned ab initio calculations are based on the global energy minimum of the APES in the appropriate electronic structures and charge states which correspond to the zero kelvin solution. However, the photo-ionisation tresholds of diamond NV centre are often observed at room temperature. Principally, the effective ionisation threshold energies may change as a function of temperature. A very characteristic example is the silicon carbide (SiC) divacancy defects. In particular, four divacancy configurations with similar electronic structures occur in the so-called 4H polytype of SiC [101103], some of them well observable at room temperature [104, 105]. We note that 4H SiC exhibits a band gap of about 3.3 eV that can host visible and near-infrared colour centres acting as qubits [11, 106, 107]. The neutral divacancy defects in 4H SiC have isovalent electronic structure to diamond NV centre as depicted in Figure 1 but the energy gaps between the levels are about twice as smaller, thus they produce near-inrared emission. It was found that upon photo-excitation of the defect it falls to a “dark” shelving state due to two-photon absorption or other complex processes. This shelving state has been finally identified as the negative charge state of the defect [108110]. Close to cryogenic temperatures, a shorter wavelength laser beam was applied to drive the divacancy back to the neutral qubit state by promoting an electron from the e level to the CBM [110]. It was later found that at elevated temperature the quenching of the fluorescence of 4H SiC divacancy (V2) defects does not occur and they remain optically stable [111]. This can be interpreted as the photo-ionisation threshold energy associated with the V2(−) → V2(0) process is decreasing upon raising the temperature.

The modelling of temperature-dependent photo-ionisation processes requires the following considerations: (i) the CBM and VBM of the host material may shift with temperature, (ii) the formation enthalpy so the charge transition level of the defect with respect to the band edges may shift with temperature via the vibration entropy, and (iii) the phonon assisted ionisation process may be activated by raising the temperature. The first effect is basically the temperature-dependent electron-phonon renormalization of the bands. This can be computed by the many-body perturbation theory on the electron-phonon coupling [112114] that was applied to 4H SiC [115]. The usual case is that the CBM and VBM shifts down and up with raising the temperature, respectively, leading to an effective decrease of the fundamental band gap. It was found by ab initio calculations that the CBM of 4H SiC shifts down by about 5 meV from zero kelvin to room temperature [115]. In other materials with low Debye-temperature, this effect could be significantly enhanced. The second effect assumes an ab initio treatment of the thermodynamic properties of solid-state defects that was thouroughly discussed in Ref. [34]. Here, the key effect is the vibration entropy correction to the formation energy of the defect,

(14)Fq(T)=i12ωi+kBTln1expωikBT,

where and kB are the reduced Planck constant and Boltzmann constant, respectively, and ω i is the frequency of the ith phonon mode in charge state q of the defect at the given T temperature. The first term in Eq. (14) is the zero-point energy. The actual values of F q (T) could differ for a given defect in various charge states (ΔF(T) = F(T) − F0(T)) which results in a shift of charge transition level, here with respect to CBM, E/0CBM. This correction may increase or decrease the effective photo-ionisation threshold. The total correction is then,

(15)E/0CBM(corr)=E/0CBMΔECBM(T)+ΔF(T).

The value of E/0CBM(corr) corresponds to the ionisation threshold energy without involvement of phonons in the photo-ionisation process [111]. By raising the temperature the phonon excited states are occupied and may contribute to the ionisation process which may be considered as an phonon-assisted optical transition between the neutral defect binding an electron in the CBM and the negatively charged defect where we neglect the interaction of the CBM elecron with the rest of the electrons. This effect may strongly contribute to the reduction of the photo-ionisation threshold at elevated temperatures. The effect is analogous with the appearance of the phonon bands at lower energies than ZPL energy in the absorption spectrum of the defect at elevated temperatures. The phonon-assisted optical spectrum can be calculated within Franck–Condon approximation at ab initio level for defective supercell models (see Refs. [30, 116] and references therein). The calculated temperature-dependent ionisation threshold energies for the divacancy centres in 4H SiC based on this theory were reported in Ref. [111]. Although, the very long wavelength acoustic phonons were not calculated in that study for the vibration entropy and the Franck–Condon theory that shed doubts on the convergence of the results, nevertheless, the equal or similar contributions of the temperature-dependent vibration entropy and Franck–Condon terms are demonstrated, and it provided explanation about the photo-stability of 4H SiC divacancy qubits at room temperature [111].

3.3 Supercell modelling of defects: extrapolation to dilute limits

In the supercell modelling of point defects the goal is to describe isolated point defects. In practice, the size of the supercell is limited up to about 10,000 atoms for KS DFT calculations due to computational resources. By applying accurate hybrid KS DFT functionals the size of the supercell is further reduced to about 1000 atoms. This size of the supercell suffices to obtain accurate results for deep defects because the defect induced wavefunctions and strain fields decay relatively fast from the core of the defect. However, shallow states with weakly localised character of defect wavefunctions such as the case of P-donor in Si, require extrapolation to dilute limits. For deep defects, certain properties also call for special treatments in the supercell formalism that we non-exhaustively list here: charge correction for charged defective supercells, acoustic phonon couplings to electron orbitals and spins, excitation and ionisation towards electronic bands of the host.

3.3.1 Charge correction schemes: a recent breakthrough

In our previous review paper [30], the charge correction for charged defective supercells was thouroughly discussed. To sketch here the problem and the possible solution we note that the introduced defect charge is neutralized by a compensating jellium charge in the supercell. The charge of the defect and the jellium background charge interact with their periodic images that goes with the leading point charge Coulomb interaction of which energy scales with L−1, where L is the edge of the simple cubic supercell. This theory already gives a hint about the expected scaling property in the correction of the total energy to the isolated defect or dilute limit. Indeed, the total energy of the charged defective supercells converges notoriously slow with the supercell size as indicated by the simple theory above. By applying charge corrections to the total energy, the convergence can be well accelerated. In the previous years (see Refs. [30, 34, 117] and references therein), many a posteriori schemes have been developed for the correction of the total energy of charged defective supercells in 3D and 2D models. In those schemes, the possible change in the character of the wavefunction due to the charge correction was not taken into account that might lead to qualitatively wrong results in notorious cases. This problem is in particular severe for slab models of crystal surfaces with negative electron affinities. Recently, a self-consistent potential charge correction (SCPC) method was developed to heal this issue which goes beyond the previous a posteriori total energy corrections of the charged defects, and they derive the KS potential associated with the charge correction and self-consistently solve the constructed KS DFT equations [117].

The SCPC method yield the corrective potential (Vcor) in an iterative manner: (i) the distribution of the extra charge in the supercell (δρ) is determined, (ii) the corresponding periodic electrostatic potential (Vper) is calculated, (iii) the potential for the same but isolated charge distribution (Viso) is determined by using open (Dirichlet) boundary conditions, and finally, (iv) Vper and Viso are used to determine the corrective potential Vcor, which is added to the total electronic potential. It should be noted that SCPC always aligns the final potential, considering the difference between the electrostatic potentials of the charged and the reference system far away from the defect position (ΔV).

The method was originally built in the VASP code but an interface has been developed to Quantum Espresso code too. The SCPC method was applied to diamond slab model with NV centre where the (100) diamond surface was terminated by hydrogen. It was shown that without SCPC method the negative charge of the defect artificially pull down the bands of the surface states, so-called image states, which results in a false electronic structure even in large supercells [117]. Here, the self-consistent correction of the potential is essential. We note that the self-consistent correction is not a must for many defects in 3D solids and a posteriori charge correction schemes can provide qualitatively good results. Despite the self-consistent nature of the correction, the SCPC method does show supercell size dependence which comes from the fact that the character of the underlying (localized) defect wavefunctions may change with supercell size due to the defect-induced strain fields and other factors, so supercell size scaling is still necessary with this method but converges faster than the a posteriori total energy correction schemes [117].

3.3.2 Embedding of long wavelength phonons in the finite size of the supercell

The long-wavelength acoustic phonons could contribute to the phonon sideband of the optical spectra of defects. Audrius Alkauskas and co-workers developed an embedding method to include the electron-phonon coupling in the optical transition [116] that they applied to the optical transition between the triplet states of diamond NV centre [79, 116] which has been also implemented in the sf-TDDFT study for the absorption spectrum of the singlet states [78]. Similar treatment is advisable for the temperature-dependent photo-ionisation spectrum [111].

3.3.3 Reconstruction of the deep-energy valence and high-energy conduction bands

Accurate absorption cross section calculation of excitation from/to deep levels to/from solid-state bands within supercell formalism requires special attention due to band folding in the Brillouin-zone [118]. This was recognized by Razinkovas and co-workers when they studied the absorption cross section for the photo-ionisation of diamond NV centre [89]: they found minigaps in the conduction band which affected the calculated absorption cross section at the given high energy excitation. They used the following technique to circumvent this problem: (i) they identified the folded band in the Brillouin-zone of the primitive diamond cell as described in Ref. [118], (ii) then they averaged out the calculated absorption cross section values closest to the corresponding k-space within the energy region of about 0.08 eV. The resulted photo-ionisation absorption spectrum is then converged well.

3.3.4 Treatment of spatially extended defect wavefunctions: beyond effective mass theory

KS DFT calculation of the properties of shallow donors in Si is computationally very challenging because of the spatial extension of the donor wavefunction. Accurate calculation requires hybrid HSE06 functional, e.g., Ref. [119], which becomes extremely difficult to carry out for sufficiently large supercell sizes. As was done for the total energy correction of charged defective supercells (see Refs. [30, 34, 117] and references therein), one can apply a strategy by studying theoretically and numerically the scaling of the given property as a function of the supercell size, and then extrapolate the result to the dilute (isolated defect) limit. In many cases, the exact scaling properties are unknown or they are too much complex because different effects (electron density and charge distribution, strain field distribution) are intertwinned and they often depend on the local electron density distribution of the defect that might change by increasing the size of the supercell.

In practice, numerical KS DFT investigations could lead to converged results which requires sufficiently large number of sampling points for defining the scaling law, i.e., ideally up to supercell size with about 10,000 atoms. This is prohibited by the accurate HSE06 functional so a typical strategy is to calculate the scaling by the affordable PBE functional and for a limited range of supercell sizes to compare the scaling with HSE06 and PBE functionals to verify the scaling function.

Swift and co-workers studied the shallow group-V dopants of Si by supercell KS DFT calculations [119]. These dopant introduces a singly occupied electron the state split from the CBM. In semiconductor physics, these states are described by the so-called effective mass theory from Kohn and Luttinger which treats the dopant potential as a positively charged Coulomb potential which binds an electron of which state can be described the linear combination of the CBM valleys. The solution of this system results in a Rydberg or hydrogenic series of excitation energies until it converges to the ionisation level. The ionisation or binding energy of the electron (E b ) can be defined as the lowest energy level ϵdonor (1s like envelop function) with respect to CBM of the perfect crystal (ϵCBM). The validity of this approximation from many-body electron-phonon point of view was briefly discussed in Section 3.2. It is known in experiments (see Ref. [119] and references therein) that the ionisation energies of various group-V dopants in Si differ. Therefore, a so-called “central cell correction” was introduced to the effective mass theory which assumes that the 1s ground state wavefunction has the largest overlap with the dopant ion’s attractive potential which will pull down its energy level with respect to the purely hydrogenic solution (enlarging the donor ionisation energy). As the potential of the dopant ion is characteristic to the dopant within short range around the dopant ion, thus the “central cell correction” will be dopant dependent. The central cell correction is merely semiempirical correction to the effective mass theory of shallow donors and acceptors, where the acceptor levels are measured with respect to VBM. In KS DFT calculations, the ionisation energies can be calculated at ab initio level free from any assumptions on the nature of the potential induced by the dopants. It can be estimated from the effective mass theory that the 1s donor wavefunction in Si will decay at around 55 Å from the position of the dopant which would require a supercell of about 64,000 atoms to accommodate the donor wavefunction without significant overlap. This definitely calls to apply a scaling procedure even for PBE DFT functional.

Swift and co-workers [119] applied scaling method for calculating E b of arsenic (As) and bismuth (Bi) donors in Si by PBE and HSE06 functionals. The calculations were carried for supercells from 64 to 2744 atoms at PBE level and for supercells from 64 to 1000 atoms at HSE06 level. They calculated E b at a given size of the supercell as

(16)Eb=ϵCBMϵdonor+eΔV,

where ΔV potential alignment between the defective and perfect supercells appears the charge correction of defects in Section 3.3.1. Interestingly, the ΔV does not show a clear monotonous decay with increasing size of the supercell which was not explained [119]. The scaling of the HSE06 data was carried out as

(17)bHSE06=bPBEfit12bPBE,δexfit+12bHSE06,δexfit,

where bPBEfit is the slope fit to the PBE binding energies, bPBE,δexfit is the slope fit to the PBE exchange splitting, and bHSE06,δexfit is the slope fit to the HSE06 exchange splitting. The exchange splitting was defined as the difference between the spin-up and spin-down eigenvalues of the donor state. Subtracting half of the exchange splitting from the binding energy yields a spin-averaged value which is corrected between the PBE and HSE06 results with assuming that the supercell-size dependence obtained in PBE for the spin-averaged case applies also to the HSE06 values. They obtained 54 meV and 67 meV binding energies for As and Bi donors, respectively, in very good agreement with the experimental data at 53.9 meV and 70.9 meV, respectively [119].

The neutral donor defects introduce S = 1/2 spin that can interact with the nuclear spins of the dopant or the proximate 29Si I = 1/2 nuclear spins which is called hyperfine interaction. Actually, the interplay between the electron and nuclear spins could represent qubit states in these systems, therefore, understanding this interaction as a function of electric field and strain is highly important (see Ref. [120] and references therein). This hyperfine interaction can be generally written as

(18)Hhyp=ŜAÎ,

where A is the hyperfine tensor and Ŝ, Î are the electron spin, nuclear spin vector operators, respectively. The Fermi-contact and the dipole-dipole terms of the hyperfine tensor can be respectively written as

(19)Aab(n)=2μ03geμBgnμnnsRS+μ04πgeμBgnμn1S×3rarbr2δabr5nsrd3r,

where nsr is the electron spin density, r is the vector between the electron spin and nuclear spin at R, gn is the nuclear g-factor, and μn is the nuclear magneton for a given nucleus n. For the ground state 1s donor wavefunction, the Fermi-contact term predominates which depends on the localization of the donor wavefunction at place of the dopant. Swift and co-workers found that (i) PBE produces too low hyperfine constants for shallow donors in Si, (ii) reliable values are obtain for supercell size of 512 atoms or larger number of atoms, and (iii) in that size range the hyperfine constant for the dopant atom scales as L−1. The final HSE06 values are 132.5 MHz and 1262 MHz for As and Bi spins, respectively, compared to the experimental data at 198.3 MHz and 1475 MHz, respectively [119]. They also found that the electric field gradient, V zz , around the dopant atom, so the quadrupole interaction strength C Q = 3eQatomV zz /4h is vanishingly small (1 MHz), where h and e are the Planck-constant and the charge of the electron, and Qatom is the nuclear electric quadrupole moment of the dopant atom. Furthermore, they also studied the strain dependence of the contact hyperfine tensor of the Bi dopant, which generally reads [120] as

(20)A/A0=1+K3(εxx+εyy+εzz)+L2(εyyεzz)2+(εxxεzz)2+(εxxεyy)2+Nεyz2+εxz2+εxy2,

where A0 is the value in the absence of strain, K coupling is responsible for the hydrostatic strain, L coupling and N coupling describes the uniaxial and shear strain effects, respectively. Swift and co-workers only focused to the hydrostatic strain in their ab initio study: they found that K scales the same in HSE06 and PBE functionals in relatively small supercells, thus the error in PBE in the difference of absolute values of hyperfine constants as a function of strain is cancelled. As a consequence, the PBE scaling can be applied for larger supercell for extrapolation to the dilute limit. They obtained K = 20.2 which is close to the experimental data at K = 19.1 (see Refs. [119, 120]).

In the afore-mentioned examples, the defect wavefunction is spatially extended in the electronic ground state. Perhaps, it is not a common knowledge among scientists coming from the quantum optics field that similarly extended wavefunctions could exist in the electronic excited states. Again, the best example is the most studied small band gap material, silicon. As silicon has a band gap of 1.215 eV at cryogenic temperatures, there is a little room to introduce multiple levels by deep defects. Defects may introduce only a single occupied deep level in the fundamental band gap of Si where the electron could be promoted from the in-gap defect level to CBM. In this case, the defect can be described as a positively charged centre which Coulombically binds an electron split from CBM. This definitely shows a similarity to the shallow donor states in Si. For example, the photoluminescence C-centre in Si shows a sharp ZPL at 789 meV for which photoluminescence excitation (PLE) measurements revealed a hydrogenic or Rydberg series of excited states [121]. Later it was shown that this type of PLE features is common for other deep optical centres of Si that was called “pseudo-donor” model [122, 123]. Recently, the pseudo-donor model of C-centre in Si has been confirmed by HSE06 calculations [124]. The neutral C i O i defect associated with this optical centre indeed produces a deep level in the fundamental band gap (see Figure 3), and the calculated ZPL energy at 750 meV by ΔSCF method agreed well with the experiment. In this calculation, 512-atom supercell was employed with the same correction in the total energy of the excited state of the defect (56 meV) as for the positively charged defect. The reasoning behind this method was the following. The excited state involves a spatially extended wavefunction. The scaling property of the total energy of the excited state was assumed to go similarly to that of the positively charged defect within the accuracy of about 50 meV because the extended electron occupying the state split from CBM could behave similarly to all the crystalline valence bands of the system which leaves a positively charged core; in other words, the pseudo-donor electron does not “shield” the positively charged core. There is a further note here about the accuracy of ΔSCF method. In a 216-atom supercell, the results from ΔSCF method and GW + BSE method were compared. It was found that GW + BSE confirmed the composition of the exciton as the deep hole state and an electon state split from CBM and the vertical excitation energies were within 11 meV. Thus, ΔSCF method can be applied for the bound exciton excitation too which is important to calculate the Stokes-shift upon excitation as no quantum mechanical force calculation has been yet implemented to GW + BSE methods. For the case of shallow donors in Si, the geometry change upon ionisation was neglected. However, this cannot be neglected for the fluorescence spectrum of deep defects. Indeed, the sharp features in the phonon sideband of the PL spectrum could be well reproduced by applying the Franck–Condon theory [124]. According to the calculations, C-centre is a potential building block of quantum repeaters in the telecom L-band [124].

Figure 3: 
C-centre in silicon. (a) The atomic structure of the CiOi defect complex, consisting of neighbouring carbon (Ci − Si1) and oxygen (Oi − Si2) split-interstitial defects associated with C-centre in Si. (b) Kohn–Sham level structure of the CiOi(0) defect ground state (1A′) and singlet (1A′′) excited state in the spin-polarised HSE06 calculation. VBM, CBM and C label the valence band maximum, the conduction band minimum, and the dangling bond orbital of the carbon atom, respectively.
Figure 3:

C-centre in silicon. (a) The atomic structure of the CiOi defect complex, consisting of neighbouring carbon (Ci − Si1) and oxygen (Oi − Si2) split-interstitial defects associated with C-centre in Si. (b) Kohn–Sham level structure of the CiOi(0) defect ground state (1A′) and singlet (1A′′) excited state in the spin-polarised HSE06 calculation. VBM, CBM and C label the valence band maximum, the conduction band minimum, and the dangling bond orbital of the carbon atom, respectively.

Another deep optical centre in Si, the W-centre, has been recently isolated as a single quantum emitter with ZPL wavelength close to the telecom region at 1218 nm (1.018 eV) (see Ref. [20] and references therein). The defect contains a complex of three self-interstitial silicon atoms. The most stable configuration, so called I3-V configuration with C3 symmetry, has been recently identified by HSE06 calculations where the relative stability changes with respect to PBE calculations [20], which we call here I3 for the sake of simplicity. The electronic structure of the neutral I3 is very interesting: it shows a single resonant a level at 73 meV below VBM (see Figure 4). At first glance, this defect may be considered as electrically and optically inactive. However, after ionisation, the unoccupied defect level emerges inside the band-gap, and the (+/0) charge transition level is at 55 meV above VBM after applying charge correction in the total energy of the positively charged defect. As the stability of the positive charge state is confirmed, the positively charged defect may Coulombically bind an electron with the state split from CBM, alas, the neutral excitation of I3 is a bound exciton with a strongly localized hole on the defect and a loosely bound electron [20]. The pseudo-donor nature of the defect was confirmed by HSE06 ΔSCF calculation in 512-atom supercell. The estimation of the ZPL energy was based on the full geometry relaxation of the electronic ground and excited states with scaling of supercell sizes,

(21)EZPL(L)=A/L+B/L3+C,

where L is the side of the simple cubic supercell, A, B, C are fitting constants, where C value corresponds to the dilute limit. L was varied between the supercells of 216-atom and 8000-atom for PBE calculations and up to 1000-atom for HSE06 calculations. It was found that the 216-atom supercell results do not fit to the trend and should be ignored. The idea of the formula in Eq. (21) is that the excited state requires charge correction. Since too few data points could be calculated at HSE06 level, it was assumed that the PBE results well reproduce the electrostatics of the problem, and the A and B fit results can be used for HSE06 data points. This procedure finally yields C = 1.102 ± 0.003 eV which is within 0.1 eV when compared to experimental data.

Figure 4: 
W-centre in silicon. (a) The atomic structure of the tri-interstitial (I3) complex associated with W-centre in Si. (b) Kohn–Sham level structure of the I3(0) defect ground state (1A) and (1A) excited state in the spin-polarised HSE06 calculation. VBM and CBM label the valence band maximum and the conduction band minimum, respectively.
Figure 4:

W-centre in silicon. (a) The atomic structure of the tri-interstitial (I3) complex associated with W-centre in Si. (b) Kohn–Sham level structure of the I3(0) defect ground state (1A) and (1A) excited state in the spin-polarised HSE06 calculation. VBM and CBM label the valence band maximum and the conduction band minimum, respectively.

The pseudo-donor or bound exciton excitation can occur in wide band gap semiconductors too. A very nice example is the so-called D I centre in 4H SiC [125, 126]. The optical activity of the defect is identified as the silicon antisite [127, 128] which is an isovalent centre with producing a deep donor level in the fundemental band gap. The defect can be positively charged and then it can Coulombically bind an electron from CBM with producing Rydberg series of excited states [125]. It can be expected that similar bound exciton excited states may be found in diamond.

Indeed, a recent study has identified Rydberg series in the optical spectrum of the diamond neutral silicon-vacancy [SiV(0)] defect in a joint experimental and theoretical study [11]. Interestingly, optical spin-polarisation and ODMR signals could be also observed through their bound exciton states [11] which makes the analysis of these states highly important as this defect can be isolated in diamond as a near-infrared single photon emitter [129].

The electronic structure of SiV(0) in diamond was already briefly described in Section 3.1.1 that we extend here before we proceed to the discussion of the bound exciton states. In SiV(0) defect, Si atom sits in the inversion centre of diamond so the defect can be rather described as a V2 defect with six carbon danglings bonds whereas the “dopant” Si ion resides in the empty space of V2 with the farthest distance from these six carbon atoms (see Figure 5). The six carbon dangling bonds create a double degenerate e u level resonant with the valence band and a double degenerate in-gap e g level occupied by two eletrons. This forms the 3A2g ground state. The usual optical activity is associated with promoting an electron from the e u level to the e g level for which the optically allowed 3E u 3A2g transition yields the ZPL energy at 1.311 eV (see Figure 6(c)). This energy is much smaller than the ionisation threshold energy at 1.53 eV which corresponds to the neutral to the negative charge transition (see Ref. [11] and references therein). It is important to notice the selection rules of optical centre with inversion symmetry that the optical transition is only allowed by changing the parity of the participating wavefunctions.

Figure 5: 
Extension of effective mass state in diamond SiV(0) defect. (a) Geometry of the SiV(0) defect in diamond as optimized by HSE06 in the ground state. The dashed circles represent the missing C atoms, i.e., adjacent vacant sites or V2. (b) The Bohr-diameter of the heavy-hole for the n = 1 or 1s effective mass state. Apparantly, the dangling bond orbitals in V2 are confined to the same space as the 1s effective mass state.
Figure 5:

Extension of effective mass state in diamond SiV(0) defect. (a) Geometry of the SiV(0) defect in diamond as optimized by HSE06 in the ground state. The dashed circles represent the missing C atoms, i.e., adjacent vacant sites or V2. (b) The Bohr-diameter of the heavy-hole for the n = 1 or 1s effective mass state. Apparantly, the dangling bond orbitals in V2 are confined to the same space as the 1s effective mass state.

Figure 6: 
Experimental signals and effective mass theory for the excited state of diamond SiV(0) defect. (a) Experimental photoluminescence excitation (PLE) and optical spin-polarisation (OSP) spectra from Ref. [11]. The OSP can be observed in the spin-polarisation of the ground state towards the m
s
 = 0 electron spin state in the electron paramagnetic resonance spectrum. (b) Scaling of the peak positions extracted from PLE in (a). The fit uses Rydberg scaling E
n
 = E
I
 − E
y
/n2 associated with the effective mass states, where n refers to the principal quantum number. Due to similar fine structures of 2p and 3p states, we fit different fine structure transitions (wiggles in the PLE and OSP curves) separately corresponding to the different coloured curves. The fitted ionization energy (E
I
) and Rydberg energy (E
y
) are 1.53 eV and 0.4 eV, respectively. The horizontal dashed line indicates the fitted ionization energy. States with “s”-like character are taken from spin-polarisation measurements, and are shown with triangles. Δ1 and Δ2 are energy deviations for 1s and 2s states compared to the fitted Rydberg scaling that involve both central cell correction and the localized phonon energy. (c) Proposed bound exciton model for the higher-lying excited states showing orbital ground and excited states and BE states at higher energies in the hole picture. The lower levels closer to the valence band maximum for electrons require higher excitation energy.
Figure 6:

Experimental signals and effective mass theory for the excited state of diamond SiV(0) defect. (a) Experimental photoluminescence excitation (PLE) and optical spin-polarisation (OSP) spectra from Ref. [11]. The OSP can be observed in the spin-polarisation of the ground state towards the m s = 0 electron spin state in the electron paramagnetic resonance spectrum. (b) Scaling of the peak positions extracted from PLE in (a). The fit uses Rydberg scaling E n = E I E y /n2 associated with the effective mass states, where n refers to the principal quantum number. Due to similar fine structures of 2p and 3p states, we fit different fine structure transitions (wiggles in the PLE and OSP curves) separately corresponding to the different coloured curves. The fitted ionization energy (E I ) and Rydberg energy (E y ) are 1.53 eV and 0.4 eV, respectively. The horizontal dashed line indicates the fitted ionization energy. States with “s”-like character are taken from spin-polarisation measurements, and are shown with triangles. Δ1 and Δ2 are energy deviations for 1s and 2s states compared to the fitted Rydberg scaling that involve both central cell correction and the localized phonon energy. (c) Proposed bound exciton model for the higher-lying excited states showing orbital ground and excited states and BE states at higher energies in the hole picture. The lower levels closer to the valence band maximum for electrons require higher excitation energy.

By increasing the excitation energy above the ZPL at 1.311 eV but below the ionisation threshold energy at 1.53 eV, one can excite the hole from the VBM which results in a SiV(−) defect plus a loosely bound hole, i.e., a bound exciton state of SiV(0). Generally, analysing the hole bound exciton spectrum has an increasing complexity over that of electron bound exciton spectrum because of the orbital degeneracy of VBM at the Γ-point which results in an effective spin–orbit interaction. A detailed description is beyond the scope of the present review paper. We rather defer the readers to the supplemental material of Ref. [11] which is now very briefly summarized here. The three-fold degenerate VBM of diamond sligthly splits due to the defect potential resulting in a1g and e g bands, where a1g band lies above e g band. The VBM splitting also affected by the spin–orbit coupling which has similar energy as the crystal field splitting induced by the defect potential [11]. The spin–orbit coupling creates light-hole, heavy-hole and a split-off hole in the VBM, where the heavy-hole has the shortest Bohr-radius effective mass state orbitals [11]. According to the theory from Thiering and Gali [11], the 1s, 2s, … effective mass states will transform as A1g, whereas 2p, 3p, … effective mass states will transform as A2u and E u , and 3d, 4d, … effective mass states will transform as A1g and E g . As a consequence, only the p-type effective mass states can be optically excited from the A2g ground state. The lowest energy 1s effective mass state may be observed with the contribution of u-type of phonons as phonon-assisted optical transition. Since the translation motion of Si ion in the void of V2 transforms with u odd-parity [130] the 1s, 2s, … as well as 3d, 4d, … effective mass states could be optically excited via the A2u quasi-local phonon mode of the Si ion which is about 43.4 meV according the PBE calculations [11]. Indeed, the 1s effective mass state was not observable in the PLE spectrum but well detectable in the optical spin-polarisation spectrum mediated by the Si-ion vibrations [11] (see Figure 6(a)).

By applying Rydberg scaling to the experimental data (Figure 6(b)), one can find that the expected binding energy of the 1s level is about Δ1 − 0.04 = 0.19 eV deeper when the experimental data is corrected with the phonon energy of the Si ion vibration.

An important observation is that the central cell correction makes the 1s level shallower (i.e., its binding energy becomes smaller) than the value of the effective mass theory, in stark contrast to the case of shallow donor and acceptor dopants in semiconductors. The central cell correction energy can reach hundreds of meV for deep defects in diamond.

The qualitative explanation behind this observation can be drawn from Figure 5: the majority of the 1s effective mass state are localized in the core region of the defect in which the localized orbitals are confined. As a consequence, the electron cloud of the localized orbitals will shield the effective attractive potential of the defect and repel the 1s effective mass state which finally shifts its energy level closer to the ionisation threshold energy. The quantative prediction of the 1s energy level calls for ab initio calculations. Unlike the case of the deep defects in Si with bound exciton excited states, the relatively short Bohr-diameter of the VBM hole 1s state makes it possible to embed the excited state wavefunction in a few-thousand atom supercell, viable at KS DFT level. From a remote distance, the SiV(−) + hole system looks completely neutral which is the case of a giant supercell by completely embeding the 1s wavefunction. In smaller sizes of supercells, the systems looks like a negatively charged defect as the hole wavefunction is completely delocalized within the applied supercells, and a total energy correction should be applied similar to the SiV(−) defect. However, the total energy correction for the SiV(−) + hole system should not be exactly the same with the SiV(−) defect as the size of the supercell is increased because the bound hole provides a screening towards the SiV(−) defect core. Thiering and Gali suggested the following formula [11],

(22)En=1corr=ALexpDL+BL3+C,

where D is the screening length which effectively screens the monopole charge induced by the defect, the s-type spherical potential. The quadropole term B may also incorporate the strain field effects too, and C is the value of the dilute limit. It is a critical issue how large could be the screening length. This can be illustrated by numerical modelling of a hydrogen atom in the simple cubic supercell with lattice constants (L) which can be calculated at Hartree–Fock level within VASP with using a soft PAW potential. The results are shown in Figure 7(d). One can clearly see that the Coulombic scaling (−1/L) deviates at sufficiently large supercells. At sufficiently large supercell size (L > 4 Å), the total energy of the system converges exponentially to a constant energy. The dilute limit (C) is found by

(23)EH(L)=ALexpDHL+C,

where the Coulomb interaction (A/L) is screened by expDHL. We note that the repulsive 1/L3 term is missing because only a single proton appears in the system. Using this fitting procedure, DH=1.90Å=3.56a0 was found, where a0=0.53Å is the Bohr-radius of the isolated free hydrogen atom for n = 1 1s state.

Figure 7: 
Excitation processes of SiV(0). (a) Charge transition level of SiV between neutral and negative charge states by means of HSE06 and PBE functionals. In the dilute L → +∞ limit, the HSE06 results (1.55 eV) agree with the experimental data at 1.53 eV. (b) n = 1 1s bound exciton excitation by means of HSE06 and PBE functionals. Here we can see that the HSE06 limit at L → +∞ with screening included can explain the experimentally data at 1.39 eV. (c) Scaling of the e
g
 ↔ e
u
 excitation process by means of HSE06 functional. (d) Total energy of the hydrogen atom in a Hartree–Fock Γ-point calculation in a simple cubic supercell as a function of the size of the supercell. (e) Schematic of the hydrogen atom in vacuum. The electron is effectively closed into a L3 box. However, it is effectively not a box as it warps around its edges due to the interaction with its periodic images. From a sufficiently large L≳4Å$L\gtrsim 4&#x30a;{\mathrm{A}}$ distance, the H-atom in the supercell can be interpreted as a free non-interacting H-atom. By adding the atomic energy of the employed soft PAW potential for H ion, which causes an artificial constant shift in the absolute total energy.
Figure 7:

Excitation processes of SiV(0). (a) Charge transition level of SiV between neutral and negative charge states by means of HSE06 and PBE functionals. In the dilute L → +∞ limit, the HSE06 results (1.55 eV) agree with the experimental data at 1.53 eV. (b) n = 1 1s bound exciton excitation by means of HSE06 and PBE functionals. Here we can see that the HSE06 limit at L → +∞ with screening included can explain the experimentally data at 1.39 eV. (c) Scaling of the e g e u excitation process by means of HSE06 functional. (d) Total energy of the hydrogen atom in a Hartree–Fock Γ-point calculation in a simple cubic supercell as a function of the size of the supercell. (e) Schematic of the hydrogen atom in vacuum. The electron is effectively closed into a L3 box. However, it is effectively not a box as it warps around its edges due to the interaction with its periodic images. From a sufficiently large L4Å distance, the H-atom in the supercell can be interpreted as a free non-interacting H-atom. By adding the atomic energy of the employed soft PAW potential for H ion, which causes an artificial constant shift in the absolute total energy.

One can conclude from the results of this simple model that the screening length is multiple times longer than the Bohr-radius of the effective mass state. This result could explain the need of simplified computational approaches on silicon defects described above because the Bohr-radius itself is already too long to be accommodated by 10,000 atom supercell, so the screening radius cannot be computed at ab initio level. In those cases, Eq. (21) was applied [20] which overcorrects because of the neglect of the screening effect. In diamond, the Bohr-radius of the effective mass hole states is short enough to observe the deviation from the formula in Eq. (21) due to the screening effect in few thousand atom supercell calculations. Thiering and Gali in Ref. [11] applied PBE functional to yield the 1s total energy by ΔSCF method including 8000-atom supercells, and the resulting screening length was fixed in the fit to HSE06 ΔSCF energies as a function of the supercell size where the maximum size was 1000-atom supercell (Figure 7(b)). The scalings of the charge transition level (Figure 7(a)) and the excitation energy of the 1s state show similarities in the range of small supercell size but a clear deviation can be observed for supercell size with >1000 atoms. That deviation is essential to obtain the accurate value, in good agreement with the experimental data. Certainly, the highly localized orbitals and the corresponding excitation energy converges much faster with the size of the supercell as demonstrated in (Figure 7(c)) which also holds for the crystal field splitting and spin–orbit coupling parameters too [11].

These recent findings cannot be found in the textbooks about semiconductor physics and can be considered as an extension of effective mass theory towards the excited states of deep defects.

The bound exciton states may be formed not just upon optical excitation but by capturing carriers by the deep defects. The capture rate can be, in particular, effective if the defect Coulombically attract the carriers. In the previous example, it can be imagined that the SiV(−) defect in its ground state captures a hole from VBM, which turns SiV(−) to SiV(0) plus a bound hole excited state, and that will decay either radiatively or non-radiatively to the ground state of SiV(0). This effect was first considered in the electroluminescence of single NV defect in diamond (see Supplemental Material in Ref. [99]). In this case, the negatively charged NV defect binds a hole, creating NV(−) plus bound hole from the VBM. That was calculated by HSE06 ΔSCF method without studying the convergence of the excitation energy [99]. In a later study [100], the interaction of two individual NV defects was investigated when the photo-ionisation of one NV leads to the emission of holes toward the neighbour NV(−) defect which can capture that hole. According to the interpretation of the measurements, a giant σ h ≈ 3 × 10−3 μm2 hole capture rate was derived for NV(−). This was rationalised by involving the formation of a bound exciton state featuring an electron localized at the NV(−) plus a bound hole from VBM. Flick in Ref. [100] calculated the total energy of this bound exciton excited state by following the recipe in Ref. [11] with a little modification of Eq. (22) to apply the screening effect also on the quadrupole term. It was found that the binding energy of the exciton is about 40 meV for the 1s state [100]. That is definitely a stable state at room temperature.

3.4 Computation of zero-field splitting for high-spin defects and the g-tensor for defects with heavy ions

The computation of magneto-optic parameters of defect qubits such as zero-field splitting (ZFS) or g-tensor is of high importance not only because that they act as a fingerprint for identification of defect qubits with unknown microscopic structure but they can govern the type of interaction with external magnetic, electric, and strain fields as well as temperature. The importance of interaction between the hyperfine tensor with the strain field was already illustrated in Section 3.3.4. Here, we briefly list the advance in the calculation of ZFS and g-tensor for defect qubits.

The high-spin (S ≥ 1) defects with axial or lower symmetry may experience the electron spin-electron spin dipole-dipole interaction that may be expressed as

(24)HSS=μ0ge2μB24πi>j3ŜirijrijŜjŜiŜjrij2rij5i>jŜiDijŜj,

where r ij = r i r j . The 3 × 3 D-tensor can be diagonalized to find the spectrum and spin eigenstates. The D tensor is associated with the two-particle spin density matrix, n2(r1, r2), which can be approximated by using the Slater-determinant of the KS wave functions ϕ of the considered system, so that n2(r1,r2)Φij(r1,r2)2, where Φij(r1,r2)=12ϕi(r1)ϕj(r2)ϕj(r1)ϕi(r2) and then

(25)Dab=12μ04πge2μB2S(2S1)i>joccupiedχij×Φij(r1,r2)2r2δab3rarbr5d3r1d3r2,

where ra,b=r1r2a,b and χ ij is either 1 or −1 for KS i, j states of the same or different spin channels, respectively. Note that in DFT the spin-polarised KS states are not spin restricted. Consequently, not only the unpaired KS states but also the rest of the occupied states can contribute to the spin density and the ZFS [131]. However, there is no guarantee in spin-polarised KS DFT methods that the final solution will be the eigenstate of the spin operator, and the discrepancy in it is called spin contamination. In practice, spin contamination may be small with using (semi)local DFT functionals but could be significant with hybrid DFT functionals. This may result in a significant error in the approximation of the two-particle spin density matrix. Biktagirov and co-workers suggested a workaround for this problem [132] which is illustrated for the S = 1 case. The idea is that if there is a spin contamination in the electronic structure then the m s = 0 spin configuration with spin-polarised DFT should produce non-zero contribution to the D-tensor, called D̄ms=S1. The bare D-tensor with m s = S spin-polarised DFT results in Dms=S. Finally, the corrected D-tensor, D̃, is

(26)D̃=S2(DD̄).

Each of these m s = S − 1 configurations can be obtained by changing the occupation of one of the half-filled KS orbitals from spin up to spin down and subsequently performing the self-consistent field calculation. It was found that for divacancy defects in 4H SiC that the calculated D-constants are at around 1.6 GHz but D̃=1.3 GHz are obtained after correction, close to the experimental data (see Ref. [132] and references therein).

The ZFS may have other contribution for S ≥ 1 systems than electron spin-electron spin dipole-dipole (DSS) interaction as given in Eq. (25). As an example, we mention here the neutral nickel-vacancy (NiV) defect in diamond which has the same structure as SiV(0) defect discussed above in this review paper. NiV(0) has also 3A2g ground state similar to that of SiV(0) with six carbon dangling bonds which constitute of the ground state electron wavefuntion. The calculated DSS = 967 MHz = 0.004 meV which is typical for the third neighbour distance of dangling bonds in diamond (see Table 2 below). However, this value is very far from the observed D = 170 GHz = 0.703 meV (see Ref. [133] and references therein). As 3A2g is an orbital singlet, thus first-order spin–orbit interaction does not enter here. However, second-order spin–orbit interaction between the ground state triplet and excited state singlet states may play a role which can be selective towards the m s = 0 state of the triplet, and it results in an effective energy shift of the m s = 0 level and opening the gap between the m s = 0 and m s = ±1 levels. To illustrate this using the first-order perturbation theory, we consider the interaction of the 3A2g and 1A1g state that are linked by the parallel component of the spin–orbit operator (ĤSO), λ z , where

(27)ĤSO=iλL̂i,xŜi,x+L̂i,yŜi,y+λzL̂i,zŜi,z.
Table 2:

Spin–orbit and spin-spin contributions to the ZFS (in MHz) of the NV centre and a set of group-IV–vacancy defects in diamond (compared to available experimental data). All the defects have S = 1 ground state with orbital singlet A2-type many-body wavefunction. The geometry and electronic structure of group-IV–vacancy defects are akin to those of SiV(0). See Ref. [132] and references therein.

Defect D SO D SS D SO+SS Experiment
NV(−) 6 2722 2728 2878
SiV(0) 480 570 1050 929
GeV(0) 1469 630 2099 2248
SnV(0) 10,763 630 11,393
PbV(0) 144,860 660 145,520

The energy gap between 3A2g and 1A1g is Θ before applying ĤSO and the spin–orbit coupling between |3A2g, m s = 0⟩ and 1A1g through λ z is λ0. In that case, the first-order perturbation theory yields that |3A2g, m s = 0⟩ level shifts downwards by λ02/Θ. This means that

(28)DSO=λ02Θ

and it will be dominant over DSS. According to HSE06 calculations (see details in Ref. [133]), Θ ≈ 0.68 eV and λ0 = 23.2 meV which results in DSO = 0.79 meV. This is much closer to the experimental data at 0.703 meV.

One can go beyond the first-order perturbation theory and consider the change in the wavefunction due to spin–orbit interaction (second-order perturbation),

(29)Ψ(1)=|3A2g,ms=0)+λ0Θ|1A1g,

which may result in a more accurate result than that by first-order perturbation theory.

In general, the problem can be rephrased by considering the total energy of the system as a function of the spin quantisation direction, Etot(nσ). In a uniaxialrt: case, the magnetic anisotropy energy is then defined as difference between Etot(nσ) obtained with nσ) parallel (z) and perpendicular (⊥) to the anisotropy direction, ESO = Etot(z) − Etot(⊥), and then the corresponding D-constant can be evaluated as DSO = ESO/S2 for integer S and DSO = ESO/(S2 − 1/4) for half-integer S. For the case of NiV(0) with S = 1, ESO should be calculated self-consistently. Etot(z) = Etot(↑↑) calculation can be carried out with the usual m s = +1 setting for the KS orbitals together with the scalar-relativistic spin–orbit interaction. Etot(⊥) calculation is a bit tricky. The two spins should be individually rotated by 90° about the x-axis. Then the total energy of the system should be calculated by scalar-relativistic spin–orbit interaction. This total energy is not identical to Etot(⊥) because the rotated-spins system is not the exact m s = 0 eigenstate but a mixture of m s = ±1 and m s = 0, so we label it by Etot(→ →). The final expression [133] is

(30)DSO=3Etot()Etot().

The self-consistent HSE06 DSO = 0.73 meV which is 0.06 meV deeper than the first-order perturbation theory value at 0.79 meV, and it brings the result closer to the experimental data at 0.703 meV. This shows that self-consistent spin–orbit calculation needed for obtaining accurate ZFS for defects consisting of heavy ions. It is interesting to note that self-consistent spin–orbit PBE calculations results in DSO = 1.35 meV which is significantly larger than the HSE06 and experimental values. We note that Θ ≈ 0.25 eV with PBE which explains the too large DSO with PBE as first-order perturbation theory showed that DSO scales inversely between the gap of the triplet and singlet levels [133].

These results clearly demonstrate [133] that the energy gap between triplet (high-spin) and singlet (low-spin) levels are highly critical in obtaining an accurate ZFS for defects which consist of heavy ions.

The disasvantage of the self-consistent spin–orbit calculations is that it can principally work for sufficiently large spin–orbit energies, usually created by heavy atoms, so that it does not fall below the numerical noise. For defects with light atoms, one has to rely on the first-order perturbation theory which was previously sketched for a special case (diamond NiV defect) as an introduction to the problem.

Biktagirov and co-workers [134] implemented the perturbation theory based method to the GIPAW-tree of the Quantum Espresso package. They apply collinear spin polarisation approximation with direction a = x, y, z. Then the SO coupling in direction a ĤaSO and b ĤbSO contributes to the total energy of the system in second-order perturbation theory as

(31)EabSO=o,s,sReΨ̄osĤaSOGs(ϵe)ĤbSOΨ̄os.

Here the sum runs over the spin channels s and s′ and the occupied states os. Thereby, |Ψ̄os are the corresponding unperturbed KS wave functions (obtained without SO coupling), and G s (ϵ e ) is the Green’s function of the empty states es′. In their implementation, explicit summation over empty states is avoided by calculating Gs(ϵ e ) by projecting the empty states onto the valence bands. By this, the approach becomes faster, numerically more stable, and almost unaffected by the gap issues quoted above for the (semi)local DFT functionals. On the other hand, the method is not so intuitive as the majority of the interaction comes from the closest high-spin low-spin energy states which cannot be directly analysed by this method.

Biktagirov and co-workers applied their method to diamond NV centre and group-IV–vacancy defects [134]. The results are listed in Table 2. As can be seen the defect with the lightest atom exhibits the smallest DSO whereas it increases orders of magnitude with going to heavier atoms. In the group-IV–vacancy defects DSS increases slightly as the heavier atoms push the neighbour carbon atoms farther from each other but the vast contribution comes from DSO except for SiV(0) and partially for GV(0) where the two contributions are similar. We note that PbV(0) shows about DSO = 145 GHz for which self-consistent DSO calculation would result in a lower value.

Although, the calculated DSO is only 6 MHz for diamond NV centre but it couples directly to the electric field unlike DSS which couples to the electric field only indirectly via changing the electron cloud so the spin density (e.g., Ref. [135]). As a consequence, the field applied along the defect’s symmetry axis, the DSO part dominantly drives the predicted Stark coefficient, 0.034 GHz Å/V, into the experimentally observed confidence interval of 0.035 ± 0.002 GHz Å/V (see Ref. [134]). The simulation was carried in a (111) diamond slab where the electric field was switched on during the calculations of DSS and DSO.

Previously, we discussed the spin level structure in the absense of external magnetic field. Nevertheless, it is highly important to understand the coupling of defect spins to external magnetic fields. The external magnetic fields could be intentionally turned on for manipulation of the qubits on one hand, and on the other hand, randomly distributed electron or nuclear spins proximate to the defect qubits could influence their longitudonal relaxation and coherence times. Here, we discuss the issue of a constant macroscopic external magnetic field interacting with the defect’s electron spin which can be generally described as

(32)Ĥ=Bμ̂=BμBgS̃,

where μ̂ is the magnetic dipole momentum operator, μB is the Bohr-magneton of the electron, g is the g-factor, and S̃ is a phenomenological pseudo-spin, which is set to the net electron spin of the system, i.e., S̃=1/2 for the Kramers-doublet defect’s electron spin. Eq. (32) has the form of the Zeeman formula for the free electron but g e = 2.0023 free electron scalar value is substituted by g. Unlike the free electron case, the defect’s electron spin feels the potential of ions which is less symmetric than spherical, the defect’s electron may have an effective angular momentum with this condition, e.g., localized on the d-orbital, which can be also influenced by the electron-phonon coupling. All of these effects are packed into a single tensor g.

This issue is illustrated on the neutral vanadium defect substituting the Si-site in 4H SiC which has become a very promising spin-to-photon interface with a quantum memory and optical emission at the telecom wavelength (see Refs. [136, 137] and references therein). The d-orbital of the vanadium ion splits due to the C3v symmetric crystal field of 4H SiC and then a double degenerate e-orbital occurs in the gap localized on the d-orbital of vanadium, occupied by a single electron. Because of the double degenerate d-orbital, one can expect an effective spin–orbit coupling between the orbital and the electron spin, where the low symmetry crystal field will reduce the effective angular momentum of the orbital called Stevens reduction factor (r) (e.g., see the origin of this effect in more detail for group-IV–vacancy defects in Ref. [138]). However, it is also known that this is an Ee Jahn–Teller system [139, 140] which can also effectively reduce the angular momentum of the electron orbit so the effective spin–orbit splitting known as Ham reduction factor (p). As a consequence, μ̂ in Eq. (32) can be written as

(33)μ̂=(μBprL̂+μBgeŜ)=μBgS̃,

where the contributions of L̂ and Ŝ are separated, so it gives an opportunity to unravel the microscopic origin of g. Csóré and Gali carried out ab initio calculations [141] to determine the r and p factors, so then g can be obtained as

(34)g=2geSz+Lzeff=2μzμB,
(35)g=μ++μ+i(μμ+)μB,

where expectation values of the ladder magnetic dipole moment operators are used (μ±) to express g. In Eq. (34) S z and Lzeff are expectation values of Ŝz and the effective angular momentum operator, L̂zeff=prL̂z, respectively.

We note that the d-orbitals may require a special attention for accurate calculation, and indeed, HSE06 DFT functional overlocalises the d state that should be corrected [142]. For the heavy-atom defects one may assume that the vast majority of the spin–orbit coupling, so the effective angular momentum, comes from the single heavy-atom. In this case, the analysis of the d-orbitals and the actual DFT wavefunctions can reveal the deviation of the d-orbitals from the spherical symmetry (see Table 3), so the r can be computed. As an example, Ψ1 and Ψ4 as well as Ψ2 and Ψ3 can be coupled in Eq. (34) where the corresponding wavefunction coefficients can be extracted for the given spin–orbit state in the KS DFT calculation. It was found that the vanadium at one site of 4H SiC feels isotropic environment with small effective angular momentum which is an order of magnitude larger in a truly axial-like environment at the other site of 4H SiC. After solving the Jahn–Teller Hamiltonian (see Refs. [30, 141]), the typical Ham reduction factor is at about 0.6 which is significant, so the electron-phonon coupling effectively reduce the interaction between the defect’s pseudo-spin and the external magnetic field parallel to the symmetry axis of the defect. The final typical values of g are around 1.9.

Table 3:

Kramers doublets formed by d-orbitals and the corresponding single and double group irreducible representations under C3v symmetry. We give widespread notations for double group irreducible representations (irreps) and also the corresponding m j = m l + m s values.

Labels Orbitals Irreps. m j
Single Double
Ψ1 |d+2,+12;|d2,12 2 E E124) ±52
Ψ2 |d+2,12;|d2,+12 2 E E325,6) ±32
Ψ3 |d+1,+12;|d1,12 2 E E325,6) ±32
Ψ4 |d+1,12;|d1,+12 2 E E124) ±12
Ψ5 |d0,+12;|d0,12 2 A 1 E124) ±12

For the calculation of g (Eq. (35)), the ladder magnetic dipole operator was considered, μ̂± that can couple state |m j ⟩ to state |m j ± 1⟩, where m j = m l + m s . However, the in-gap defect states transform as either E1/2 (linear combination of Ψ1 and Ψ4) or E3/2 (linear combination of Ψ2 and Ψ3) with the m j values given in Table 3. Consequently, μ̂± cannot couple neither Ψ1 and Ψ4, nor Ψ2 and Ψ3 therefore g = 0 in each case. Deviation from g = 0 may occur due to secondary effects.

The same method was applied to the nickel defects in diamond where the NiV(−) was identified by first principles calculations as an excellent qubit candidate analogous to the group-IV-vacancy qubits in diamond [133] which has an optical emission at about 1.4 eV and a highly anisotropic g-tensor. In the literature, the NE4/AB1 EPR centre with S = 1/2 spin and relatively isotropic g-tensor with g = 2.0027(2) and g = 2.0923(2) was previously associated with NiV(−) which is linked to the 1.72-eV optical centre (see Ref. [143] and references therein). Clearly, the NE4/AB1 centre should be associated with another nickel-related defect in diamond. Thiering and Gali tentatively assigned Nis(Ns)3(0) defect to this centre which has an unpaired electron on the a1 orbital strongly localized on Ni 3d orbitals. In this case, the g-tensor is modified from the free electron value because of the orbital moment of the Ni 3d states as explained above for vanadium defects in 4H SiC. This justifies to calculate the total orbital moment (L̂x,y,z) of the defect within the PAW sphere of the ions where the largest contribution comes from the Ni ion. Because of C3v symmetry L̂x=L̂y=L̂. The main components of the g-tensor can be given as g=ge+2L̂z and g=ge+2L̂. Finally, g = 2.0058 and g = 2.0942 are obatined, in good agreement with the experimental data.

3.5 Spin-phonon coupling: temperature-dependent longitudonal spin relaxation time and magneto-optical parameters

3.5.1 Longitudonal spin relaxation time

A key parameter of qubits is the longitudonal spin relaxation time which is the characteristic time of flipping the spin, and it is labeled by T1 in the literature. This sets the absolute limit of the spin coherence time, i.e., the characteristic quantum information processing operation time of the qubit. It is of high importance to understand the underlying microscopic processes. In nuclear spin physics, the origin of spin flipping was identified as the interaction between phonons and the spin which is manifested as a highly temperature-dependent phenomenon; therefore, it is also often called spin-lattice relaxation time. As T1 often exponentially decay with elevating the temperature it is imperative to characterise T1 as a function of temperature, in order to explore the applicability of qubits as sensor probes in biology which requires ambient conditions. Our review paper only focusses to the recent advances on defect qubits, in particular, on S ≥ 1 defect qubits.

Spin flipping processes require such interaction Hamilton operator which contains spin shift operators. It can be easily recognized that the spin-spin dipole-dipole interaction in Eq. (24) contains single and double spin shift operators, e.g., ŜxŜz and ŜxŜx, respectively. Therefore, if the defect qubit’s spin interact with other spin species then it causes a spin flip of the defect qubit. The strength of dipole-dipole interaction goes with inverse cube of the distance between the spins. This interaction is weakly dependent from temperature and it highlights that the longitudonal spin relaxation time is not necessarily a spin-phonon interaction. In practice, this type of T1 process becomes only important at elevated temperatures for diamond NV centre if the concentration of defect spins is relatively high, e.g., 4–8 particle per million (ppm) in diamond [62].

Ivády applied the cluster-correlation expansion (CCE) [144, 145] to model the interaction of the central diamond NV centre’s electron spin with other electron spins such as the environmental NV centres, nitrogen donor spins (labelled as P1 EPR centre), and the 13C nuclear spins also as a function of the external magnetic field and strain [62]. The CCE approach will be shortly discussed in the next chapter.

We note that another study only restricted this investigation to the bath of 13C nuclear spins but taking only the dipolar interactions into account with far 13C nuclear sites [146]. However, the bath of 13C with Fermi-contact hyperfine terms cannot be ignored for accurate simulations which calls for ab initio simulations.

Ivády calculated the hyperfine tensors for 13C isotopes by HSE06 DFT method in a 1728-atom simple cubic supercell [62]. Since the hyperfine tensors should be determined at large distances from the defect site this required a special approach in order to avoid finite size effect problems. Ivády utilised a real space grid combined with the PAW method to calculate hyperfine tensors. The Fermi contact term, dipole-dipole interaction within the PAW sphere, and core polarisation corrections are calculated within the PAW formalism from the convergent spin density. The dipolar hyperfine contribution from spin density localized outside the PAW sphere is calculated by using a uniform real space grid. This procedure enabled to obtain hyperfine coupling tensors excluding effects from periodic replicas of the spin density due to the periodic boundary condition. Additionally, hyperfine tensors for atomic sites outside the boundaries of the supercell were calculated by neglecting Fermi contact interactions with achieving a smooth transition in the hyperfine constants at the boundary of the two approaches.

Ivády used the extended Lindbladian equation in order to simulate the spin dynamics of the central spin and its relaxation rate 1/T1 in materials where the electron spin density cannot be ignored beside the nuclear spin bath [62]. In this model the total spin relaxation rate (1/T1tot) is

(36)1T1tot=1T1P1+1T1NVbasal+1T1NVpar+1T113C,

where “NV-basal” and “NV-par” label such NV centres in the environment which have other and parallel symmetry axis with that of the central NV centre. Finally, it was found that the environmental NV centres have a dominant effect on the spin relaxation rate [62]. At special setting of the magnetic fields, either ground state level anticrossing (GSLAC) or excited state level anticrossing (ESLAC), the relaxation rate is accelerated because the P1 centres and the nuclear spins can easily induce spin flip-flop processes that were otherwise protected by the energy gap between the electron spin levels of the NV centre. At GSLAC, the spin-polarisation of the electron spin and the coupled nuclear spin also changes that can be observed by the change of the PL intensity as the external magnetic field is swept around the GSLAC region [147]. This modelling also rationalised the photo-electric read-out process of the single 14N nuclear spin of the NV centre at ESLAC condition [148].

T1P1 time was further investigated in detail [149] in which they considered the microscopic structure of the P1 centres in diamond, namely, the strong Jahn-Teller distortion will generate four different symmetry ⟨111⟩ axes in diamond, and the 14N nitrogen hyperfine tensor to the P1 centre’s electron spin adapts to these orientations. As a consequence, the spin flip-flop processes between the P1 pairs are reduced with respect to the case of unrealistic aligned P1 centres [149]. With taking the microscopic structure of the P1 centres into account, the calculated spin relaxation times of the diamond NV centre exhibits a clear linear dependence on P1 concentrations on a log scale with a slope of −1.06, in excellent agreement with some experimental data (see Ref. [149] and references therein).

This theory was also applied to the divacancy qubit in 4H SiC by considering other divacancy spins (S = 1), negatively charged Si-vacancy spins (S = 3/2), nitrogen donor spins (S = 1/2), as well as 13C and 29Si I = 1/2 nuclear spins in the environment, also as a function of the external magnetic field [150]. It was found that the cross-relaxation accelerate spin flip-flop rates again in the region of GLSAC and ESLAC magnetic fields for each considered environmental spins. At zero magnetic field a simple relation was found for the interaction between N-donor and the central divacancy spin,

(37)1T1=βC2,

where β = 1.6 × 10−35 Hz/cm−6 and C is the concentration of the N-donor. It is noted that nitrogen implanted samples the distribution of nitrogen donor is not homogeneous, and then the “concentration” should be considered near the target divacancy spin which was created as a result of the implantation [150]. The theory was also employed to the Si-vacancy S = 3/2 qubit in 4H SiC [151]. In this case, Bulancea-Lindvall and co-workers considered the interaction the Si-vacancy qubit spin with S = 1/2 defects, e.g., N-donors. Si-vacancy in 4H SiC has minor ZFS, thus at a given external magnetic field with similar Zeeman shifts the two spin systems can be effectively coupled by dipole-dipole interaction unlike the case of divacancy with S = 1 spin and high ZFS (1.3 GHz). On the other hand, 29Si nuclear spins at natural abundance (4.5%) exhibit a considerable hyperfine splitting when they interact with the Si-vacancy qubit spin which will supress the cross-relaxation between the Si-vacancy qubit spin and the S = 1/2 spins in the environment. By isotope purification (reducing the concentration of 29Si isotopes in 4H SiC), the cross-relaxation so the spin flip-flop process accelerates and the spin relaxation time of the Si-vacancy qubit significantly reduces when the concentration of S = 1/2 spins are high (10151018 cm−3). This surprising result was unraveled by these simulations with using ab initio hyperfine tensors in the parametrisation of the interaction Hamiltonian.

High quality materials with single defect spins and low concentration of nuclear spins do not experience spin flip-flop events with electron spins in the environment, and the flip-flop processes caused by nuclear spins are only observable at special conditions (e.g., external magnetic field is set to GSLAC condition). In this case, the spin-phonon coupling is responsible for T1 and it becomes strongly temperature dependent. In the case of molecules, vibrations are indeed responsible behind the spin flipping process. In a recent review paper [152], the ab initio theory and its application to molecules have been presented in detail. The formulas and basic equations apply to defect qubits too which are not reiterated here in detail.

Regarding the temperature dependence of T1 = 1/Γ of defect qubits, the most studied one is the diamond NV centre [31, 153], [154], [155], and the first ab initio results have been reported for this defect qubit because the theories could be well tested on the accurately recorded experimental data. In recent studies [31, 156], the m s = ±1 levels of 3A2 ground state was split by a small external magnetic field aligned parallel to the symmetry axis of the defect (145 MHz), and they could observe both single-flip rates (|0⟩ ↔ |±1⟩, labelled by Ω) and double-flip rates (|+1⟩ ↔ |−1⟩, labelled by γ) and found that γ > Ω at any observed temperature (T > 125 K) and γ ≈ 2 Ω at room temperature. In the temperature region of 125 K and 400 K, the rates increase from 1 Hz to 200 Hz. At this temperature region only phenomenological theory was considered using an empirical model in which the high-temperature behavior is characterized by a term that scales with temperature as T5 [153, 154], which may arise due to Raman scattering of low-energy acoustic phonons which are weakly coupled to the spin via first-order interactions. However, insights from ab initio simulations should verify this model including the magnitude of the double-flip rates.

The spin relaxation rate may be expressed as

(38)Γ=Γ0+Γ1(1)(T)+Γ1(2)(T)+Γ2(1)(T)+,

where the superscript refers to the order of the spin-phonon interaction (i.e., terms with superscript 1 or 2 are linear or quadratic in the atom displacements respectively) and the subscript refers to the order in perturbation theory. Γ0 is a sample-dependent constant term arising from spin-spin interactions that was discussed above. Γ1(1) describes the absorption or emission of a single resonant phonon by the spin. Because the ZFS energy of the NV ground-state triplet is small in comparison to typical phonon energies in diamond, this process is only relevant at subkelvin temperatures [155]. Γ2(1) corresponds to the Raman-scattering of low-energy acoustic phonons via first-order interactions. However, we will present below [156] unlike for other spin systems (e.g., several coordination compounds as recently shown in Ref. [157]), first-order spin-phonon interactions provide only negligible contributions to Raman scattering for the NV centre in diamond. The major effect comes from Γ1(2)(T) so the quadratic displacements of ions.

In order to calculate Γ1(2)(T), the spin-phonon matrix elements should be obtained by exploiting the dependence of the D-tensor on the normal coordinates (Q) as

(39)D(R)=D(R=0)+iDQiR=0Qi+12ijDQiQjR=0QiQj,

where the coefficients in Eq. (39) were extracted from VASP PBE calculations as implemented by Thiering and Gali [156]. In order to evaluate the second-order derivatives, only the diagonal terms were considered which satisfy i = j and distort the C3v symmetric atomic positions by all degenerate e x , e y phonon modes of the supercell by (ΔR)2=0.1Åa.m.u.. The second-order spin-flipping matrix elements V+0ll and V+ll then determine the D-tensor according to the symmetry-adapted expression in which only the quadratic terms are expresse here as

(40)1313+23i3V00iiRi2+lV+ll×11Xl2Yl2+112XlYl+l2V+0ll11Xl2Yl2+112XlYl

where R i , X l and Y l are the dimensionless coordinates (not normal coordinates) for the phonon modes at energy ℏω i or ℏω l . We note that while the index l only covers the e modes once, the index i covers all a1, a2, e x , e y modes and thus runs over the e modes twice. Eq. (40) can be employed to transform it into the spin-phonon interaction V̂. By expicitly writing only the quadratic terms, it reads as

(41)V̂=SDS=DŜz213S(S+1)+i3V00iiŜz213Ŝ(Ŝ+1)R̂i2+lV+llŜx2Ŝy2X̂l2Ŷl2+ŜxŜy+ŜyŜxX̂lŶl+ŶlX̂l+l2V+0llŜxŜz+ŜzŜxX̂l2Ŷl2+ŜyŜz+ŜzŜyX̂lŶl+ŶlX̂l,

where S=(S)=ŜxŜyŜz. The dimensionless coordinates are expanded in terms of the phonon creation and annihilation operators: R̂i=bi+bi/2 and {X̂,Ŷ}l=b{X,Y}l+b{X,Y}l/2.

As a next step for defining the equations for the rates, one can apply RPA for these processes so it is assumed that the consequtively absorbed phonons are not coherent. Furthermore, one can further simplify the equations by considering the fact that the ZFS energy of the diamond NV centre is much smaller than the typical phonon energies coupled to the spin. The key matrix elements are the first-order spin-phonon coupling coefficients Vmsmsl from the first-order spin-phonon interaction V̂(1)=lmsmsVmsmslal+al and the second-order spin-phonon coupling coefficients V̂(2)=llmsmsVmsmsllal+al(al+al) from Eq. (41). The equations involve Dirac δ functions for conserving the energy in the process. However, finite size DFT supercell calculations do not produce continuous phonon density of states, so it is required to use a Gaussian convolution for Vmsmsl and Vmsmsll which results in first-order Fmsms(1)(ω)=l|Vmsmsl|2δ(ωωl) and second-order Fmsms(2)(ω,ω)=ll|Vmsmsll|2δ(ωωl)δωωl spectral functions in the continuum limit, respectively [156].

Finally, the appropriate spin relaxation rates are

(42)Γ2msms(1)(T)=4πms0d(ω)nB(ω)[nB(ω)+1]Fmsms(1)(ω)Fmsms(1)(ω)(ω)2

and

(43)Γ1msms(2)(T)=4π00d(ω)d(ω)nB(ω)nB(ω)+1Fmsms(2)(ω,ω)δ(ωω),

respectively. The temperature dependence enter via the Bose–Einstein occupation function (nB) of the phonons at ω, ω′ energies. In the current implementation [156], ℏω′ = ℏω constraint was employed in Eq. (43) that also enforces l = l′ and so Vmsmsll diagonal matrix elements were considered in the reported ab initio calculations [156].

The numerical ab initio calculations provided very slow rates for Γ2msms(1). The basic reason behind this observation is that ℏω)2 in the denominator completely suppresses Fmsms(1) because the largest values of Fmsms(1) are typical in the range of 1–10 MHz whereas the relevant phonon frequencies are in the order of 10 THz. This can be rationalised by noting that the first-order interaction strength is roughly DΔu/a, where D is the zero field splitting (2.8 GHz), Δu is the atomic displacement, and a is the nearest neighbor distance in diamond. In contrast, contributions to Raman scattering from second-order interactions scale quadratically with the second-order interaction strength, approximately Du/a)2. Thus the ratio between the first-order and second-order contributions is on the order of (D/ℏω)2 ∼ 10−8, indicating that Raman scattering due to first-order interactions can be neglected for the diamond NV centre and other S = 1 defects in diamond [31].

The calculated second-order spin-phonon coupling coefficients are depicted in Figure 8. Two broad peaks can be observed at certain phonon frequencies that are associated with the motion of the carbon dangling bonds so the spin density (see Ref. [156] and references therein). As the frequencies of the effective phonons are much higher than the thermal energy of the measurement temperature this will scale as Orbach-process. Ab initio simulation revealed that two effective phonon frequencies exist, thus Γ = 1/T1 can be described as a double Orbach-process, where the higher effective frequency plays a role at elevated temperatures [156]. The theory also well describes the double-flip transition and the appropriate rate equations, and both processes are double Orbach-processes. The double Orbach-process parameters could be well fitted to the experimental data with providing 68.2 meV and 167 meV effective phonon frequencies (gray lines in Figure 9) which agree well with the features of the calculated spin-phonon spectral functions. Insights from theory provided a physically well motivated model for the spin-lattice relaxation times of diamond NV centre. The calculated rates are depicted in Figure 9. The agreement between theory and experiment for γ is very good whereas a larger discrepancy is observed for Ω. It was hypothesized that the discrepancy in the predicted Ω is due to the exclusion of combinations of modes for which ll′, as combinations of modes with different symmetries likely account for significant matrix elements associated with pairs of different spin operators, which correspond to the single-quantum transitions. It was also discussed that F00(2)(ω,ω) in Figure 8 plays an important role in the phonon-assisted decoherence process, also showing a double Orbach-process character, which has not yet been recognized in previous works [156]. This will be discussed in the next chapter.

Figure 8: 
Ab initio calculation of the second-order spin-phonon coupling coefficients (thin lines) and spectral function (thick lines) for a single diamond NV centre in a 512-atom supercell. NV spin-phonon dynamics are characterised by the magnitudes of the matrix elements ŜzŜ+${\hat{S}}_{z}{\hat{S}}_{+}$ (blue), Ŝ+2${\hat{S}}_{+}^{2}$ (red), and Ŝz2−13Ŝ2${\hat{S}}_{z}^{2}-\frac{1}{3}{\hat{S}}^{2}$ (black), which cause single-flip relaxation, double-flip relaxation, and dephasing respectively. The spectral function display peaks near the values of 68.2(17) and 167(12) meV extracted from the fit of the two-phonon Orbach-process model to the experimental data (gray lines and ± 1σ intervals, see text).
Figure 8:

Ab initio calculation of the second-order spin-phonon coupling coefficients (thin lines) and spectral function (thick lines) for a single diamond NV centre in a 512-atom supercell. NV spin-phonon dynamics are characterised by the magnitudes of the matrix elements ŜzŜ+ (blue), Ŝ+2 (red), and Ŝz213Ŝ2 (black), which cause single-flip relaxation, double-flip relaxation, and dephasing respectively. The spectral function display peaks near the values of 68.2(17) and 167(12) meV extracted from the fit of the two-phonon Orbach-process model to the experimental data (gray lines and ± 1σ intervals, see text).

Figure 9: 
Comparison between ab initio theory and experiment. (a) Dotted lines show relaxation rates obtained by evaluating Eq. (43) with the ab initio second-order spin-phonon spectral function shown in Figure 8. Dashed lines show fit of the analytical model (see text) to the experimental data with sample-dependent constants set to zero. (b) Ratio of the ab initio relaxation rates to the analytical model rates. In the phonon-limited regime (gray line) the ab initio theory underestimates the experimentally measured relaxation rates by approximately 16% for γ and a factor of 8 for Ω at room temperature.
Figure 9:

Comparison between ab initio theory and experiment. (a) Dotted lines show relaxation rates obtained by evaluating Eq. (43) with the ab initio second-order spin-phonon spectral function shown in Figure 8. Dashed lines show fit of the analytical model (see text) to the experimental data with sample-dependent constants set to zero. (b) Ratio of the ab initio relaxation rates to the analytical model rates. In the phonon-limited regime (gray line) the ab initio theory underestimates the experimentally measured relaxation rates by approximately 16% for γ and a factor of 8 for Ω at room temperature.

3.5.2 Temperature shifts of magneto-optical parameters

Understanding the temperature shifts of magneto-optical parameters of defect qubits is of high importance in various aspects. One of the most obvious issues is the temperature sensing with defect qubits at the nanoscale which requires temperature characterisation of the basic magnetic parameters. Again, the diamond NV centre is the most investigated defect qubits in this regard (see Ref. [30] and references therein). As an example, the temperature dependence of the ZFS of the diamond NV centre was modelled by the thermal expansion [158] which results in an increase in the distance between carbon dangling bonds so the decrease in the ZFS parameter (D-constant). However, the obtained coupling coefficient was much lower than the experimental data. Recently, Tang and co-workers pointed out [159] that the thermal expansion model covers a “third order” effect as the measured magneto-optical entity ν will be a statistical average of the phonon mode distribution as

(44)ν=ν0+iνQi0Qi+122νQi20Qi2,

where {Q i } are the normal coordinates of the phonons wihtin quasi-harmonic approximation. In the Born-Oppenheimer approximation, the global energy minimum results in ⟨Q i ⟩ = 0 (forces are zero) and it becomes non-zero because of violating of the harmonic approximation, i.e., anharmonicity of the phonons. This can be taken into account in the thermal expansion of the lattice. However, the second term is then expected to be dominating. As Qi2=Miωi(nB(ωi)+12) with nB Bose–Einstein occupation function Eq. (44) can be expressed as

(45)ν=ν0(a)+i122νQi2MiωinB(ωi)+12,

where M i is the mode-specific effective mass and ν0 is calculated at the lattice constant a which corresponds to the thermal expansion at the given temperature. The spectral function was defined as

(46)Si(ω)=j1Mj2νQj2δ(ωωi).

The second derivative was calculated numerically as implemented in VASP in a 128-atom face centred cubic diamond supercell with 3 × 3 Monkhorst–Pack k-point sampling and DFT PBE functional [159]. The theory was applied the D-constant, the hyperfine constant A zz and the quadrupole moment Q zz of the 14N, and the ZPL energy of the diamond NV centre [159]. It was found that temperature dependence of the D-constant could be well reproduced where the dynamical effects play a major role, although, the termal expansion effect cannot be neglected. The ZPL energy shifts were well reproduced too by this theory. However, it is unexpected that the calculated spectral function for D-constant does not show up a peak at around 70 meV phonon energy in the study of Tang and co-workers [159], which is quite visible in the F00(2) spectral function in Figure 8 from Thiering and Gali in Ref. [156]. The temperature shift of Q zz is also well-reproduced (see Ref. [159] and references therein) and good agreement can be anticipated for A zz in comparison to the few experimental data points available to date. All-in-all, this theory seems to be highly promising with good predictive power after achieving the convergent parameters in the ab initio simulations.

3.6 Ab initio theory of coherence of defect spins in solids

The coherence time of the defect qubits’ spin corresponds to a decoherence of the transverse nuclear spin magnetisation which is generally labelled by T2. It is also called spin-spin relaxation time as usually the interaction of the defect qubit’s electron spin with the nuclear spin bath limits its value in high quality (small electron spin bath) materials. These nuclear spins do not precess with the same frequency in real materials which can lead to a distribution of resonance frequencies around the ideal. Over time, this distribution can lead to a dispersion of the tight distribution of magnetic spin vectors, and loss of signal. This is called the dephasing time, labelled as T2*, which is associated with static or slowly varying inhomogeneities in a spin system. T2* is the characteristic decay time of a free-induction-decay (FID) measurement, wherein a series of Ramsey sequencies (π/2 pulse – τπ/2 pulse) are performed with varying free-precession interval τ, and an exponential decay is observed. Dephasing from fields that are static over the measurement duration can be reversed by application of a π pulse halfway through the free-precession interval invented by Erwin L. Hahn [160]. In this protocol, the π pulse alters the direction of spin precession, such that the phase accumulated due to static fields during the second half of the sequence cancels the phase from the first half. The spin phase is then refocused which appears as an echo signal in the spin resonance spectrum, i.e., Hahn-echo signal. The decay of this echo signal, due to magnetic fields that fluctuate over the course of the measurement sequence, is characterized by the coherence time T2. In practice, T2 exceeds T2* by orders of magnitude. We further note that the phonon induced spin relaxation, longitudonal relaxation time T1, as an incoherent process places hard limits on the maximum achievable coherence times, where the theoretical limit is T2 = 2T1.

The calculation of the hard limit of T2 at a given temperature and electon/nuclear spin bath requires the calculation of T1 as described in the previous chapter. For S = 1 defect qubit’s spin, the |±1⟩ spin states may split due to the low-symmetry of the defect or an external small constant magnetic field. In recent studies [31, 156] it has been found for the diamond NV centre that the double-flip transition (γ rate) is even faster than single-flip transition (Ω rate) induced by phonons. However, γ was neglected in previous studies (see Ref. [156] and references therein), thus T2 have likely been overestimates which should be rewritten as

(47)T2,max(SF)=23Ω+γ

for a superposition in the {|0⟩, |±1⟩} single-flip subspace and

(48)T2,max(DF)=1Ω+γ.

for a superposition in the {|−1⟩, |+1⟩} double-flip subspace [161, 162].

If spin-lattice relaxation does not interfere then the central spin (qubit’s electron spin) and the (nuclear) spin bath interaction and their dynamics should be simulated for obtaining the spin dephasing and decoherence times where the simulations should consider how the qubit’s spin is controlled and driven for yielding T2* and T2 times. This is a highly complex process because the qubit’s control starts with the initialisation of the qubit, i.e., spin-polarisation of the qubit’s electron spin and the spin-polarisation may be transferred to the nuclear spin bath depending on the qubit control protocol. At certain magnetic fields, some nuclear spin levels could drive into resonance to the electron spin levels, e.g., GSLAC condition, which can significantly change the dynamics of the defect qubit’s electron spin coupled to the spin bath.

For defect qubits, the CCE approach has been successfully applied that was originally invented by Yang and Liu [144, 145]. It was originally applied to calculating the pure dephasing of the diamond NV centre’s electron spin in the large detuning regime. However, the central spin flip must be considered when the energy relaxation of the diamond NV centre is involved in the nearly-resonant regime, i.e., GSLAC condition [146]. Often, this is called generalised CCE or gCCE approach. To briefly sketch the problem and the neccessity of approximations, an open system S that consists of a central spin s0 and a number of bath spins s i , where i = {1 … n}. The master equation of the open system S to obtain the density matrix ρ can be written as

(49)ρ̇S=i[H0,ρS]+E(ρS),

where the Hamiltonian H0 can be written as

(50)H0=h0+i=1n(hi+h0i),

where h0 is the Hamiltonian of the central spin, h i is the Hamiltonian of the coupled spin s i , and h0i describes the coupling of the central spin and the bath spin s i . The last term on the right-hand side of Eq. (49) accounts for environmental effects that are not included in S, through the Lindbladian E that we already discussed above. One can define h0 and h0i such a way which include the driving fields or other external fields like a constant magnetic field or strain field.

The size of the problem, i.e., the dimension of the Hilbert space, increases exponentially with n, which makes an exact solution unfeasible for large n. To model the dynamics of S it is divided into a cluster CN of overlapping cluster systems, where N is the order of the cluster approximation as illustrated in Figure 10. All of the subsystems are artificially coupled then together through a modified Lindbladian superoperator

(51)L(bil)=lḃilTrCilCilρci×CilρciCil12ρciCilCil+CilCilρci

added to the master equation of each subsystem, where and C0l and C il are Lindblad operators. We consider C0l and C il operators that describe solely spin flip-flop transitions of the central spin. Here b il are time-dependent rates determined from the flip-flops occurring in subsystem other than i. The Lindbladian formalism ensures that all the spin flip-flops occurring in the different subsystems is carried out in all subsystems. This way the central qubit replicas evolve identically in all subsystems. Due to the extended Lindbladian, spin momentum is conserved no longer in the subsystems but in the whole cluster approximation. This way the cluster approximation together with the Lindbladian coupling describes the dynamics of the whole qubit-spin bath system approximately. Considering the dynamics, the main approximation of the method is the neglect of the intra spin bath coupling and entanglement that may affect the dynamics of the central qubit through spin diffusion as well as constructive and destructive interference that can give rise to echo signals and dark states, respectively. These limitations can, however, be systematically lifted by including more and more environmental spins in the subsystems of the cluster approximation.

Figure 10: 
Cluster approximations of a many-spin system S$\mathcal{S}$. (a) S$\mathcal{S}$ consists of a central spin s0 and number of n coupled spins s
i
 that couple only to the central spin s0. (b) First-order cluster approximation of S$\mathcal{S}$ (CCE1) that comprises n + 1 cluster systems c0 and c
i
. c0 includes the central spin s0 only, while c
i
 for i ≠ 0 includes a pair of spins, s0 and s
i
. (c) Second-order cluster approximation of S$\mathcal{S}$ (CCE2) that comprises n/2 + 1 cluster systems c
I
, where each cluster system contains s0 and two coupled spins s
I
 where 1 ≤ I ≤ n/2. c0 includes solely the central spin s0.
Figure 10:

Cluster approximations of a many-spin system S. (a) S consists of a central spin s0 and number of n coupled spins s i that couple only to the central spin s0. (b) First-order cluster approximation of S (CCE1) that comprises n + 1 cluster systems c0 and c i . c0 includes the central spin s0 only, while c i for i ≠ 0 includes a pair of spins, s0 and s i . (c) Second-order cluster approximation of S (CCE2) that comprises n/2 + 1 cluster systems c I , where each cluster system contains s0 and two coupled spins s I where 1 ≤ In/2. c0 includes solely the central spin s0.

In the CCE approach, the corresponding Hilbert-space can be significantly truncated that are coupled to each other in which the density matrix of the central spin can be consistently calculated. As we increase the order of expansion, the results should converge to the theoretical limit, in good analog to the CI expansion method for approaching the accurate correlation energy of the many-electron system. For instance, Yang and co-workers found for the diamond NV centre in the nuclear spin bath of remote 13C nuclear spins [146] that gCCE4 and gCCE5 results agree, thus gCCE4 can considered as absolutely convergent for this particular system.

Seo and co-workers [163] applied CCE method for diamond NV centre and divacancy qubits in 4H SiC for understanding the spin dynamics between the qubit’s electron spin and the nuclear spin bath with assuming natural abundance of 13C at 1.1% and 29Si at 4.5%. The authors ignored the Fermi-contact term in the hyperfine interaction between the electron spin and the nuclear spins. A constant magnetic field was applied in the simulation. We note that because not the gCCE method was applied in Ref. [163], therefore, these simulations could not well describe the spin dynamics for the magnetic fields at GSLAC and ESLAC conditions of the systems, which results in a rapid decrease of the coherence times. This was later done by gCCE method for divacancy qubits in 4H SiC [164]. Seo and co-workers found that CCE2 level of theory well converges with the afore-mentioned conditions and the radius of the spin bath at around 50 Å from the defect qubit’s spin provides convergent results. At the CCE2 level, the distance between interacting nuclear spins was set to 8 Å which converged well [163]. It was found that ensemble averages over 50 samples are good enough to produce numerically converged results. For the magnetic fields above 30 mT they found a simple relation between the T2* times and the concentrations of 13C and 29Si isotopes in the SiC crystal showing that the 13C and 29Si nuclear spins are completely decoupled due to the different Zeeman-splitting for these two species which hinders spin flip-flop processes between these types of nuclear spins. As a consequence, the spin dephasing times of divacancy qubits in 4H SiC are equal or greater to that of diamond NV centre despite the larger density of I = 1/2 nuclear spins in SiC than that in diamond [163]. This conclusion was earlier achieved in a similar study applied to the Si-vacancy qubit (S = 3/2) in 4H SiC [165] as we quote below. For a given cut-off distance R c , the number of total nuclear spins (NC + NSi) surrounding the Si-vacancy centre in SiC is about 2.6 times as large as that of the diamond NV centre, due to the larger 29Si abundance. The reasons for the T2* time of Si-vacancy centres not being reduced by the same factor are as follows: (i) the C-Si bond length 1.89 Å in SiC is larger than the C–C bond length 1.54 Å in diamond, which implies the volume density of nuclear spins is reduced by a factor of (1.89/1.54)3 = 1.8; and (ii) about 80% of the nuclear spins in the bath are 29Si, which have smaller gyromagnetic ratios than 13C (|γC/γSi| ≈ 1.3) and, as a result, produce weaker hyperfine fluctuations. These two factors compensate the larger natural abundance of the 29Si, and results in similar T2* times of defect qubits in SiC and NV centres in diamond.

In a subsequent publication, the spin dephasing times were calculated for hypothetical defect spins with no Fermi-contact hyperfine interaction with moderate magnetic fields (0.1–0.5 mT) in 2D materials [166] by applying the same method. The subject of this study was later extended to 12,000 materials in which both spin dephasing (free-induction decay) and spin coherence (Hahn-echo) times were considered at large constant magnetic fields (e.g., 5 T) [13]. With these simulation conditions, Seo and co-workers [166] found for hBN and molybdenum-disuplhide (MoS2) materials that the spin dephasing time in bulk hBN should be around 18 μs for natural abundance of 11B and 10B isotopes whereas it is about 1.18 ms in MoS2. They attributed the orders of magnitude difference to partially to the variant of nuclear spin density in the two materials and partially to the relatively small gyromagnetic constant of Mo isotopes [166]. By replacing all 11B by 10B should result in >2 factor improvement in the spin dephasing time in hBN, according to the model. Kanai and co-workers [13] concluded that SiC and Si could reach T2 ≈ 5 ms beating even diamond (4 ms). They identified various chalgogenides which have very long intrinsic spin-spin related T2 times at magnetic field of 5 T such as CeO2 (179 ms), CaO (77 ms), α-quartz (8.5 ms), wurtzite ZnO (4.2 ms), and MgO (1.33 ms). It should be mentioned that the T2* time can be long in a common 2D material, WS2 (11 ms). These results should be interpreted with the caveat that the temperature-dependent spin-lattice relaxation times and other limiting factors are not included in that study, thus these results might be valid at relatively low operation temperatures for many materials.

Although, the results on spin dephasing and coherence times with hypothetical defects are somewhat indicative for classifying materials, the defect qubit’s spin relaxation properties may crucially depend on the local environment induced by the defect in terms of ZFS, strain fields and spin density distribution. The last entity is in particular important for materials with dense nuclear spin bath. In that case, the Fermi-contact term in the hyperfine tensor is dominant, and such an effect cannot be fully neglected even in diamond or SiC with relatively dilute nuclear spin densities.

For S = 1 defects the case of GSLAC or ESLAC condition was already mentioned where a simple external parameter, magnitude and direction of the constant magnetic field, may drastically change the coherence properties of the defect qubits, e.g., the interplay between the actual D-constant of the defect spin and the strength of the constant magnetic field [62, 147, 148, 164]. Another interesting example is the so-called clock-transition quantum optics protocol which may be realized by low-symmetry defect spins, e.g., basal divacancy defects in 4H SiC, where the |+1 z ⟩ and |−1 z ⟩ levels naturally splits (see Refs. [164, 167] and references therein). Combining the E = 18.4 MHz splitting with a strong longitudinal splitting (D = 1.334 GHz), the ZFS tensor leads to an avoided crossing of electron spin levels at zero magnetic field from which a clock transition emerges. The qubit levels at the clock transition correspond to |+=(1/2)(|1z+|1z) and |0⟩ = |0 z ⟩ (e.g., see Ref. [167]). The frequency of clock transitions is insensitive to magnetic fields to first order, thus increasing protection from the nuclear bath induced decoherence. Onizhuk and co-workers proved by gCCE theory that the clock transition can indeed elongate the coherence times and the coherence times can be further elongated with opening the gap between |+1⟩ and |−1⟩ levels which indeed occurs for the other basal plane divacancy qubit with E = 82.0 MHz and D = 1.222 GHz ZFS parameters [164]. In experiments, the fluctuating electric fields can lead to a serious decoherence for |+⟩ state, so the charges should be depleted to observe the predicted improvement in the coherence times that was achieved by applying an external electric field to the system [164, 168]. In the simulations, the fully convergent results were achieved at gCCE4 level.

Another important defect spin is the negatively charged boron-vacancy (VB) in hBN that was already mentioned in our review paper. The defect has three nitrogen dangling bonds with large (47 MHz) hyperfine coupling to the electron spin (S = 1) with ZFS in the GHz region both for the ground and excited states (see Ref. [27] and references therein). In other words, the three 14N nuclear spins are strongly coupled to the electron spin for which the Fermi-contact term is dominant. This explains a recent observation [24] that the spin coherence time of VB in hBN is much shorter at the condition of close to zero magnetic fields than that previously anticipated, and it is below 0.1 μs, c.f., 18 μs in Ref. [166] for a hypothetical defect. Ivády and Gali carried out gCCE simulations on this system at the experimental external magnetic field of 14 mT24. It was found that gCCE3 level is convergent. The maximum distance between the defect spin and the nuclear spin should be around 10 Å whereas the maximum distance between the nuclear spins should be around 7.5 Å, in order to achieve convergent results. The simulated spin echo decay curves of the centre in h11BN and h10BN are fitted with a stretched exponential function, exp(t/T2)n, leading to T2 = 92 ns and T2 = 115 ns for h11BN and h10BN, respectively, with an exponent n ≈ 1.35. The dependence of the coherence time with the isotopic content exhibits a linear increase with the 10B abundance. This effect results from the reduced nuclear gyromagnetic ratio of 10B that weakens the hyperfine interaction and the boron nuclear spin flip-flop rate, both of which has a positive impact on the coherence time of the central spin. In a recent study it has been shown [169] that the long T2 ≈ 27 μs can be retained for VB in hBN at giant external magnetic fields, e.g., 3 T, which suppresses the strong electron-nuclear spin couplings.

At moderate external magnetic fields, it is challenging to observe the Rabi-oscillation of the VB electron spin because of the intrinsically short coherence times. Ivády and Gali showed by gCCE simulations [25] that if the microwave field is tuned at the centre hyperfine peak or the m I = ±1 hyperfine peaks of the strongly coupled three 14N spins then the Rabi-oscillation of the 4-spin VB becomes observable, and even 10 MHz detuning significantly suppresses the amplitude of the Rabi-oscillation. In experiments [25], one can see a multiple-frequency oscillation, in which a beat is clearly recognized, and it is superposed on another slow oscillation. The results could be interpreted by the results from gCCE simulations: the nearest neighbour 14N nuclear spins are driven by the microwave field at the given magnetic field (44 mT) and the observed multiple frequencies in Ramsey interference correspond to spin-rotation frequencies in rotating frame on three hyperfine levels. The three frequencies in the spectrum can be identified as the detuning between the microwave field and the centre hyperfine level, as well as the m I = ±1 levels [25]. The spin bath of the 4-spin VB system was further simulated by taking 127 14N and 127 11B nuclear spins into account in which the HSE06 hyperfine tensor and electric field gradient tensors were applied. In the simulations, an effective spin-polarisation transfer could be observed towards these 14N nuclear spins whereas the spin-polarisation towards 11B is small. Hence, it was concluded that the neighbour 14N nuclear spins are responsible for the modulation of the Rabi oscillation, including the decay of the background beyond 0.2 μs, and rest of the spin bath is responsible for the decay of the Rabi oscillation [25].

In pulsed electron spin resonance measurements, the 4-spin nature of VB was further confirmed [26]. The measurements were carried out in a W-band (94 GHz) microwave resonator which brings the electron spin resonance frequencies at around 3.5 T. At this high magnetic fields, the T2 time of VB extends to 15.1 μs observed by electron spin echo (ESE) measurements [26], in good agreement with the recent gCCE3 simulations [169]. The decay curve reveals its oscillatory behavior especially pronounced at the very beginning of the transient curve. Such oscillations refer to electron spin echo envelop modulation (ESEEM) and manifest the presence of coherent coupling of the VB electron spin with magnetic moments of nuclei available in the hBN lattice. The observed beating frequencies corresponds to the nuclear magnetic resonance frequencies (consisting of the combination of hyperfine and quadrupole interactions) of the three nearest neighbour 14N nuclear spins. PBE DFT calculations confirmed that no other 14N near the vacancy could produce such nuclear spin resonances [26] which makes the 4-spin VB model consistent. Finally, it was found that the optical nuclear spin-polarisation at the GSLAC and ESLAC conditions of the external magnetic fields at respective 124 mT and 74 mT can be efficiently carried out towards the neighbour 14N nuclear spins and it can be coherently driven at ≈5 MHz which is much faster than the appropriate nuclear spin resonance frequency [27]. To our knowledge, the spin dephasing and spin-echo simulations for these strongly coupled systems have not yet been reported.

4 Summary

In this paper, we reviewed the recent advances on ab initio theory on defect qubits. A strong emphasis was put on the calculation of excited states, photo-ionisation thresholds and optical excitation spectra also as a function of temperature. A novel theory has been developed on the effective mass states of the excited states of deep defects. Major breakthroughs have been presented on the calculation spin dynamics of the defect qubits which converted the phenomenological description of the spin relaxation times to fully ab initio solution.


Corresponding author: Ádám Gali, Wigner Research Centre for Physics, PO. Box 49, Budapest H-1525, Hungary; and Department of Atomic Physics, Institute of Physics, Budapest University of Technology and Economics, Műegyetem rakpart 3., Budapest H-1111, Hungary, E-mail: .

Funding source: European Commission

Award Identifier / Grant number: 862721

Funding source: European Commission

Award Identifier / Grant number: 101046911

Funding source: National Research, Development and Innovation Office of Hungary

Award Identifier / Grant number: KKP129866

Funding source: National Research, Development and Innovation Office of Hungary

Award Identifier / Grant number: 2022-2.1.1-NL-2022-00004

Funding source: European Commission and National Research, Development, and Innovation Office of Hungary

Award Identifier / Grant number: QuantERA II MAESTRO

  1. Author contributions: All the authors have accepted responsibility for the entire content of this submitted manuscript and approved submission.

  2. Research funding: A.G. acknowledges the Hungarian NKFIH grants No. KKP129866 of the National Excellence Program of Quantum-coherent materials project, the EU QuantERA II MAESTRO project, EU HE FETOPEN project QuanTELCO (Grant No. 862721), the EU HE EIC Pathfinder project QuMicro (Grant No. 101046911) as well as the Quantum Information National Laboratory sponsored by the Innovation and Cultural Ministry (Grant No. 2022-2.1.1-NL-2022-00004).

  3. Conflict of interest statement: The authors declare no conflicts of interest regarding this article.

References

[1] J. N. Eckstein and J. Levy, “Materials issues for quantum computation,” MRS Bull., vol. 38, p. 783, 2013. https://doi.org/10.1557/mrs.2013.210.Search in Google Scholar

[2] G. Wolfowicz, F. J. Heremans, C. P. Anderson, et al.., “Quantum guidelines for solid-state spin defects,” Nat. Rev. Mater., vol. 6, p. 906, 2021. https://doi.org/10.1038/s41578-021-00306-y.Search in Google Scholar

[3] B. E. Kane, “A silicon-based nuclear spin quantum computer,” Nature, vol. 393, p. 133, 1998. https://doi.org/10.1038/30156.Search in Google Scholar

[4] A. Morello, J. J. Pla, F. A. Zwanenburg, et al.., “Single-shot readout of an electron spin in silicon,” Nature, vol. 467, p. 687, 2010. https://doi.org/10.1038/nature09392.Search in Google Scholar PubMed

[5] A. Gruber, A. Dräbenstedt, C. Tietz, L. Fleury, J. Wrachtrup, and C. von Borczyskowski, “Scanning confocal optical microscopy and magnetic resonance on single defect centers,” Science, vol. 276, p. 2012, 1997. https://doi.org/10.1126/science.276.5321.2012.Search in Google Scholar

[6] E. Bourgeois, A. Jarmola, P. Siyushev, et al.., “Photoelectric detection of electron spin resonance of nitrogen-vacancy centres in diamond,” Nat. Commun., vol. 6, p. 8577, 2015. https://doi.org/10.1038/ncomms9577.Search in Google Scholar PubMed PubMed Central

[7] P. Siyushev, M. Nesladek, E. Bourgeois, et al.., “Photoelectrical imaging and coherent spin-state readout of single nitrogen-vacancy centers in diamond,” Science, vol. 363, p. 728, 2019. https://doi.org/10.1126/science.aav2789.Search in Google Scholar PubMed

[8] J. J. L. Morton, A. M. Tyryshkin, R. M. Brown, et al.., “Solid-state quantum memory using the 31P nuclear spin,” Nature, vol. 455, p. 1085, 2008. https://doi.org/10.1038/nature07295.Search in Google Scholar

[9] F. Jelezko, T. Gaebel, I. Popa, M. Domhan, A. Gruber, and J. Wrachtrup, “Observation of coherent oscillation of a single nuclear spin and realization of a two-qubit conditional quantum gate,” Phys. Rev. Lett., vol. 93, p. 130501, 2004. https://doi.org/10.1103/physrevlett.93.130501.Search in Google Scholar

[10] J. R. Weber, W. F. Koehl, J. B. Varley, et al.., “Quantum computing with defects,” Proc. Natl. Acad. Sci., vol. 107, p. 8513, 2010. https://doi.org/10.1073/pnas.1003052107.Search in Google Scholar PubMed PubMed Central

[11] G. Zhang, Y. Cheng, J. P. Chou, and A. Gali, “Material platforms for defect qubits and single-photon emitters,” Appl. Phys. Rev., vol. 7, p. 031308, 2020. https://doi.org/10.1063/5.0006075.Search in Google Scholar

[12] A. M. Ferrenti, N. P. de Leon, J. D. Thompson, and R. J. Cava, “Identifying candidate hosts for quantum defects via data mining,” Npj Comput. Mater., vol. 6, p. 126, 2020. https://doi.org/10.1038/s41524-020-00391-7.Search in Google Scholar

[13] S. Kanai, F. J. Heremans, H. Seo, et al.., “Generalized scaling of spin qubit coherence in over 12,000 host materials,” Proc. Natl. Acad. Sci., vol. 119, p. e2121808119, 2022. https://doi.org/10.1073/pnas.2121808119.Search in Google Scholar PubMed PubMed Central

[14] O. L. Hebnes, M. E. Bathen, Ø. S. Schøyen, S. G. Winther-Larsen, L. Vines, and M. Hjorth-Jensen, “Predicting solid state material platforms for quantum technologies,” Npj Comput. Mater., vol. 8, p. 207, 2022. https://doi.org/10.1038/s41524-022-00888-3.Search in Google Scholar

[15] J. Davidsson, V. Ivády, R. Armiento, and I. A. Abrikosov, “ADAQ: automatic workflows for magneto-optical properties of point defects in semiconductors,” Comput. Phys. Commun., vol. 269, p. 108091, 2021. https://doi.org/10.1016/j.cpc.2021.108091.Search in Google Scholar

[16] F. Bertoldo, S. Ali, S. Manti, and K. S. Thygesen, “Quantum point defects in 2D materials - the QPOD database,” Npj Comput. Mater., vol. 8, p. 56, 2022. https://doi.org/10.1038/s41524-022-00730-w.Search in Google Scholar

[17] W. Redjem, A. Durand, T. Herzig, et al.., “Single artificial atoms in silicon emitting at telecom wavelengths,” Nat. Electron., vol. 3, p. 738, 2020. https://doi.org/10.1038/s41928-020-00499-0.Search in Google Scholar

[18] M. Hollenbach, M. Hollenbach, M. Hollenbach, et al.., “Engineering telecom single-photon emitters in silicon for scalable quantum photonics,” Opt. Express, vol. 28, p. 26111, 2020. https://doi.org/10.1364/oe.397377.Search in Google Scholar PubMed

[19] P. Udvarhelyi, B. Somogyi, G. Thiering, and A. Gali, “Identification of a telecom wavelength single photon emitter in silicon,” Phys. Rev. Lett., vol. 127, p. 196402, 2021. https://doi.org/10.1103/physrevlett.127.196402.Search in Google Scholar

[20] Y. Baron, A. Durand, P. Udvarhelyi, et al.., “Detection of single W-centers in silicon,” ACS Photonics, vol. 9, p. 2337, 2022. https://doi.org/10.1021/acsphotonics.2c00336.Search in Google Scholar

[21] D. B. Higginbottom, A. T. K. Kurkjian, C. Chartrand, et al.., “Optical observation of single spins in silicon,” Nature, vol. 607, p. 266, 2022. https://doi.org/10.1038/s41586-022-04821-y.Search in Google Scholar PubMed

[22] A. Gritsch, L. Weiss, J. Früh, S. Rinner, and A. Reiserer, “Narrow optical transitions in erbium-implanted silicon waveguides,” Phys. Rev. X, vol. 12, p. 041009, 2022. https://doi.org/10.1103/physrevx.12.041009.Search in Google Scholar

[23] A. Ulanowski, B. Merkel, and A. Reiserer, “Spectral multiplexing of telecom emitters with stable transition frequency,” Sci. Adv., vol. 8, p. eabo4538, 2022. https://doi.org/10.1126/sciadv.abo4538.Search in Google Scholar PubMed PubMed Central

[24] A. Haykal, R. Tanos, N. Minotto, et al.., “Decoherence of VB− spin defects in monoisotopic hexagonal boron nitride,” Nat. Commun., vol. 13, p. 4347, 2022. https://doi.org/10.1038/s41467-022-31743-0.Search in Google Scholar PubMed PubMed Central

[25] W. Liu, V. Ivády, Z. P. Li, et al.., “Coherent dynamics of multi-spin VB− center in hexagonal boron nitride,” Nat. Commun., vol. 13, p. 5713, 2022. https://doi.org/10.1038/s41467-022-33399-2.Search in Google Scholar PubMed PubMed Central

[26] F. F. Murzakhanov, G. V. Mamin, S. B. Orlinskii, et al.., “Electron–nuclear coherent coupling and nuclear spin readout through optically polarized VB– spin states in hBN,” Nano Lett., vol. 22, p. 2718, 2022. https://doi.org/10.1021/acs.nanolett.1c04610.Search in Google Scholar PubMed

[27] X. Gao, S. Vaidya, K. Li, et al.., “Nuclear spin polarization and control in hexagonal boron nitride,” Nat. Mater., vol. 21, p. 1024, 2022. https://doi.org/10.1038/s41563-022-01329-8.Search in Google Scholar PubMed

[28] P. Deák, P. Udvarhelyi, G. Thiering, and A. Gali, Nat. Commun., vol. 14, p. 361, 2023. https://doi.org/10.1038/s41467-023-36090-2.Search in Google Scholar PubMed PubMed Central

[29] G. Thiering and A. Gali, “Theory of the optical spin-polarization loop of the nitrogen-vacancy center in diamond,” Phys. Rev. B, vol. 98, p. 085207, 2018. https://doi.org/10.1103/physrevb.98.085207.Search in Google Scholar

[30] A. Gali, “Ab initio theory of the nitrogen-vacancy center in diamond,” Nanophotonics, vol. 8, p. 1907, 2019. https://doi.org/10.1515/nanoph-2019-0154.Search in Google Scholar

[31] M. C. Cambria, A. Gardill, Y. Li, A. Norambuena, J. R. Maze, and S. Kolkowitz, “State-dependent phonon-limited spin relaxation of nitrogen-vacancy centers,” Phys. Rev. Res., vol. 3, p. 013123, 2021. https://doi.org/10.1103/physrevresearch.3.013123.Search in Google Scholar

[32] W. Kohn and L. J. Sham, “Self-consistent equations including exchange and correlation effects,” Phys. Rev. Lett., vol. 140, p. A1133, 1965. https://doi.org/10.1103/physrev.140.a1133.Search in Google Scholar

[33] P. E. Blöchl, “Projector augmented-wave method,” Phys. Rev. B, vol. 50, p. 17953, 1994. https://doi.org/10.1103/physrevb.50.17953.Search in Google Scholar PubMed

[34] C. Freysoldt, B. Grabowski, T. Hickel, et al.., “First-principles calculations for point defects in solids,” Rev. Mod. Phys., vol. 86, p. 253, 2014. https://doi.org/10.1103/revmodphys.86.253.Search in Google Scholar

[35] C. E. Dreyer, A. Alkauskas, J. L. Lyons, A. Janotti, and C. G. Van de Walle, “First-principles calculations of point defects for quantum technologies,” Annu. Rev. Mater. Res., vol. 48, p. 1, 2018. https://doi.org/10.1146/annurev-matsci-070317-124453.Search in Google Scholar

[36] G. Onida, L. Reining, and A. Rubio, “Electronic excitations: density-functional versus many-body green’s-function approaches,” Rev. Mod. Phys., vol. 74, p. 601, 2002. https://doi.org/10.1103/revmodphys.74.601.Search in Google Scholar

[37] Y. Ma, M. Rohlfing, and A. Gali, “Excited states of the negatively charged nitrogen-vacancy color center in diamond,” Phys. Rev. B, vol. 81, p. 041204, 2010. https://doi.org/10.1103/physrevb.81.041204.Search in Google Scholar

[38] S. Choi, M. Jain, and S. G. Louie, “Mechanism for optical initialization of spin in NV−center in diamond,” Phys. Rev. B, vol. 86, p. 041202, 2012. https://doi.org/10.1103/physrevb.86.041202.Search in Google Scholar

[39] G. Kresse and J. Hafner, “Ab initio molecular dynamics for liquid metals,” Phys. Rev. B, vol. 47, p. 558, 1993. https://doi.org/10.1103/physrevb.47.558.Search in Google Scholar PubMed

[40] G. Kresse and J. Furthmüller, “Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set,” Phys. Rev. B, vol. 54, p. 11169, 1996. https://doi.org/10.1103/physrevb.54.11169.Search in Google Scholar PubMed

[41] M. Bockstedte, F. Schütz, T. Garratt, V. Ivády, and A. Gali, “Ab initio description of highly correlated states in defects for realizing quantum bits,” Npj Quantum Mater., vol. 3, p. 31, 2018. https://doi.org/10.1038/s41535-018-0103-6.Search in Google Scholar

[42] A. V. Krukau, O. A. Vydrov, A. F. Izmaylov, and G. E. Scuseria, “Influence of the exchange screening parameter on the performance of screened hybrid functionals,” J. Chem. Phys., vol. 125, p. 224106, 2006. https://doi.org/10.1063/1.2404663.Search in Google Scholar PubMed

[43] P. Giannozzi, S. Baroni, N. Bonini, et al.., “Quantum espresso: a modular and open-source software project for quantum simulations of materials,” J. Phys. Condens. Matter, vol. 21, p. 395502, 2009. https://doi.org/10.1088/0953-8984/21/39/395502.Search in Google Scholar PubMed

[44] H. Ma, M. Govoni, and G. Galli, “Quantum simulations of materials on near-term quantum computers,” Npj Comput. Mater., vol. 6, p. 85, 2020. https://doi.org/10.1038/s41524-020-00353-z.Search in Google Scholar

[45] W. Pfäffle, D. Antonov, J. Wrachtrup, and G. Bester, “Screened configuration interaction method for open-shell excited states applied to NV centers,” Phys. Rev. B, vol. 104, p. 104105, 2021. https://doi.org/10.1103/physrevb.104.104105.Search in Google Scholar

[46] R. Resta, “Thomas-Fermi dielectric screening in semiconductors,” Phys. Rev. B, vol. 16, p. 2717, 1977. https://doi.org/10.1103/physrevb.16.2717.Search in Google Scholar

[47] G. Davies, “Dynamic Jahn–Teller distortions at trigonal optical centres in diamond,” J. Phys. C Solid State Phys., vol. 12, p. 2551, 1979. https://doi.org/10.1088/0022-3719/12/13/019.Search in Google Scholar

[48] N. B. Manson, K. Beha, A. Batalov, et al.., “Assignment of the NV0575-nm zero-phonon line in diamond to a2E-2A2transition,” Phys. Rev. B, vol. 87, p. 155209, 2013. https://doi.org/10.1103/physrevb.87.155209.Search in Google Scholar

[49] L. Muechler, D. I. Badrtdinov, A. Hampel, J. Cano, M. Rösner, and C. E. Dreyer, “Quantum embedding methods for correlated excited states of point defects: case studies and challenges,” Phys. Rev. B, vol. 105, p. 235104, 2022. https://doi.org/10.1103/physrevb.105.235104.Search in Google Scholar

[50] N. Sheng, C. Vorwerk, M. Govoni, and G. Galli, “Green’s function formulation of quantum defect embedding theory,” J. Chem. Theory Comput., vol. 18, p. 3512, 2022. https://doi.org/10.1021/acs.jctc.2c00240.Search in Google Scholar PubMed

[51] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, “Maximally localized Wannier functions: theory and applications,” Rev. Mod. Phys., vol. 84, p. 1419, 2012. https://doi.org/10.1103/revmodphys.84.1419.Search in Google Scholar

[52] V. Ivády, R. Armiento, K. Szász, E. Janzén, A. Gali, and I. A. Abrikosov, “Theoretical unification of hybrid-DFT and DFT + U methods for the treatment of localized orbitals,” Phys. Rev. B, vol. 90, p. 035146, 2014. https://doi.org/10.1103/physrevb.90.035146.Search in Google Scholar

[53] K. Haule, “Exact double counting in combining the dynamical mean field theory and the density functional theory,” Phys. Rev. Lett., vol. 115, p. 196403, 2015. https://doi.org/10.1103/physrevlett.115.196403.Search in Google Scholar

[54] J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized gradient approximation made simple,” Phys. Rev. Lett., vol. 77, p. 3865, 1996. https://doi.org/10.1103/physrevlett.77.3865.Search in Google Scholar

[55] J. H. Skone, M. Govoni, and G. Galli, “Nonempirical range-separated hybrid functionals for solids and molecules,” Phys. Rev. B, vol. 93, p. 235106, 2016. https://doi.org/10.1103/physrevb.93.235106.Search in Google Scholar

[56] G. Thiering and A. Gali, “The (eg ⊗ eu) ⊗ Eg product Jahn–Teller effect in the neutral group-IV vacancy quantum bits in diamond,” Npj Comput. Mater., vol. 5, p. 18, 2019. https://doi.org/10.1038/s41524-019-0158-3.Search in Google Scholar

[57] G. Thiering and A. Gali, Semiconductors and Semimetals, vol. 103, C. E. Nebel, I. Aharonovich, N. Mizuochi, and M. Hatano, Eds., Cambridge, MA, Academic Press, 2020, pp. 1–36.Search in Google Scholar

[58] G. Barcza, V. Ivády, T. Szilvási, et al.., “DMRG on top of plane-wave Kohn–Sham orbitals: a case study of defected boron nitride,” J. Chem. Theory Comput., vol. 17, p. 1143, 2021. https://doi.org/10.1021/acs.jctc.0c00809.Search in Google Scholar PubMed

[59] V. Ivády, G. Barcza, G. Thiering, et al.., “Ab initio theory of the negatively charged boron vacancy qubit in hexagonal boron nitride,” Npj Comput. Mater., vol. 6, p. 41, 2020. https://doi.org/10.1038/s41524-020-0305-x.Search in Google Scholar

[60] M. Abdi, J. P. Chou, A. Gali, and M. B. Plenio, “Color centers in hexagonal boron nitride monolayers: a group theory and ab initio analysis,” ACS Photonics, vol. 5, p. 1967, 2018. https://doi.org/10.1021/acsphotonics.7b01442.Search in Google Scholar

[61] A. Gottscholl, M. Kianinia, V. Soltamov, et al.., “Initialization and read-out of intrinsic spin defects in a van der Waals crystal at room temperature,” Nat. Mater., vol. 19, p. 540, 2020. https://doi.org/10.1038/s41563-020-0619-6.Search in Google Scholar PubMed

[62] V. Ivády, “Longitudinal spin relaxation model applied to point-defect qubit systems,” Phys. Rev. B, vol. 101, p. 155203, 2020. https://doi.org/10.1103/PhysRevB.101.155203.Search in Google Scholar

[63] J. R. Reimers, A. Sajid, R. Kobayashi, and M. J. Ford, “Convergence of defect energetics calculations,” J. Phys. Chem. C., vol. 124, p. 21178, 2020. https://doi.org/10.1021/acs.jpcc.0c06445.Search in Google Scholar

[64] J. R. Reimers, J. Shen, M. Kianinia, et al.., “Photoluminescence, photophysics, and photochemistry of the VB− defect in hexagonal boron nitride,” Phys. Rev. B, vol. 102, p. 144105, 2020. https://doi.org/10.1103/physrevb.102.144105.Search in Google Scholar

[65] A. Pershin, G. Barcza, O. Legeza, and A. Gali, “Highly tunable magneto-optical response from magnesium-vacancy color centers in diamond,” Npj Quantum Inf., vol. 7, p. 99, 2021. https://doi.org/10.1038/s41534-021-00439-6.Search in Google Scholar

[66] T. Lühmann, R. John, R. Wunderlich, J. Meijer, and S. Pezzagna, “Coulomb-driven single defect engineering for scalable qubits and spin sensors in diamond,” Nat. Commun., vol. 10, p. 4956, 2019. https://doi.org/10.1038/s41467-019-12556-0.Search in Google Scholar PubMed PubMed Central

[67] G. H. Booth, A. Grüneis, G. Kresse, and A. Alavi, “Towards an exact description of electronic wavefunctions in real solids,” Nature, vol. 493, p. 365, 2013. https://doi.org/10.1038/nature11770.Search in Google Scholar PubMed

[68] W. Hu and G. K. L. Chan, “Excited-state geometry optimization with the density matrix renormalization group, as applied to polyenes,” J. Chem. Theory Comput., vol. 11, p. 3000, 2015. https://doi.org/10.1021/acs.jctc.5b00174.Search in Google Scholar PubMed

[69] E. Runge and E. K. U. Gross, “Density-functional theory for time-dependent systems,” Phys. Rev. Lett., vol. 52, p. 997, 1984. https://doi.org/10.1103/physrevlett.52.997.Search in Google Scholar

[70] M. E. Casida, Recent Advances in Density Functional Theory, Singapore, World Scientific, 1995, p. 155.10.1142/9789812830586_0005Search in Google Scholar

[71] A. Gali, “Time-dependent density functional study on the excitation spectrum of point defects in semiconductors,” Phys. Status Solidi B, vol. 248, p. 1337, 2011. https://doi.org/10.1002/pssb.201046254.Search in Google Scholar

[72] A. Gali, “Excitation spectrum of point defects in semiconductors studied by time-dependent density functional theory,” J. Mater. Res., vol. 27, p. 897, 2012. https://doi.org/10.1557/jmr.2011.431.Search in Google Scholar

[73] I. I. Vlasov, A. A. Shiryaev, T. Rendler, et al.., “Molecular-sized fluorescent nanodiamonds,” Nat. Nanotechnol., vol. 9, p. 54, 2014. https://doi.org/10.1038/nnano.2013.255.Search in Google Scholar PubMed

[74] A. Gali, E. Janzén, P. Deák, G. Kresse, and E. Kaxiras, “Theory of spin-conserving excitation of the N−V−Center in diamond,” Phys. Rev. Lett., vol. 103, p. 186404, 2009. https://doi.org/10.1103/physrevlett.103.186404.Search in Google Scholar PubMed

[75] Y. Jin, M. Govoni, G. Wolfowicz, et al.., “Photoluminescence spectra of point defects in semiconductors: validation of first-principles calculations,” Phys. Rev. Mater., vol. 5, p. 084603, 2021. https://doi.org/10.1103/physrevmaterials.5.084603.Search in Google Scholar

[76] Z. Li and W. Liu, “Theoretical and numerical assessments of spin-flip time-dependent density functional theory,” J. Chem. Phys., vol. 136, p. 024107, 2012. https://doi.org/10.1063/1.3676736.Search in Google Scholar PubMed

[77] D. Casanova and A. I. Krylov, “Spin-flip methods in quantum chemistry,” Phys. Chem. Chem. Phys., vol. 22, p. 4326, 2020. https://doi.org/10.1039/c9cp06507e.Search in Google Scholar PubMed

[78] Y. Jin, M. Govoni, and G. Galli, “Vibrationally resolved optical excitations of the nitrogen-vacancy center in diamond,” Npj Comput. Mater., vol. 8, p. 238, 2022. https://doi.org/10.1038/s41524-022-00928-y.Search in Google Scholar

[79] L. Razinkovas, M. W. Doherty, N. B. Manson, C. G. Van de Walle, and A. Alkauskas, “Vibrational and vibronic structure of isolated point defects: the nitrogen-vacancy center in diamond,” Phys. Rev. B, vol. 104, p. 045303, 2021. https://doi.org/10.1103/physrevb.104.045303.Search in Google Scholar

[80] E. Monino and P. F. Loos, “Spin-conserved and spin-flip optical excitations from the bethe–salpeter equation formalism,” J. Chem. Theory Comput., vol. 17, p. 2852, 2021. https://doi.org/10.1021/acs.jctc.1c00074.Search in Google Scholar PubMed PubMed Central

[81] B. A. Barker and D. A. Strubbe, “Spin-flip Bethe–Salpeter equation approach for ground and excited states of open-shell molecules and defects in solids,” 2022, arXiv:2207.04549 [cond-mat, physics:physics].Search in Google Scholar

[82] X. Andrade, D. Strubbe, U. D. Giovannini, et al.., “Real-space grids and the Octopus code as tools for the development of new simulation approaches for electronic systems,” Phys. Chem. Chem. Phys., vol. 17, p. 31371, 2015. https://doi.org/10.1039/c5cp00351b.Search in Google Scholar PubMed

[83] N. Tancogne-Dejean, M. J. T. Oliveira, X. Andrade, et al.., “Octopus, a computational framework for exploring light-driven phenomena and quantum dynamics in extended and finite systems,” J. Chem. Phys., vol. 152, p. 124119, 2020. https://doi.org/10.1063/1.5142502.Search in Google Scholar PubMed

[84] M. Zacharias, C. E. Patrick, and F. Giustino, “Stochastic approach to phonon-assisted optical absorption,” Phys. Rev. Lett., vol. 115, p. 177401, 2015. https://doi.org/10.1103/physrevlett.115.177401.Search in Google Scholar

[85] E. Cannuccia, B. Monserrat, and C. Attaccalite, “Theory of phonon-assisted luminescence in solids: application to hexagonal boron nitride,” Phys. Rev. B, vol. 99, p. 081109, 2019. https://doi.org/10.1103/physrevb.99.081109.Search in Google Scholar

[86] F. Libbi, P. M. M. de Melo, Z. Zanolli, M. J. Verstraete, and N. Marzari, “Phonon-assisted luminescence in defect centers from many-body perturbation theory,” Phys. Rev. Lett., vol. 128, p. 167401, 2022. https://doi.org/10.1103/physrevlett.128.167401.Search in Google Scholar PubMed

[87] N. Aslam, G. Waldherr, P. Neumann, F. Jelezko, and J. Wrachtrup, “Photo-induced ionization dynamics of the nitrogen vacancy defect in diamond investigated by single-shot charge state detection,” New J. Phys., vol. 15, p. 013064, 2013. https://doi.org/10.1088/1367-2630/15/1/013064.Search in Google Scholar

[88] E. Bourgeois, E. Londero, K. Buczak, et al.., “Enhanced photoelectric detection of NV magnetic resonances in diamond under dual-beam excitation,” Phys. Rev. B, vol. 95, p. 041402, 2017. https://doi.org/10.1103/physrevb.95.041402.Search in Google Scholar

[89] L. Razinkovas, M. Maciaszek, F. Reinhard, M. W. Doherty, and A. Alkauskas, “Photoionization of negatively charged NV centers in diamond: theory and ab initio calculations,” Phys. Rev. B, vol. 104, p. 235301, 2021. https://doi.org/10.1103/physrevb.104.235301.Search in Google Scholar

[90] P. Siyushev, H. Pinto, M. Vörös, A. Gali, F. Jelezko, and J. Wrachtrup, “Optically controlled switching of the charge state of a single nitrogen-vacancy center in diamond at cryogenic temperatures,” Phys. Rev. Lett., vol. 110, p. 167402, 2013. https://doi.org/10.1103/physrevlett.110.167402.Search in Google Scholar PubMed

[91] D. Wirtitsch, S. Wachter, G. Reisenbauer, et al.., “Exploiting ionization dynamics in the nitrogen vacancy center for rapid, high-contrast spin, and charge state initialization,” Phys. Rev. Res., vol. 5, p. 013014, 2023. https://doi.org/10.1103/physrevresearch.5.013014.Search in Google Scholar

[92] J. C. Slater, “The theory of complex spectra,” Phys. Rev., vol. 34, p. 1293, 1929. https://doi.org/10.1103/physrev.34.1293.Search in Google Scholar

[93] E. U. Condon, “The theory of complex spectra,” Phys. Rev., vol. 36, p. 1121, 1930. https://doi.org/10.1103/physrev.36.1121.Search in Google Scholar

[94] J. R. Maze, A. Gali, E. Togan, et al.., “Properties of nitrogen-vacancy centers in diamond: the group theoretic approach,” New J. Phys., vol. 13, p. 025025, 2011. https://doi.org/10.1088/1367-2630/13/2/025025.Search in Google Scholar

[95] D. M. Irber, F. Poggiali, F. Kong, et al.., “Robust all-optical single-shot readout of nitrogen-vacancy centers in diamond,” Nat. Commun., vol. 12, p. 532, 2021. https://doi.org/10.1038/s41467-020-20755-3.Search in Google Scholar PubMed PubMed Central

[96] Q. Zhang, Y. Guo, W. Ji, et al.., “High-fidelity single-shot readout of single electron spin in diamond with spin-to-charge conversion,” Nat. Commun., vol. 12, p. 1529, 2021. https://doi.org/10.1038/s41467-021-21781-5.Search in Google Scholar PubMed PubMed Central

[97] P. Deák, B. Aradi, T. Frauenheim, E. Janzén, and A. Gali, “Accurate defect levels obtained from the HSE06 range-separated hybrid functional,” Phys. Rev. B, vol. 81, p. 153203, 2010. https://doi.org/10.1103/physrevb.81.153203.Search in Google Scholar

[98] E. Londero, E. Bourgeois, M. Nesladek, and A. Gali, “Identification of nickel-vacancy defects by combining experimental and ab initio simulated photocurrent spectra,” Phys. Rev. B, vol. 97, p. 241202, 2018. https://doi.org/10.1103/physrevb.97.241202.Search in Google Scholar

[99] N. Mizuochi, T. Makino, H. Kato, et al.., “Electrically driven single-photon source at room temperature in diamond,” Nat. Photonics, vol. 6, p. 299, 2012. https://doi.org/10.1038/nphoton.2012.75.Search in Google Scholar

[100] A. Lozovoi, H. Jayakumar, D. Daw, et al.., “Optical activation and detection of charge transport between individual colour centres in diamond,” Nat. Electron., vol. 4, p. 717, 2021. https://doi.org/10.1038/s41928-021-00656-z.Search in Google Scholar

[101] N. T. Son, P. Carlsson, J. Ul Hassan, et al.., “Divacancy in 4H-SiC,” Phys. Rev. Lett., vol. 96, p. 055501, 2006. https://doi.org/10.1103/physrevlett.96.055501.Search in Google Scholar PubMed

[102] A. Gali, A. Gällström, N. T. Son, and E. Janzén, “Theory of neutral divacancy in SiC: a defect for spintronics,” Mater. Sci. Forum, vols. 645–648, p. 395, 2010. https://doi.org/10.4028/www.scientific.net/msf.645-648.395.Search in Google Scholar

[103] W. F. Koehl, B. B. Buckley, F. J. Heremans, G. Calusine, and D. D. Awschalom, “Room temperature coherent control of defect spin qubits in silicon carbide,” Nature, vol. 479, p. 84, 2011. https://doi.org/10.1038/nature10562.Search in Google Scholar PubMed

[104] A. L. Falk, B. B. Buckley, G. Calusine, et al.., “Polytype control of spin qubits in silicon carbide,” Nat. Commun., vol. 4, p. 1819, 2013. https://doi.org/10.1038/ncomms2854.Search in Google Scholar PubMed PubMed Central

[105] Q. Li, J. F. Wang, F. F. Yan, et al.., “Room-temperature coherent manipulation of single-spin qubits in silicon carbide with a high readout contrast,” Natl. Sci. Rev., vol. 9, p. nwab122, 2022. https://doi.org/10.1093/nsr/nwab122.Search in Google Scholar PubMed PubMed Central

[106] V. Ivády, J. Davidsson, N. Delegan, et al.., “Stabilization of point-defect spin qubits by quantum wells,” Nat. Commun., vol. 10, p. 5607, 2019. https://doi.org/10.1038/s41467-019-13495-6.Search in Google Scholar PubMed PubMed Central

[107] A. Csóré and A. Gali, Wide Bandgap Semiconductors for Power Electronics, Hoboken, New Jersey, USA, John Wiley & Sons Ltd, 2021, pp. 503–528.10.1002/9783527824724.ch17Search in Google Scholar

[108] D. A. Golter and C. W. Lai, “Optical switching of defect charge states in 4H-SiC,” Sci. Rep., vol. 7, p. 13406, 2017. https://doi.org/10.1038/s41598-017-13813-2.Search in Google Scholar PubMed PubMed Central

[109] G. Wolfowicz, C. P. Anderson, A. L. Yeats, et al.., “Optical charge state control of spin defects in 4H-SiC,” Nat. Commun., vol. 8, p. 1876, 2017. https://doi.org/10.1038/s41467-017-01993-4.Search in Google Scholar PubMed PubMed Central

[110] B. Magnusson, N. T. Son, A. Csóré, et al.., “Excitation properties of the divacancy in 4H-SiC,” Phys. Rev. B, vol. 98, p. 195202, 2018. https://doi.org/10.1103/physrevb.98.195202.Search in Google Scholar

[111] A. Csóré, I. G. Ivanov, N. T. Son, and A. Gali, “Fluorescence spectrum and charge state control of divacancy qubits via illumination at elevated temperatures in 4H silicon carbide,” Phys. Rev. B, vol. 105, p. 165108, 2022. https://doi.org/10.1103/physrevb.105.165108.Search in Google Scholar

[112] A. Marini, “Ab InitioFinite-temperature excitons,” Phys. Rev. Lett., vol. 101, p. 106405, 2008. https://doi.org/10.1103/physrevlett.101.106405.Search in Google Scholar

[113] E. Cannuccia and A. Marini, “Effect of the quantum zero-point atomic motion on the optical and electronic properties of diamond and trans-polyacetylene,” Phys. Rev. Lett., vol. 107, p. 255501, 2011. https://doi.org/10.1103/physrevlett.107.255501.Search in Google Scholar PubMed

[114] A. Gali, T. Demján, M. Vörös, G. Thiering, E. Cannuccia, and A. Marini, “Electron–vibration coupling induced renormalization in the photoemission spectrum of diamondoids,” Nat. Commun., vol. 7, p. 11327, 2016. https://doi.org/10.1038/ncomms11327.Search in Google Scholar PubMed PubMed Central

[115] E. Cannuccia and A. Gali, “Thermal evolution of silicon carbide electronic bands,” Phys. Rev. Mater., vol. 4, p. 014601, 2020. https://doi.org/10.1103/physrevmaterials.4.014601.Search in Google Scholar

[116] A. Alkauskas, B. B. Buckley, D. D. Awschalom, and C. G. V. de Walle, “First-principles theory of the luminescence lineshape for the triplet transition in diamond NV centres,” New J. Phys., vol. 16, p. 073026, 2014. https://doi.org/10.1088/1367-2630/16/7/073026.Search in Google Scholar

[117] M. Chagas da Silva, M. Lorke, B. Aradi, et al.., “Self-consistent potential correction for charged periodic systems,” Phys. Rev. Lett., vol. 126, p. 076401, 2021. https://doi.org/10.1103/physrevlett.126.076401.Search in Google Scholar PubMed

[118] V. Popescu and A. Zunger, “Extracting eversus k→effective band structure from supercell calculations on alloys and impurities,” Phys. Rev. B, vol. 85, p. 085201, 2012. https://doi.org/10.1103/physrevb.85.085201.Search in Google Scholar

[119] M. W. Swift, H. Peelaers, S. Mu, J. J. L. Morton, and C. G. Van de Walle, “First-principles calculations of hyperfine interaction, binding energy, and quadrupole coupling for shallow donors in silicon,” Npj Comput. Mater., vol. 6, p. 181, 2020. https://doi.org/10.1038/s41524-020-00448-7.Search in Google Scholar

[120] J. Mansir, P. Conti, Z. Zeng, et al.., “Linear hyperfine tuning of donor spins in silicon using hydrostatic strain,” Phys. Rev. Lett., vol. 120, p. 167701, 2018. https://doi.org/10.1103/physrevlett.120.167701.Search in Google Scholar PubMed

[121] J. Wagner, K. Thonke, and R. Sauer, “Excitation spectroscopy on the 0.79-eV (C) line defect in irradiated silicon,” Phys. Rev. B, vol. 29, p. 7051, 1984. https://doi.org/10.1103/physrevb.29.7051.Search in Google Scholar

[122] J. Wagner, A. Dörnen, and R. Sauer, “Donorlike excited states of the thermally induced 0.767-eV (Pline) defect in oxygen-rich silicon,” Phys. Rev. B, vol. 31, p. 5561, 1985. https://doi.org/10.1103/physrevb.31.5561.Search in Google Scholar PubMed

[123] J. H. Svensson, B. Monemar, and E. Janzén, “Pseudodonor electronic excited states of neutral complex defects in silicon,” Phys. Rev. Lett., vol. 65, p. 1796, 1990. https://doi.org/10.1103/physrevlett.65.1796.Search in Google Scholar PubMed

[124] P. Udvarhelyi, A. Pershin, P. Deák, and A. Gali, “An L-band emitter with quantum memory in silicon,” Npj Comput. Mater., vol. 8, p. 262, 2022. https://doi.org/10.1038/s41524-022-00957-7.Search in Google Scholar

[125] T. Egilsson, J. P. Bergman, I. G. Ivanov, A. Henry, and E. Janzén, “Properties of the D1 bound exciton in 4H−SiC,” Phys. Rev. B, vol. 59, p. 1956, 1999. https://doi.org/10.1103/physrevb.59.1956.Search in Google Scholar

[126] L. Storasta, F. H. C. Carlsson, S. G. Sridhara, et al.., “Pseudodonor nature of the DI defect in 4H−SiC,” Appl. Phys. Lett., vol. 78, p. 46, 2001. https://doi.org/10.1063/1.1334907.Search in Google Scholar

[127] A. Gali, P. Deák, P. Ordejón, N. T. Son, E. Janzén, and W. J. Choyke, “Aggregation of carbon interstitials in silicon carbide: a theoretical study,” Phys. Rev. B, vol. 68, p. 125201, 2003. https://doi.org/10.1103/physrevb.68.125201.Search in Google Scholar

[128] T. A. G. Eberlein, R. Jones, S. Öberg, and P. R. Briddon, “Density functional theory calculation of the DI optical center in SiC,” Phys. Rev. B, vol. 74, p. 144106, 2006. https://doi.org/10.1103/physrevb.74.144106.Search in Google Scholar

[129] B. C. Rose, D. Huang, Z. H. Zhang, et al.., “Observation of an environmentally insensitive solid-state spin defect in diamond,” Science, vol. 361, p. 60, 2018. https://doi.org/10.1126/science.aao0290.Search in Google Scholar PubMed

[130] E. Londero, G. Thiering, L. Razinkovas, A. Gali, and A. Alkauskas, “Vibrational modes of negatively charged silicon-vacancy centers in diamond from ab initio calculations,” Phys. Rev. B, vol. 98, p. 35306, 2018. https://doi.org/10.1103/PhysRevB.98.035306.Search in Google Scholar

[131] M. J. Rayson and P. R. Briddon, “First principles method for the calculation of zero-field splitting tensors in periodic systems,” Phys. Rev. B, vol. 77, p. 035119, 2008. https://doi.org/10.1103/physrevb.77.035119.Search in Google Scholar

[132] T. Biktagirov, W. G. Schmidt, and U. Gerstmann, “Spin decontamination for magnetic dipolar coupling calculations: application to high-spin molecules and solid-state spin qubits,” Phys. Rev. Res., vol. 2, p. 022024, 2020. https://doi.org/10.1103/physrevresearch.2.022024.Search in Google Scholar

[133] G. Thiering and A. Gali, “Magneto-optical spectra of the split nickel-vacancy defect in diamond,” Phys. Rev. Res., vol. 3, p. 043052, 2021. https://doi.org/10.1103/physrevresearch.3.043052.Search in Google Scholar

[134] T. Biktagirov and U. Gerstmann, “Spin-orbit driven electrical manipulation of the zero-field splitting in high-spin centers in solids,” Phys. Rev. Res., vol. 2, p. 023071, 2020. https://doi.org/10.1103/physrevresearch.2.023071.Search in Google Scholar

[135] A. L. Falk, P. V. Klimov, B. B. Buckley, et al.., “Electrically and mechanically tunable electron spins in silicon carbide color centers,” Phys. Rev. Lett., vol. 112, p. 187601, 2014. https://doi.org/10.1103/physrevlett.112.187601.Search in Google Scholar PubMed

[136] L. Spindlberger, A. Csóré, G. Thiering, et al.., “Optical properties of vanadium in 4H silicon carbide for quantum technology,” Phys. Rev. Appl., vol. 12, p. 014015, 2019. https://doi.org/10.1103/physrevapplied.12.014015.Search in Google Scholar

[137] G. Wolfowicz, C. P. Anderson, B. Diler, O. G. Poluektov, F. J. Heremans, and D. D. Awschalom, “Vanadium spin qubits as telecom quantum emitters in silicon carbide,” Sci. Adv., vol. 6, p. eaaz1192, 2020. https://doi.org/10.1126/sciadv.aaz1192.Search in Google Scholar PubMed PubMed Central

[138] G. Thiering and A. Gali, “Ab initio magneto-optical spectrum of group-IV vacancy color centers in diamond,” Phys. Rev. X, vol. 8, p. 021063, 2018. https://doi.org/10.1103/physrevx.8.021063.Search in Google Scholar

[139] I. Bersurker, The Jahn–Teller Effect, Cambridge, UK, Cambridge University Press, 2006.Search in Google Scholar

[140] F. S. Ham, “Dynamical Jahn–Teller effect in paramagnetic resonance spectra: orbital reduction factors and partial quenching of spin-orbit interaction,” Phys. Rev., vol. 138, p. A1727, 1965. https://doi.org/10.1103/physrev.138.a1727.Search in Google Scholar

[141] A. Csóré and A. Gali, “Ab initio determination of pseudospin for paramagnetic defects in SiC,” Phys. Rev. B, vol. 102, p. 241201, 2020. https://doi.org/10.1103/physrevb.102.241201.Search in Google Scholar

[142] V. Ivády, I. A. Abrikosov, E. Janzén, and A. Gali, “Role of screening in the density functional applied to transition-metal defects in semiconductors,” Phys. Rev. B, vol. 87, p. 205201, 2013. https://doi.org/10.1103/physrevb.87.205201.Search in Google Scholar

[143] W. Gehlhoff and R. N. Pereira, “The identity of the AB1 and NE4 electron paramagnetic resonance spectra in high-pressure high-temperature diamond,” J. Phys.: Condens. Matter, vol. 14, p. 13751, 2002. https://doi.org/10.1088/0953-8984/14/50/304.Search in Google Scholar

[144] W. Yang and R. B. Liu, “Quantum many-body theory of qubit decoherence in a finite-size spin bath,” Phys. Rev. B, vol. 78, p. 085315, 2008. https://doi.org/10.1103/physrevb.78.085315.Search in Google Scholar

[145] W. Yang and R. B. Liu, “Quantum many-body theory of qubit decoherence in a finite-size spin bath. II. Ensemble dynamics,” Phys. Rev. B, vol. 79, p. 115320, 2009. https://doi.org/10.1103/physrevb.79.115320.Search in Google Scholar

[146] Z. S. Yang, Y. X. Wang, M. J. Tao, et al.., “Longitudinal relaxation of a nitrogen-vacancy center in a spin bath by generalized cluster-correlation expansion method,” Ann. Phys., vol. 413, p. 168063, 2020. https://doi.org/10.1016/j.aop.2019.168063.Search in Google Scholar

[147] V. Ivády, H. Zheng, A. Wickenbrock, et al.., “Photoluminescence at the ground-state level anticrossing of the nitrogen-vacancy center in diamond: a comprehensive study,” Phys. Rev. B, vol. 103, p. 035307, 2021. https://doi.org/10.1103/physrevb.103.035307.Search in Google Scholar

[148] M. Gulka, D. Wirtitsch, V. Ivády, et al.., “Room-temperature control and electrical readout of individual nitrogen-vacancy nuclear spins,” Nat. Commun., vol. 12, p. 4421, 2021. https://doi.org/10.1038/s41467-021-24494-x.Search in Google Scholar PubMed PubMed Central

[149] H. Park, J. Lee, S. Han, S. Oh, and H. Seo, “Decoherence of nitrogen-vacancy spin ensembles in a nitrogen electron-nuclear spin bath in diamond,” Npj Quantum Inf., vol. 8, p. 95, 2022. https://doi.org/10.1038/s41534-022-00605-4.Search in Google Scholar

[150] O. Bulancea-Lindvall, N. T. Son, I. A. Abrikosov, and V. Ivády, “Dipolar spin relaxation of divacancy qubits in silicon carbide,” Npj Comput. Mater., vol. 7, p. 213, 2021. https://doi.org/10.1038/s41524-021-00673-8.Search in Google Scholar

[151] O. Bulancea-Lindvall, M. T. Eiles, N. T. Son, I. A. Abrikosov, and V. Ivády, “Isotope purification induced reduction of spin relaxation and spin coherence times in semiconductors,” 2022, arXiv:2205.05105 [cond-mat].10.1103/PhysRevApplied.19.064046Search in Google Scholar

[152] A. Lunghi, “Spin-phonon relaxation in magnetic molecules: theory, predictions and insights,” 2022, arXiv:2202.03776 [cond-mat].10.1007/978-3-031-31038-6_6Search in Google Scholar

[153] A. Jarmola, V. M. Acosta, K. Jensen, S. Chemerisov, and D. Budker, “Temperature- and magnetic-field-dependent longitudinal spin relaxation in nitrogen-vacancy ensembles in diamond,” Phys. Rev. Lett., vol. 108, p. 197601, 2012. https://doi.org/10.1103/physrevlett.108.197601.Search in Google Scholar

[154] A. Norambuena, E. Muñoz, H. T. Dinani, et al.., “Spin-lattice relaxation of individual solid-state spins,” Phys. Rev. B, vol. 97, p. 094304, 2018. https://doi.org/10.1103/physrevb.97.094304.Search in Google Scholar

[155] J. Gugler, T. Astner, A. Angerer, J. Schmiedmayer, J. Majer, and P. Mohn, “Ab initio calculation of the spin lattice relaxation time T1 for nitrogen-vacancy centers in diamond,” Phys. Rev. B, vol. 98, p. 214442, 2018. https://doi.org/10.1103/physrevb.98.214442.Search in Google Scholar

[156] M. C. Cambria, A. Norambuena, H. T. Dinani, et al.., “Temperature-dependent phonon-induced relaxation of the nitrogen-vacancy spin triplet in diamond,” 2022, arXiv:2209.14446 [cond-mat, physics:quant-ph].Search in Google Scholar

[157] A. Lunghi, “Toward exact predictions of spin-phonon relaxation times: an ab initio implementation of open quantum systems theory,” Sci. Adv., vol. 8, p. eabn7880, 2022. https://doi.org/10.1126/sciadv.abn7880.Search in Google Scholar PubMed PubMed Central

[158] V. Ivády, T. Simon, J. R. Maze, I. A. Abrikosov, and A. Gali, “Pressure and temperature dependence of the zero-field splitting in the ground state of NV centers in diamond: a first-principles study,” Phys. Rev. B, vol. 90, p. 235205, 2014. https://doi.org/10.1103/physrevb.90.235205.Search in Google Scholar

[159] H. Tang, A. R. Barr, G. Wang, P. Cappellaro, and J. Li, “First-principles calculation of the temperature-dependent transition energies in spin defects,” 2022, arXiv:2205.02791 [cond-mat, physics:quant-ph].Search in Google Scholar

[160] E. L. Hahn, “Spin echoes,” Phys. Rev., vol. 80, p. 580, 1950. https://doi.org/10.1103/physrev.80.580.Search in Google Scholar

[161] B. A. Myers, A. Ariyaratne, and A. C. B. Jayich, “Double-quantum spin-relaxation limits to coherence of near-surface nitrogen-vacancy centers,” Phys. Rev. Lett., vol. 118, p. 197201, 2017. https://doi.org/10.1103/physrevlett.118.197201.Search in Google Scholar PubMed

[162] A. Gardill, M. Cambria, and S. Kolkowitz, “Fast relaxation on qutrit transitions of nitrogen-vacancy centers in nanodiamonds,” Phys. Rev. Appl., vol. 13, p. 034010, 2020. https://doi.org/10.1103/physrevapplied.13.034010.Search in Google Scholar

[163] H. Seo, A. L. Falk, P. V. Klimov, K. C. Miao, G. Galli, and D. D. Awschalom, “Quantum decoherence dynamics of divacancy spins in silicon carbide,” Nat. Commun., vol. 7, p. 12935, 2016. https://doi.org/10.1038/ncomms12935.Search in Google Scholar PubMed PubMed Central

[164] M. Onizhuk, K. C. Miao, J. P. Blanton, et al.., “Probing the coherence of solid-state qubits at avoided crossings,” PRX Quantum, vol. 2, p. 010311, 2021. https://doi.org/10.1103/prxquantum.2.010311.Search in Google Scholar

[165] L. P. Yang, C. Burk, M. Widmann, S. Y. Lee, J. Wrachtrup, and N. Zhao, “Electron spin decoherence in silicon carbide nuclear spin bath,” Phys. Rev. B, vol. 90, p. 241203, 2014. https://doi.org/10.1103/physrevb.90.241203.Search in Google Scholar

[166] M. Ye, H. Seo, and G. Galli, “Spin coherence in two-dimensional materials,” Npj Comput. Mater., vol. 5, p. 44, 2019. https://doi.org/10.1038/s41524-019-0182-3.Search in Google Scholar

[167] K. C. Miao, A. Bourassa, C. P. Anderson, et al.., “Electrically driven optical interferometry with spins in silicon carbide,” Sci. Adv., vol. 5, p. eaay0527, 2019. https://doi.org/10.1126/sciadv.aay0527.Search in Google Scholar PubMed PubMed Central

[168] A. Bourassa, C. P. Anderson, K. C. Miao, et al.., “Entanglement and control of single nuclear spins in isotopically engineered silicon carbide,” Nat. Mater., vol. 19, p. 1319, 2020. https://doi.org/10.1038/s41563-020-00802-6.Search in Google Scholar PubMed

[169] J. Lee, H. Park, and H. Seo, “First-principles theory of extending the spin qubit coherence time in hexagonal boron nitride,” Npj 2D Mater. Appl., vol. 6, p. 60, 2022. https://doi.org/10.1038/s41699-022-00336-2.Search in Google Scholar

Received: 2022-11-21
Accepted: 2023-01-23
Published Online: 2023-02-01

© 2023 the author(s), published by De Gruyter, Berlin/Boston

This work is licensed under the Creative Commons Attribution 4.0 International License.

Downloaded on 5.3.2024 from https://www.degruyter.com/document/doi/10.1515/nanoph-2022-0723/html?lang=en
Scroll to top button