Charge-State Stability of Color Centers in Wide-Bandgap Semiconductors

The NV$^-$ color center in diamond has been extensively investigated for quantum sensing, computation, and communication applications. Nonetheless, charge-state decay from the NV$^-$ to its neutral counterpart the NV$^0$ detrimentally affects the robustness of the NV$^-$ center and remains to be fully overcome. In this work, we provide an $ab~initio$ formalism for accurately estimating the rate of charge-state decay of color centers in wide-bandgap semiconductors. Our formalism employs density functional theory calculations in the context of thermal equilibrium. We illustrate the method using the transition of NV$^-$ to NV$^0$ in the presence of substitutional N [see Z. Yuan $et~al$., PRR 2, 033263 (2020)].


I. INTRODUCTION
Color centers in wide-bandgap semiconductor hosts have garnered significant interest for potential applications in quantum computation [1][2][3], communication [1,4,5], and sensing [6].The NV − color center in diamond in particular, consisting of a single substitutional N atom adjacent to a C vacancy with an additional negative charge, has enjoyed several research advances including the realization of a coherence time on the order of seconds [7] and entanglement of NV − pairs over a distance greater than a kilometer [8].A drawback in the utilization of NV − centers for computation, communication, and sensing applications, however, is their tendency to revert to the neutral state after optical initialization into the singly negatively charged state [9].The implied large hole-capture cross section of optically activated NV − centers in diamond has been investigated both experimentally and theoretically [9][10][11], where the theory has included a bound-exciton model [10] and semi-classical Monte Carlo simulations [11].Arriving at a first-principles description of carrier capture that can apply in thermal equilibrium or under the action of an external field is still the subject of much work [11][12][13][14][15][16].Herein we provide such a framework which is accurate and demonstrates that ionizing-dopant concentrations along with the electronic structure of the color center of interest and of the ionizing dopant are crucial for the determination of the expected timescale for hole capture by the ionized color centers.To provide an illustration of the method, we investigate charge transfer rates for NV − centers [17][18][19][20][21][22][23][24][25][26] in the presence of N in diamond.

II. COMPUTATIONAL METHODS
We used VASP [65][66][67][68] for our formation energy calculations with the screened HSE06 hybrid functional for exchange and correlation [69,70].We terminated our calculations when the forces in the atomic-position relaxations dropped below a threshold of 10 −2 eV•Å −1 .The wavefunctions were expanded in a planewave basis with a cutoff energy of 430 eV and the size of the supercell was 512 atoms (4 × 4 × 4 multiple of the conventional unit cell).We employed Γ-point integration to evaluate energies.The elements used in our calculations and the associated ground-state structures and values of their chemical potentials are: N (β hexagonal close-packed structure, −11.39 eV/atom) and C (diamond structure, −11.28 eV/atom).Only ground-state formation energies were computed in this work.

A. Thermally driven charge transfer
As in earlier work [71], we consider the expected rate of transfer for an electron associated with a defect in a crystal, but no longer in the fully dilute limit.Suppose the electron has a definite momentum at each point in time hk(t).In the non-relativistic limit, we can compute the rate associated with the transfer of the electron with effective mass m * across a distance |∆r| between ionized defects A and D as Above, we have used the definition of the Dirac delta function to convey that the rate is obtained by considering the average projection of the velocity onto the desired displacement of the electron divided by the distance |∆r| at the time the electron has completed the trajectory.
The likelihood that the electron has enough initial kinetic energy to overcome the ionizing potential and recover the neutral system is given by a Fermi-Dirac distribution.We consider the donor level E F D to be associated with an ionized donor-dopant D and the acceptor level E F A to be associated with an ionized color center A. We will show below that the required initial kinetic energy K initial is equal to the difference between the adiabatic charge-transition levels E F D and E F A for defects D and A respectively, up to a correction of order k B T .We can therefore write the expected rate of charge transfer as where T is the temperature and k B is the Boltzmann constant.The ionization probability is modeled through the Fermi-Dirac distribution since we consider thermal processes for the charge recombination and not processes involving optical excitation.

B. Introducing electric forces
The calculation of the built-in electric field between defect pairs, which counteracts the electromagnetic potential that results in their ionization, requires the evaluation of electrostatic potentials.
In order to compute the total energies of the systems from which the electrostatic potentials will be obtained, we solve for the eigenvalues of the Hamiltonian [72] where r i and m e denote the position and rest mass of electron i, respectively, and R I , Z I , and M I are the position, valence charge, and rest mass of ion I, respectively.We apply the Born-Oppenheimer approximation to decouple the electronic and ionic degrees of freedom and solve for the electronic degrees of freedom given static ionic positions.If the eigenvalues for the electronic Hamiltonian have zero dispersion in reciprocal space, the expectation value of the velocity of the electrons will be zero implying that we can neglect the kinetic terms.Therefore, for defects introduced into a crystal in the dilute limit where the dispersion of the donor and acceptor levels is vanishingly small, our Hamiltonian effectively captures Coulomb or electrostatic potentials.67,71,[73][74][75][76][77][78][79] calculated in the dilute limit for defects capture electrostatic potentials as a result of the aforementioned arguments and would therefore allow for the determination of built-in fields.In ∆H f (X q , {µ X i }, E F ), X is the defect species for which the formation energy is being calculated, q is the charge state of X, {µ X i } is the set of chemical potentials for the constituents of X, and E F is the Fermi level.In order to consistently determine the potentials at the locations of the defect species and to be in line with standard conventions, we use the neutral system in each case as the reference for the zero of the energy.

