Methods to accelerate high-throughput screening of atomic qubit candidates in van der Waals materials

Research in the field of layered materials has grown quickly since the development of the high-quality sample yielding a scotch-tape exfoliation technique.^1,2 For instance, the quantum spin Hall effect up to 100 K was observed in monolayer tungsten ditelluride (WTe$2$⁠).³ Another important discovery was the observation of a correlated insulator state by tuning the twist degree of freedom in bilayer structures of graphene due to the presence of flat bands near zero Fermi energy at a twist angle of about $1.1°$⁠, which upon electrostatic doping yields superconducting states with a critical temperature of up to 1.7 K.^4,5 Research in the field of point defect qubit candidates has also grown rapidly, particularly since the detection of single negatively charged nitrogen–vacancy (N$V−$⁠) color centers in diamond.⁶ These qubit candidates consist of defects in the crystal structure involving substitutional or interstitial atoms and/or vacancies and act as single-photon sources as well as sources of electronic spin. Desirable characteristics for such qubits include indistinguishability of emitted photons, negligible spectral diffusion, and long spin coherence time. State-of-the-art results include second-long coherence time of N$V$ centers,⁷ km-scale entanglement of two N$V$s,⁸ and the discovery of spectrally stable germanium–vacancy (Ge$V$⁠) and silicon–vacancy (Si$V$⁠) point defects.^9,10 The properties of spin qubits can be tuned by changing the crystal structure and the nature of the defect. Engineering spin qubits with properties that are desirable for applications in sensing, quantum computing, and quantum communications are at the forefront of research interest. Nevertheless, identifying the right combination of crystal structures and defects is enormously challenging due to the combinatorially large number of potential candidates that need to be explored using experiments or first-principles calculations.

The layered material hexagonal boron nitride (⁠$h$-BN) is of particular interest among layered materials as it has a wide bandgap,¹¹ enabling its ability to encapsulate other layered materials and to host a variety of point defects, giving rise to optical transitions lying within the range of its bandgap.^2,12–15 The issue with $h$-BN as a host for a qubit candidate employing electronic spins is that the atomic nuclei of boron and nitrogen have spins as well, which can cause spin-decoherence of the electronic spin state. Such an argument invites consideration of SiS$2$ as a host, which theory has predicted to exist in a layered form between 2.8 and 3.5 GPa for the space group P$21$/c.¹⁶ The SiS$2$ host would be diamond-like or silicon carbide (SiC)-like in that the atoms that constitute it have a low natural abundance of isotopes with nuclear spin, which can be further suppressed via isotope purification in growth, diamond and SiC being systems that have been extensively studied as hosts for qubit candidates.^9,10,17–32 The point defect we will investigate in SiS$2$ has the further advantage of exhibiting inversion symmetry, which eliminates the issue of an electric dipole moment making it susceptible to external noise and local fields that cause broadening of transitions.

Given the growing number of experimental investigations, the need to rapidly identify promising point defect qubit candidates is being recognized. Indeed, recent work³³ has looked into automating the process of characterizing point defects through a code named ADAQ, to aid in the more rapid identification of defect signatures from experiments, and a push toward high-throughput point defect calculations was also present in the earlier PyCDT code.³⁴ In identifying point defects in experiments, knowledge of the zero-phonon line (ZPL) transition is essential. We present in this paper two methods that can be incorporated into high-throughput searches for point defect qubit candidates. The first is a quick new method for estimating the error associated with computing the energy of the zero-phonon line (ZPL) transition of a point defect using Janak’s theorem, and the second is a new method for rapid convergence of excited state self-consistent field (SCF) calculations. These new tools are an important step in accelerating the high-throughput discovery of materials that can support spin qubits for applications in sensing, quantum communication, and quantum computing.

Janak’s theorem states that³⁵ 