The formation energies ∆H
The potential due to each charged defect consequently becomes the difference between the energies of the charged-defect containing system and the neutral system divided by the compensating background charge.The potential associated with a defect X with charge q is then simply and similarly for a defect Y with charge −q the potential follows from (4) upon the substitution X → Y and q → −q, with r X and r Y denoting the locations of the respective defects.Here, E F is set to zero since its inclusion in the expression for the formation energy subtracts out the energy associated with adding charge to the defect, which is no longer necessary if we are employing the electrostatic energy corresponding to a given charged defect.Above we have also treated the defects in the dilute limit so that the compensating background charge of the computational supercell can be treated as a point charge relative to the entire crystal.The built-in field at the location of the defect Y is then given by where ∆r = r Y − r X .
We demonstrate the equivalence with our earlier work [71] as follows.In our earlier work [71], we provided the expression where r A denoted the location of an acceptor A, r D denoted the location of a donor D, ∆r = r A − r D , and E V indicated the energy of the valence band extremum.We had found that Setting q = 1 implicitly in Eq. ( 5), we find agreement between Eq. ( 5) and Eq. ( 6) demonstrating the desired equivalence between the 'local' Fermi level and built-in electric-field formulations.

C. Computing the momentum due to electric forces
At every point in time the charge must have enough energy to sustain motion along the path between the defects.In order to compute the necessary momentum for an electron moving between the defects, we recognize that the momentum must be sufficient to overcome the ionizing electromagnetic potential energy.We note, however, that as the charge approaches defect Y from defect X a built-in field will emerge to cancel the ionizing potential.In order to capture the emergence of the built-in field, we employ an energy conservation argument.Effectively, our argument is that the sum of the potential energy leading to the ionization of the defect pair and the kinetic energy of the electron must be conserved, so that U + K = U(K = 0).Thus, one finds Above, the electromagnetic potential energy, U(K = 0), reflects the relative formation energies associated with the placement of the electron that will travel from defect X to defect Y.A more general potential would allow for the inclusion of arbitrary external fields.The value of the potential energy landscape is defined in the manner given in Eq. ( 9) at the location of defect X so that its value at the location of defect Y can be set to zero.
In order to explicitly determine the required momentum as a function of time, we apply the second law of motion and employ a discrete approximation given the small distances with the origin at the location of defect Y to obtain The non-relativistic limit has been applied above, which is justified by the fact that the maximum speed an electron can attain according to our calculations is less than 0.002c, where c is the speed of light.Once the electron arrives at defect Y from defect X, the kinetic energy would be dissipated in a process akin to the Mössbauer effect so that there is no need to produce a large initial change in momentum.Therefore, given an initial speed corresponding to a kinetic energy E thermal = k B T in a random direction at an angle θ with respect to the radial direction After writing hk = m * dr dt , the relation has been used to obtain Eq. (11).We neglect where we have used the definition of the donor and acceptor levels [71].
In previous work [71] we considered transport via extended states, resulting in capture cross sections that would have finite extent if the wavevectors of the charges were allowed to evolve under the influence of the corresponding electron or hole forces.Without wavevector evolution, the capture cross sections would be points.The wavevector evolution employed in this work for charges in extended states under the influence of forces is consistent with the experimental finding of a large hole-capture cross section for optically activated NV − centers [9][10][11].Inclusion of the dθ dt and d 2 θ dt 2 terms would allow for further refinement of the capture cross sections.Therefore, if the electron is located at defect Y with an initial kinetic energy E F D − E F A up to a correction of order k B T , then the electron will be able to return to defect X.For a nonzero gradient U(K = 0) > 0 and zero initial velocity, subsequent r will be less than r 0 = |∆r|, requiring a velocity given by applying the negative root in Eq. (11).We note that the required kinetic energy can equivalently be viewed as the energy needed to excite the excitonic system from its ground state corresponding to the ionized defect pair where the charges are bound to the defects to its excited stated corresponding to the defect pair where the charges are free from the defects and the defects are neutral, in which case the momentum simply follows from the electric force between the charges.
We can average the reciprocal of Eq. ( 2) over all possible initial velocity directions, given an initial thermal energy of E thermal = k B T , which yields ⟨τ⟩ ≈ 1 In more detail, the barrier at zero temperature (in the absence of phonons) is given by E F D − E F A .At finite temperature (in the presence of phonons), the barrier is shifted by −E initial = −k B T cos 2 (θ )•sgn(cos(θ )).Averaging this shift over all θ (from θ = 0 to θ = π) should give zero for the net shift of the barrier.The barrier is in the argument of an exponential, however, so that the contributions when sgn(cos(θ )) is negative carry more weight in the average than the contributions when sgn(cos(θ )) is positive.Thus, the effective barrier (the barrier we would obtain if we replaced the timescale expression that averages over θ with a timescale expression that employs a single fixed value for the barrier) is higher than the barrier obtained by neglecting the thermal correction.Indeed, at a population fraction of approximately 50% for NV − centers, the upper and lower estimates for the charge-state decay time computed without the thermal correction yield values that are approximately 3.98 times smaller than the values reported in Fig. 1.Therefore, the room temperature (300 K) contribution of phonons is to increase the barrier by approximately k B T ln(3.98)≈ 0.036 eV.
Care must be applied if θ > π/2 in which case we must first integrate with positive velocity from r = r 0 to r = r 0 exp(k B T • cos 2 (θ )/U(K = 0)) and then back from r = r 0 exp(k B T • cos 2 (θ )/U(K = 0)) to r = r 0 with negative velocity before performing the integral between r = r 0 and r = 0 with negative velocity.
In order to obtain the fraction of color centers that have undergone charge-state decay for a given timescale ⟨τ⟩, we compute the probability that a color center will have undergone chargedecay as Above, | • | denotes the volume of the enclosed region, B r 0 [0] is the closed ball of radius r 0 centered at the origin, l X = n where n X is the concentration of the ionizing dopant X, and d max is the maximum implantation depth, since we consider the case where d max < l X /2 in this work.The probability reflects the fact that since ⟨τ⟩ is monotonic in r 0 the fraction of color centers having undergone charge state decay for a given ⟨τ⟩ corresponding to a given r 0 will be given by the fraction of color centers that are separated from their ionizing dopant by a distance of r 0 or less.We evaluate Eqs. ( 12) and ( 13) for r 0 uniformly distributed between r 0 = a and r 0 = where a is the lattice constant of the conventional unit cell of diamond (a = 3.549 Å [71]).