where $E$ is the total energy, $ni$ is the orbital occupation of the $i$th orbital, and $ϵi$ is the corresponding eigenvalue of the orbital. This theorem is similar to Koopmans’ theorem for the Hartree–Fock theory³⁶ in that it relates energy differences under a change in the number of electrons to orbital energies. The theorem is in fact a density-functional theory (DFT) version of a theorem originally proved by Slater for his $X$–$α$ method, which was introduced in an early attempt to account for both exchange and correlation in electronic structure calculations.^37,38 As Janak’s theorem does not imply that $ΔE=ϵiΔni$ for finite $Δni$⁠, an estimate of the error associated with using the theorem to determine excitation energies where orbital occupations change by integer amounts is necessary. Applying the theorem and the method of error estimation is inherently faster than computing the energy of the excited state for a given point defect, as the error calculation amounts to terminating the excited state calculation before full convergence. We also show that excited state calculations for systems with small strain converge faster if the charge density is initialized by appropriately mixing the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) from the ground state calculation, leaving the contribution to the charge density from the remaining orbitals unchanged, as compared to initialization of the charge density from a superposition of atomic charges.

In this work, we first outline the computational methods in Sec. II. In Sec. III, we include results related to (i) ZPL estimations for silicon monovacancies in $4H$-SiC (see Sec. III A) and (ii) the singly negatively charged calcium vacancy in SiS$2$ (see Sec. III B, which also includes stability calculations) and a discussion of these results, ending with a summary of our conclusions.

The purpose of this work is to extend the notion of a mixing parameter³⁹ to the initialization of the charge density and to bolster the validity of using Janak’s theorem³⁵ to calculate the zero-phonon line (ZPL) by demonstrating success for the case of the singly negatively charged calcium vacancy in SiS$2$⁠, which we investigate as a potential qubit candidate. Our approach builds upon previous work employing Janak’s theorem to calculate the electronic properties of excited states of atoms, molecules, and solids.^40,41 In this work, we introduce an innovation by employing Janak’s theorem for the calculation of defect properties and additionally providing a lowest order estimate of the error associated with using the theorem for an integral change in the occupation of the single-particle states. We argue that the calculations in excess of the ground state calculation in the $Δ$SCF method are not always necessary, which directly follows from Janak’s theorem in the limit where the change in occupation is infinitesimal. To motivate the argument, we compute the error between the $Δ$SCF method and the approach of using Janak’s theorem to the lowest order under a change in occupation.

We consider the operator from the single-particle equations in DFT,

where $n(r)$ is the particle number density, the first term represents the kinetic energy of noninteracting quasiparticles with electron mass, the second term represents the external potential, the third term represents the Hartree potential for the Coulomb interaction between the quasiparticles, and the last term represents the exchange-correlation potential recapturing the fermionic and many-body nature of electronic interactions.

The single-particle equation for the single-particle state $ϕi(n)$ with eigenvalue $ϵi(n)$ then reads

Under a change in occupation, let $n′(r)$ be the new number density and let $ϕi(n′)$ be the new $i$th single-particle state such that

is satisfied.

The equation for the total energy from DFT is⁴² 

If we let

and

then we can approximate the $Δ$SCF result as

Therefore, if we only take the difference in ground state eigenvalues, the associated error is

We can choose to obtain $n′(r)$ at any arbitrary iteration in the full constrained-occupation calculation for the excited state. Therefore, in the limit where $n′(r)$ is obtained at the final iteration of the full constrained-occupation calculation for the excited state, the error becomes exact. In performing the full constrained-occupation calculation for the excited state, we have also explored initializing the charge density by mixing the HOMO and LUMO by varying amounts.

To obtain the defect levels and total energies, we performed first-principles DFT calculations for the various defects in $4H$-SiC and SiS$2$ using the VASP code^43–45 and the QUANTUM ESPRESSO code^46,47 for $Δ$SCF calculations. In VASP, atomic structures were first converged using the generalized gradient approximation (GGA) for the exchange-correlation energy of electrons, as parameterized by Perdew, Burke, and Erzenhof (PBE)⁴⁸ and then, for the calculation of defect levels of the uncompressed structure, using the screened hybrid functional of Heyd, Scuseria, and Ernzerhof (HSE) with the original parameters (0.2 Å$−1$ for screening and 25% for mixing).^49,50 The terms in $F[n(r)]$ were obtained for the ground state and for the state with the changed occupation, where a non-SCF calculation was performed until convergence keeping the charge density fixed in the latter case. The eigenvalues, $ϵi(n′)$⁠, were also provided by the non-SCF calculation. Code from the work by Feenstra et al.⁵¹ was used to change orbital occupations. In QUANTUM ESPRESSO, we performed $Δ$SCF calculations to investigate the SiS$2$ system using PAW pseudopotentials⁴⁵ with a 108-atom supercell with gamma-point integration. Modified source code was used to alter the charge density in QUANTUM ESPRESSO. The different values of strain investigated in this work were for biaxial strain. For strain values other than the uncompressed structure, the lattice parameters along the strain directions were rescaled from the lattice parameters of the uncompressed structure. The out-of-plane lattice parameter was held fixed. Constant volume relaxations were then performed for the strained structures as described above. Additional computational details may be found in the  Appendix.

For first-principles phonon calculations, performed using Phonopy,⁵² the atomic positions in the stoichiometric conventional unit cell were relaxed until the magnitude of the Hellmann–Feynman forces was smaller than $10−6$ eV$Å−1$ with a Monkhorst–Pack grid of $6×6×2$ k-points. Supercells containing 108 atoms (⁠$3×3×1$ multiple of the conventional unit cell) with appropriately scaled k-point grids and a cutoff energy of 500 eV were then constructed from the stoichiometric conventional unit cell and used to obtain the force constants to compute the phonon dispersion.

The formation energies of the calcium vacancy in SiS$2$ in various charge states were calculated according to the formula^53,54

where $q$ denotes the charge state, $Edef(q)$ is the total energy for the defect supercell with charge state $q$⁠, $E0$ is the total energy for the stoichiometric neutral supercell, $μi$ is the chemical potential of atom $i$⁠, $ni$ is a positive (negative) integer representing the number of atoms added (removed) from the system relative to the stoichiometric cell, $EVBM$ is the absolute position of the valence band maximum, $EF$ is the position of the Fermi level with respect to the valence band maximum (generally treated as a parameter), and $Ecorr(q)$ is a correction term to account for the finite size of the supercell when performing calculations for charged defects.⁵⁵ This correction term does not simply treat the charged defect as a point charge but rather considers the extended charge distribution. The chemical potentials of all the reference elements used in our calculations are listed as follows as a function of their crystal structure and total energy per atom: Si (diamond structure, $−5.42$ eV/atom); S (the total energy of a gas phase S$8$ molecule was calculated and the sublimation enthalpy was then subtracted,^56,57$−4.20$ eV/atom); and Ca (face-centered cubic structure, $−2.00$ eV/atom). Consideration of Si-rich or S-rich preparation conditions was made similar to the previous work.³⁰ 

In order to leverage the vast body of existing knowledge in the field, we utilized electronic databases of two-dimensional materials known as the 2D Materials Encyclopedia⁵⁸ in a top-down filtering approach and the Computational 2D Materials Database (C2DB)^59,60 to find the most promising host materials. The entries from the databases were first selected based on the presence of centrosymmetry in their space groups, yielding 4822 potential candidates with duplicates discarded. Fulfillment of this criterion allows for the possibility of creating in the host a defect without a dipole moment that could interact with local fields. The following filters were then applied (the number of candidates at each filtration step, discarding duplicates, is listed in parentheses):

1. PBE bandgap above 2 eV (985). The optical transitions of the qubit candidate implanted in the host must not introduce interference from the electronic states of the host.²¹ The value is taken to be at least as large as the energy of the optical transitions of well-known point defects such as the N$V−$ center in diamond and singly negatively charged silicon monovacancies in $4H$-SiC. Given that the PBE⁴⁸ functional generally underestimates the value of the bandgap, such a value allows for the defect levels for transitions similar in energy to those of negatively charged silicon monovacancies in $4H$-SiC and the N$V−$ center in diamond to be separated from the band edges.

2. Exfoliation energy below 80 meV/atom (441). The material must be easy to exfoliate and, therefore, to fabricate. The value is taken to be commensurate with exfoliation energies of common layered materials such as MoS$2$ with space group P$6¯$m2.