IV. RESULTS AND DISCUSSION
A. Experimentally investigating charge-state decay of NV − The details of the experiment of Yuan et al. [9] investigating charge-state instability of nearsurface NV − centers in diamond are as follows.A green pulse was used to initialize NV centers in the negative charge state.The NV − centers were then left under darkness for a variable delay time following which they were read out using a charge-state-selective orange pulse.The samples used included a diamond sample labeled A that was implanted with an implantation dose of 5 × 10 8 cm −2 at an energy of 3 keV and that was polished, pre-etched, and 12 C enriched and a diamond sample labeled F that was implanted with an implantation dose of 1 × 10 9 cm −2 at an energy of 1.5 keV and that was polished and pre-etched.They found that charge-state decay for sample F occurred on a timescale from 11-300 ms, while for sample A much less decay was observed out to 1 s.The accelerated charge conversion for sample F was attributed to the availability of electron traps near the surface, in particular boron impurities.In the following, we therefore concern ourselves with the theoretical determination of the charge-state decay rate in sample A.

B. Effect of defect species and concentration on charge-conversion timescales
In order to determine timescales for charge transfer between ionized defect species, we apply Eqs. ( 12) and ( 13) for an ionized color center with acceptor level E F A transferring an electron to an ionized donor dopant with donor level E F D .The acceptor and donor levels serving as E F A and E F D , measured relative to the valence band maximum, are given by ε NV (0/−) ≈ 2.8 eV and ε N C (0/+) ≈ 3.6 eV, respectively.An effective mass of m * ≈ 1.48m e is obtained from fitting the bandstructure in Ref. [71].We account for the fact that the charge-state decay of individual NV − centers measured in the Yuan et al. experiment was to a steady-state relative population greater than 0.5 [9], where 0.5 is the value corresponding to local pinning of the Fermi level at ε NV (0/−), by introducing a shift in E F A such that the local Fermi level E F A would produce a relative NV − population that would correspond to the experimentally measured value.Explicitly, for a final relative population of p, the shift is ∆E F A = k B T ln (1−p 0 )p p 0 (1−p) , where p 0 = 0.5.This result is obtained by considering the relative populations of the charge states of a single NV defect.Since we are considering the relative populations for a single defect, we can drop the contribution from configurational entropy.Therefore, the relative population of the NV − state is given by where we have assumed that the relative populations of the charge states other than 0 and −1 are negligible since the Fermi level is pinned near the ε NV (0/−) charge-transition level and have used and write E F = ∆E F + ε NV (0/−), then it follows immediately that The maximum shift is ∆E F A ≈ 0.07 eV and the minimum shift is ∆E F A ≈ 0.03 eV.The corresponding values of p are 0.93 and 0.76, respectively [9].We also account for the effect of the surface, which we assume to have an ether-like termination [80], so that a donor level measured relative to the valence band maximum of 3.4 eV [80,81] is induced.For NV − centers near the surface the number of ionized surface donors should be commensurate with the number of ionized bulk donors [81], so we can employ an effective E F D given by averaging the donor-level values (E F D → (E F D +3.4 eV)/2).Such averaging would also bring our earlier results in better agreement with experimental results at shallow implantation depths [71,81].The result of these corrections produces good agreement with experiment (see Fig. 1).
Explicitly, the averaging of the bulk N C defect level, E F D , and the ether-like surface termination defect level of 3.4 eV follows from assuming that on average only one N C defect and one ether-like surface defect will equilibrate after each excitation of the sample needed to perform a measurement and before the NV itself has time to equilibrate with the dopant from which it obtained its charge.
The assumption is justified by that fact that in the model proposed in this work, which considers the limit that is not fully dilute (meaning that charge does not need to fully enter the band edges to travel between defects), as well as in the model from Ref. [71], where charge must enter the band edges to travel between defects, the Fermi-level equilibration timescale will be the approximately the same for the NV -N C defect pair as for the ether-like defect-N C defect pair.This result follows from the fact that in the case of Ref. [71] the energy that governs the rate of equilibration is the difference between the energy of the conduction band minimum and the energy of the N C defect level (the defect levels of the ether-like surface defect and the NV can be neglected since they are both lower in energy than the N C defect level and the defects therefore act as acceptors for which the Femi-Dirac distribution is approximately equal to 1).In the case of this work, since the energy of the N C defect level is higher than that of either the NV or the ether-like surface defect, the factor of the Fermi-Dirac distribution yields approximately 1 for the ionization process and the speeds at which the charge travels between the N C and the NV and between the N C and the ether-like surface defect differ by a factor on the order of unity since the speeds evolve as the square root of energies that differ by a factor of order unity.By contrast, if the measurement timescale is such that more than one ether-like surface defect is allowed to equilibrate with more than one N C defect, then the value of the Fermi level is obtained from more general charge conservation considerations (Eq. (36) in Ref. [74]).
FIG. 1. Fraction of NV − population remaining as a function of time, calculated via Eqs.( 12) and ( 13) with T = 300 K. We model charge as undergoing transfer from NV − to N + C .The ionizing dopant concentration is n X = 5×10 8 cm −2 d max , using an estimated maximum implantation depth of d max = 3.5 nm • keV −1 E imp where E imp is the implantation energy [81].We employ E imp = 3 keV [9].The light-blue shaded region represents the result of introducing a shift in E F A of 0.03 eV ≲ ∆E F A ≲ 0.07 eV and of applying the transformation E F D → (E F D + 3.4 eV)/2.The orange dashed line indicates experimental results for the timescale associated with charge-state decay of four representative individual NV − centers in sample A with final relative populations between p = 0.76 and p = 0.93 [9].Since the four individual NV − centers were representative, the fraction of the NV − population remaining after the charge-state decay of each center can lie anywhere between 0 and 1.