3. Does not contain an atom with a nuclear spin (6). Decoherence caused by the interaction of nuclear spins with the electronic ones of the qubit candidate must be minimized. We set the cutoff at $<5%$ natural abundance of isotopes containing a nuclear spin using as reference $4H$-SiC, which is the less restrictive among the two commonly studied systems of diamond and $4H$-SiC with defects showing long coherence times at room temperature.

The filtering process yielded SiS$2$ with space group P$21$/c as the host candidate with the lowest decomposition energy. The other host candidates left at the final stage of filtration were CaO with space group P4/nmm, SO$2$ with space group P$21$/c, CO with space group Cmme, Si$2$O$3$ with space group P$21$/m, and S$4$O$9$ with space group P$3¯$⁠. To our knowledge, none of these compounds have been used as qubit hosts. We note that similar criteria have been used in the work of Ferrenti et al.,⁶¹ though in that work, the condition regarding the natural abundance of isotopes with a nuclear spin was relaxed to a cutoff of $<50%$ and bulk materials were considered, yielding more known candidate hosts. Given the possibility of employing isotope purification, we invite researchers wishing to follow up on our work to relax the cutoff if the above candidates do not prove readily amenable to synthesis.

Working from the literature on defects in diamond and $4H$-SiC, we investigated vacancies, germanium–vacancy complexes, and lead–vacancy complexes as possible qubit candidates in the SiS$2$ host. However, these all break inversion symmetry when they are relaxed with spin polarization. On the other hand, the singly negatively charged calcium–vacancy complex is able to preserve inversion symmetry upon relaxation with spin polarization due in part to the large size of the calcium atom.

We begin by demonstrating the accuracy of Janak’s theorem for silicon monovacancies in $4H$-SiC. This system is an example of a success of the theorem though for another prominent system, the N$V−$ in diamond, Janak’s theorem alone is not sufficient to produce good agreement between experiment and theory. The lack of agreement in that case can be explained by the significant difference between the ground and excited state wavefunctions⁶² at a location of a significant change in the external or ionic potential (at the N atom in that case), which, from the single-particle equations, will in general lead to larger energy changes under a change in occupation compared to differences elsewhere. We can motivate this argument explicitly using what is referred to as “effective-mass theory,”^38,63 where the potential introduced by the impurity is considered to be a small perturbation to the crystal potential. Let $Vdefect(r)$ be the change in the external or ionic potential of the crystal due to the presence of the defect and expand the defect wavefunction of interest, $ϕdefect(r)$⁠, in the basis of the eigenfunctions, $ψk(r)$⁠, of the perfect crystal,

The wavefunction will obey the single-particle equation

where $Hcrystal$ is the single-particle hamiltonian for the perfect crystal so that

Using the orthogonality of the crystal wavefunctions and the expansion for $ϕdefect(r)$⁠, we then have

Thus, for perturbations to the defect wavefunction that equally change the projection of the defect wavefunction onto some perfect crystal eigenfunction (i.e., that result in the same $βk$ for some $k$⁠), the perturbation that will, in general, cause the greatest change in energy is the one where the other coefficients are such that the defect wavefunction changes most at the location of the defect potential. As we will see in Sec. III B, the partial charge densities associated with the orbitals changing occupation can be a quick way to identify when Janak’s theorem is likely to fail.

For the lattice parameters of the stoichiometric hexagonal unit cell of $4H$-SiC using the HSE06 functional, we find $a=3.08$ Å and $c=10.04$ Å. These values are in good agreement with the experimental values⁶⁴ of $a=3.07$ Å and $c=10.05$ Å and other theoretical values⁶⁵ of $a=3.07$ Å and $c=10.05$ Å. The structure of the $VSi$ is shown in Fig. 1. We observe that the expected $C3v$ symmetry is weakly broken from the distances displayed in Table I.

From the work of Soykal et al.,⁶⁶ we know that for a spin-polarized system, the approach of considering holes is equivalent to the approach of considering electrons. Defect levels calculated using the HSE06 functional for the singly negatively charged silicon monovacancy $VSi−$ in $4H$-SiC with $S=3/2$ at the two inequivalent $h$ and $k$ sites are illustrated schematically in Fig. 2 and are provided in Tables II and III. We note that the zero of energy is arbitrary as it depends on the pseudopotential, but, as we consider differences in eigenvalues, the value of the zero is not important. From the work of Soykal et al.,⁶⁶ in the hole picture, the ground state manifold is composed of the states, $∥uexey+iu¯e¯xe¯y⟩/2$⁠, $∥uexey−iu¯e¯xe¯y⟩/2$⁠, $∥uexe¯y+ue¯xey+u¯exey⟩/3$⁠, $∥u¯e¯xey+u¯exe¯y+ue¯xe¯y⟩/3$⁠, while the excited state manifold in the hole picture is composed of the states, $∥vexey+iv¯e¯xe¯y⟩/2$⁠, $∥vexey−iv¯e¯xe¯y⟩/2$⁠, $∥vexe¯y+ve¯xey+v¯exey⟩/3$⁠, and $∥v¯e¯xey+v¯exe¯y+ve¯xe¯y⟩/3$⁠. Above, $u$ and $v$ are single-particle orbitals transforming as the $A1$ irreducible representation of the $C3v$ point group, while $ex$ and $ey$ transform as the $x$ and $y$ components of the two-dimensional $E$ irreducible representation of the $C3v$ point group. The overbar denotes the minority spin state. To calculate the ZPL, we take the lowest energy hole states, which from Tables II and III we see must be the $∥u¯e¯xey+u¯exe¯y+ue¯xe¯y⟩/3$ (excited) and $∥v¯e¯xey+v¯exe¯y+ve¯xe¯y⟩/3$ (ground) states. Explicitly, in the following, we denote the eigenvalue for the majority spin single-particle orbital $γ$ as $ϵγ$ and for the minority spin single-particle orbital $γ¯$ as $ϵγ¯$⁠. Then, the ZPL we obtain for the $k$ site $VSi−$ is $(2(ϵv¯−ϵu¯)/3+(ϵv−ϵu)/3)=1.34$ eV (the remaining terms cancel), while the ZPL we obtain for the $h$ site $VSi−$ is $(2(ϵv¯−ϵu¯)/3+(ϵv−ϵu)/3)=1.43$ eV (where again, the remaining terms cancel) in excellent agreement with the experimental values of 1.35 and 1.44 eV, respectively.⁶⁷ The slight underestimation of the ZPL values may, therefore, be due to the slightly smaller HSE06 bandgap (3.18 eV compared to about 3.2 eV for experiment⁶⁸). We note, however, that other theoretical calculations yield 1.44 eV for the $k$ site and 1.54 eV for the $h$ site using the HSE06 functional.⁶⁹ We believe that error compensation in taking the difference of many eigenvalues may be causing the greater accuracy of our approximation to the $Δ$SCF method using Janak’s theorem than the $Δ$SCF method itself, though later theoretical work shows better agreement with our work and with experiment.⁷⁰ 

We note that, as the many-body states from the work of Soykal et al.⁶⁶ have been defined as Slater determinants of products of single-particle orbitals, the only requirement for the correctness of the calculation outlined above is that these single-particle orbitals be orthogonal, which is consistent with the orbitals constructed in the work of Soykal et al.⁶⁶ Though after performing structural relaxation with the HSE06 functional the expected $C3v$ symmetry is weakly broken, we have used the assignment and energetic ordering of the orbitals from the work of Soykal et al.⁶⁶ as such an ordering and assignment is consistent with obtaining the lowest energy spin-conserving excitation of the $S=3/2$ states that our calculations produce. Ultimately, these results demonstrate that in cases where the wavefunctions for the states that will change occupation do not appreciably differ at locations of a significant change in the ionic potential and where the change in the ionic potential can be viewed as a small perturbation, there is no need to calculate excited state properties to obtain the ZPL as values derived entirely from the ground state calculation can provide excellent agreement with experiment. These conditions are satisfied for the single vacancy as none of the defect wavefunctions are appreciable at the location of the removal of the atom and a single vacancy is a perturbation of less than $1%$ of the integrated ionic potential for the supercell size.