V. CONCLUSIONS
We have shown that the precise electronic structure of an ionizing-dopant species and of an ionized color center are highly relevant to the charge-state decay characteristics of the ionized color center.The concentrations of dopants and color centers are also integral to the elucidation of charge-transfer rates within a semiconductor sample.A key implication of our results is that, in order to mitigate charge-state decay for ionized color centers in semiconductors, the color center should be chosen such that its charge-transition level lies much lower in energy than the donor level of the ionizing dopant.

− 1 .
Above, ∆t = 0 r 0 dr dr dt −1 and E initial = k B T cos 2 (θ ) • sgn(cos(θ )).We account for phonon excitations by the introduction of the thermal correction E initial .The introduction of the correction of order k B T influences the barrier for charge transfer in a manner analogous to the effect of phonons on the absorption energy of a fluorescent defect, namely the averaged timescale ⟨τ⟩ behaves as if the defect pair were experiencing a larger barrier for charge transfer than without the thermal correction.
R.K.D. gratefully acknowledges financial support from the Princeton Presidential Postdoctoral Research Fellowship and from the National Academies of Science, Engineering, and Medicine Ford Foundation Postdoctoral Fellowship program.We additionally acknowledge support by the STC Center for Integrated Quantum Materials, NSF Grant No. DMR-1231319.Finally, we wish to acknowledge insightful feedback from Nathalie P. de Leon and fruitful discussions with Pengning