Based on the results of the calculations for silicon monovacancies in $4H$-SiC outlined above, Janak’s theorem shows promise for rapid estimation of ZPL values. Indeed, we will see below that the accuracy of Janak’s theorem is maintained for the strain (⁠$ε$⁠) values with the two lowest total energies, namely, the uncompressed and the $ε=+2%$ structures, for the singly negatively charged calcium vacancy in SiS$2$ and that the error calculations we have outlined in Sec. II A provide a clear indication of when Janak’s theorem fails.

As alluded to above, we now turn to demonstrating the continued accuracy of Janak’s theorem for the $ε$ values with the two lowest total energies for the singly negatively charged calcium vacancy in SiS$2$⁠. We additionally show that the associated error consistently identifies the larger discrepancies between the results of the theorem and the results of $Δ$SCF calculations. We also introduce results showing faster convergence of excited state SCF calculations when the charge density is initialized by mixing the HOMO and LUMO compared to when it is initialized from a superposition of atomic charges. For the lattice parameters of the stoichiometric unit cell SiS$2$ with space group P$21$/c using the PBE functional, we find $a=5.93$ Å, $b=8.13$ Å, $α=90°$⁠, $β=102.57°$⁠, $γ=90°$ and a monolayer thickness of $3.65$ Å with a vacuum of $16.94$ Å. Using the HSE06 functional, we find $a=5.88$ Å, $b=8.04$ Å, $α=90°$⁠, $β=103.11°$⁠, $γ=90°$ and a monolayer thickness of $3.59$ Å with a vacuum of $16.67$ Å. The lattice parameters with strain are included in Table IV. Structures for the singly negatively charged calcium vacancy for five strain values are found in Fig. 3. The values in Table V show that compressive strain tends to destroy the $Ci$ symmetry present in the uncompressed, $ε=+2%$ and $ε=+5%$ structures, due to buckling of the system.

The approximation to the ZPL using Janak’s theorem and using the $Δ$SCF method for different in-plane $ε$ can be found in Fig. 4, where the corresponding eigenvalues can be found in Table VI. The difference between the LUMO and HOMO eigenvalues is used in the application of Janak’s theorem. We observe little variation of the approximation to the ZPL with tensile $ε$⁠, but more variation with compressive $ε$⁠, which can also distort the structure. Since, due to Poisson’s ratio, the application of pressure should lead to tensile in-plane $ε$⁠, we do not expect the ZPL value to change significantly from what we have predicted in experimentally realizable structures. The error calculations perform best for the uncompressed (error of $0.0128$ eV) and $ε=+2%$ (error of $0.0249$ eV) structures, which have the smallest energy differences, obtaining the correct sign of the error as well, but fail to capture the correct sign and are much too large for the remaining $ε$ values (errors of $−17.4096$ eV, $−0.3898$ eV, and $−0.6123$ eV for $ε=−5%,−2%$⁠, and $+5%$⁠, respectively). We note that for $ε=−5%$⁠, the Ca atom departed from the position that preserved inversion symmetry and that the strain value was, therefore, not similar to the other structures. The larger errors nonetheless correctly indicate more significant disagreement between the $Δ$SCF result and the result from using Janak’s theorem. We note that based on a statistical learning based prediction of the bulk modulus,^71,72 the structures that should have the two lowest total energies given the requirement of an applied pressure of 2.8–3.5 GPa are the uncompressed and $ε=+2%$ structures, which is verified by QUANTUM ESPRESSO total energy calculations. Naturally, the uncompressed structure has the lower energy of the two structures.

As shown in Fig. 5, these structures with the two lowest total energies also responded best to performing the excited state calculation by replacing the initialization of the charge density from a superposition of atomic charges with the charge density of the ground state calculation, where the HOMO contribution is modified according to

$0≤β≤1$⁠. We note that a ground state relaxation is needed for both the standard method or default initialization and the approach of initializing the charge density from a mixture of the HOMO and LUMO orbitals for the excited state SCF calculation. This prerequisite calculation, therefore, does not change the net increase or decrease in the number of iterations for one method over another. Given that there may be instances where a ground state relaxation is not performed, we report the number of iterations required to converge the ground state SCF calculation for completeness. For $ε=−5%,−2%,+2%,+5%$ and the uncompressed structure, the number of iterations required were 29, 35, 34, 31, and 35, respectively. The value of the bandgap for the uncompressed structure using the PBE functional is 3.55 eV, and the defect has total spin $S=1/2$⁠. The $Δ$SCF calculations and the differences in ground state eigenvalues were from QUANTUM ESPRESSO.

Based on the computational efficiency of our method for determining the error associated with using Janak’s theorem demonstrated in Fig. 6, we argue that the approach of using Janak’s theorem should be acceptable for screening large numbers of potential point defect single-photon emitter candidates for desired ZPL values. Our approach to the initialization of the charge density also shows promise. A key point is that the search would focus on point defect candidates with total spin $S=1/2$⁠; otherwise, more involved group theoretic arguments would be required to determine the correct many-body hole wavefunction as in the work of Soykal et al.⁶⁶ Other specific systems, such as triplet systems,⁷³ will also be the focus of future work. The bandgap for the uncompressed structure using the HSE06 functional is 4.75 eV, and we again find that the defect has total spin $S=1/2$⁠. Using QUANTUM ESPRESSO, the HSE06 calculation yields a $Δ$SCF result of 0.1705 eV and a difference in ground state eigenvalues of 0.2676 eV. The error calculated as described in Sec. II yields 0.3008 eV. For hosts where the accuracy of hybrid functionals has been demonstrated, clearly, their continued use would be recommended. However, one cannot generally state that hybrid functionals will be more accurate than semilocal ones. Indeed, Johari and Shenoy⁷⁴ find that their PBE bandgap for monolayer MoS$2$ underestimates the experimental bandgap of 1.8 eV by just 0.12 eV, while Kuc et al.⁷⁵ and Ataca and Ciraci⁷⁶ found that the PBE0 and HSE06 hybrid functionals overestimate this bandgap by approximately 1 and 0.45 eV, respectively. In our work, we have, therefore, included both the PBE and HSE06 ZPL transitions for the uncompressed structure to reduce the possibility of missing the transition, noting that the transition energy does not vary significantly with strain for low total energy strain configurations. We recommend that any application of our methods using a semilocal functional is also performed with a hybrid functional.

The DFT partial charge densities obtained for ground hole (LUMO) and excited hole (HOMO) states and differences (LUMO $−$ HOMO) are shown in Fig. 7. The differences show the closeness of certain LUMO–HOMO pairs. More importantly, Fig. 7 shows that Janak’s theorem succeeds when the differences are not appreciable at the location of the greatest change in the ionic potential of the perfect lattice (at the location of the Ca atom) and does worse otherwise. The structure of the defect is four missing atoms, two silicon and two sulfur, replaced by a single calcium atom such that the resulting point defect is inversion-symmetric, as confirmed by the partial charge densities in Fig. 7. This inversion symmetry is broken, however, for in-plane $ε=−5%$⁠. The lowest formation energy as a function of Fermi level is displayed in Fig. 8 for the uncompressed structure. The plot demonstrates that in sulfur-rich preparation conditions, the introduction of a calcium vacancy actually stabilizes the SiS$2$ structure with space group P$21$/c. Given the fact that the singly negatively charged calcium vacancy only exists for a very limited range of Fermi level values, it would be necessary to pin the Fermi level at the appropriate value by doping the system. Using a $43$Ca atom instead of a $40$Ca atom could also lead to coupling between the nuclear spin and the electronic spin to implement a long-lived quantum memory realized with the nuclear spin.⁷⁷ 

We now turn to addressing the possibility of non-radiative decay with phonon emission upon excitation of the singly negatively charged calcium vacancy. Unlike in diamond where optical phonon modes with energies of about 0.17 eV can be found,⁷⁸ the energies of the phonon modes in SiS$2$ with space group P$21$/c are all below 20 THz as shown in Fig. 9, corresponding to an energy of about 0.08 eV. This result indicates that non-radiative decay would require multi-phonon processes, which would be rare.

We note that unlike diamond that is routinely synthesized at high pressure but shows no dynamical instability at ambient condition, our calculations show (see Fig. 9) that SiS$2$ with space group P$21$/c does exhibit a dynamical instability at ambient condition and would, therefore, need to be maintained at the 2.8 GPa required for its synthesis to be stable. We would like to make the following points: Experiments at high pressure (⁠$>102$ GPa) are routinely performed. In particular, the discovery of room-temperature superconductors at high pressure (e.g., Snider et al.⁷⁹ at 267 GPa) will undoubtedly propel further technical advances in this area. Experiments of emitters at high pressure (e.g., NV centers in diamond) have also been performed (e.g., Lesik et al.,⁸⁰ Yip et al.,⁸¹ and Hsieh et al.⁸²), showing such experiments to be within current feasibility. While it is not ideal that the proposed candidate material requires high pressure, it does not remove one of the key elements of the paper—the tools for the rapid screening of defect qubits. In summary, the case of the singly negatively charged calcium vacancy in SiS$2$ confirms that using Janak’s theorem and the associated error estimation for ZPL calculations is quick and clearly signals when the theorem is applicable. We have also shown that the constrained-occupation calculation for the excited state can be completed in some cases with fewer SCF iterations by carefully initializing the charge density.

In conclusion, we have developed a novel computational framework for high-throughput screening of the ZPL of atom-like qubits with applications for materials discovery, quantum sensing, and quantum computing. We propose the use of Janak’s theorem for ZPL calculations and provide a quick new method for estimating the associated error. We also outline a new method for reducing the number of iterations required to converge excited state SCF calculations, which is compared to the default procedure for carrying out such calculations. Our ZPL results are consistent with the previous experimental work where available and suggest that considering only ground state eigenvalues in the hole approach is a computationally efficient way of calculating optical excitation energies of color centers for screening purposes. Our novel method for initializing the charge density to reduce the number of iterations required for excited state calculations also works well for systems with small strain. Finally, we propose the new singly negatively charged calcium vacancy in SiS$2$⁠, which has the advantage of being inversion-symmetric and of containing a low density of nuclear spins that would lead to poor coherence properties of spins hosted in the material. However, this new proposed system would need to be maintained at the 2.8 GPa required for its synthesis to ensure stability. We also note that our rapid screening methods can create training data for machine learning and serve in technologies for predicting defect structures with desirable properties.

We thank Efthimios Kaxiras of Harvard University and Georgios A. Tritsaris for useful discussions. R.K.D. gratefully acknowledges financial support from the Princeton Presidential Postdoctoral Research Fellowship. We also acknowledge support by the STC Center for Integrated Quantum Materials (NSF Grant No. DMR-1231319). This work used computational resources of the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation (NSF) (Grant No. ACI-1548562),⁸³ on Stampede2 at TACC through allocation (No. TG-DMR120073), and of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231.

The data from DFT calculations that support the findings of this study and the modified source code for our charge density implementation are available from the corresponding author upon reasonable request.

We provide in this appendix details of the computations. The atomic positions were relaxed until the magnitude of the Hellmann–Feynman forces was smaller than $10−4$ eV$$Å$−1$ on each atom without spin polarization and subsequently until the magnitude of the Hellmann–Feynman forces was smaller than $10−2$ eV$$Å$−1$ on each atom with spin polarization to obtain defect levels, and, for the stoichiometric conventional unit cell, the lattice parameters were concurrently relaxed. The wavefunctions were expanded on a plane wave basis with a cutoff energy of 500 eV for all systems, and a Monkhorst–Pack grid of $6×6×2$ k-points was used for integrations in the reciprocal space for SiS$2$ and a gamma centered grid of $4×4×2$ k-points was used for integrations in the reciprocal space of $4H$-SiC. The relaxed lattice parameters of the stoichiometric unit cell were then used for the supercell structures. Formation energies and defect levels were calculated using a supercell of 108 atoms for SiS$2$ (⁠$3×3×1$ multiple of the stoichiometric unit cell) with appropriately scaled k-point grids. For $4H$-SiC, a supercell of 576 atoms (⁠$6×6×2$ multiple of the stoichiometric unit cell) was used.