Carbon defect qubit in two-dimensional WS2

Identifying and fabricating defect qubits in two-dimensional semiconductors are of great interest in exploring candidates for quantum information and sensing applications. A milestone has been recently achieved by demonstrating that single defect, a carbon atom substituting sulphur atom in single layer tungsten disulphide, can be engineered on demand at atomic size level precision, which holds a promise for a scalable and addressable unit. It is an immediate quest to reveal its potential as a qubit. To this end, we determine its electronic structure and optical properties from first principles. We identify the fingerprint of the neutral charge state of the defect in the scanning tunnelling spectrum. In the neutral defect, the giant spin-orbit coupling mixes the singlet and triplet excited states with resulting in phosphorescence at the telecom band that can be used to read out the spin state, and coherent driving with microwave excitation is also viable. Our results establish a scalable qubit in a two-dimensional material with spin-photon interface at the telecom wavelength region.


SUPPLEMENTARY NOTE 1: ASSESSMENT OF DENSITY FUNCTIONAL THEORY
The two-dimensional nature of layered transition metal dichalgonides (TMDCs) have large exciton binding energy. Experiment measures the exciton binding energy is about 0.37 eV [1] for WSe 2 and 0.32-0.71 eV in WS 2 [2,3]. To accurately consider the exciton effect and calculate the fundamental band gap, GW approximation with Bethe-Salpeter equation (BSE) should be used. Previous studies indicate that the quasiparticle band gap of WS 2 is about 2.8 − 2.9 eV and the exciton binding energy is about 0.6 eV [4][5][6]. However, GW calculation is extremely time consuming and it is usually out of reach for defective supercell systems. To reach a compromise, in this paper we slightly modify the parameters of DFT HSE06 functional through mixing parameter α and the range separation parameter µ to yield reasonable band gap [7,8], as shown in Supplementary Figure 1. With α = 0.40 and µ = 0.10 A −1 , the calculated bandgap is 2.76 eV with SOC included which is close to the experimental value. We call this functional HSE throughout the paper.
The SOC effect for the pristine WS 2 is evaluated in Supplementary Figure 2. Without SOC included, the calculated HSE bandgap is 3.13 eV with the parameters above. The SOC splits the double degenerate VBM and the gap is reduced by 0.37 eV. We find that the change of VBM contributes the most to the change in the value of band gap upon turning SOC. The VBM shifts upward about 0.32 eV while the CBM shifts downward about 0.05 eV.
Supplementary Figure 1. HSE SOC calculation for the band gap as a function of the parameters of HSE functionals. The experimental optical band gap is about 1.85 − 1.99 eV [9,10], however, to reproduce the fundamental band gap (∼ 2.7 eV) due to the large exciton energy, the mixing and screening parameters need to be modified.
Supplementary Figure 2. The calculated HSE band gap of primitive WS2 (a) without and (b) with SOC included. The black dash line indicates the energy referenced to the VBM without SOC included. We aligned the energy level in (b) to see the relative shift caused by SOC. The SOC effect shrinks the gap to 2.76 eV and splits the VBM by about 0.61 eV as denoted by red circle. SUPPLEMENTARY NOTE 2: SUBSTRATE EFFECT During experiment, the Fermi level of WS 2 is aligned by the graphene substrate so it is crucial to determine the influence of substrate on the electronic structure of WS 2 . We use 5 × 5 supercell of graphene and 4 × 4 supercell of WS 2 . The graphene layer is stratched to match WS 2 , corresponding to 1% strain on graphene which is tolerable. The band structures with and without graphene substrate are shown in Supplementary Figure 3. The defect bands lie at the same position with and without the graphene substrate included.
Beside the graphene substrate, we also consider hexagonal boron nitride (hBN) as substrate and capping layer since it could protect the C S defect and avoid hydrogen reattachment. Similar to graphene substrate, we use 5 × 5 supercell of hBN layers on top and bottom and 4 × 4 supercell of WS 2 as a middle layer. With the hBN included, the heterostructure exhibits type-I band alignment, as shown in Supplementary Figure 4. The defect level of the neutral charge state only shifts by 0.06 eV. With larger models, the energy difference can reduce to 0.01 eV. We demonstrate this by applying 8 × 8 supercell of WS 2 and 8 × 8 supercell of hBN sandwich layer structure. In this case, the geometry optimization was achieved within the affordable PBE functional and the electronic structure was calculated by HSE at that geometry as a single self-consistent calculation. We conclude that the weak van der Waals interaction between the layers minutely influences the electronic structure of WS 2 , therefore hBN would be an excellent protective layer.
Supplementary Figure 3. (a) Graphene substrate (below) and WS2 monolyer (top) heterostructure. The simulated PBE band structure of (b) freestanding WS2 layer with C − S is almost unchanged compared with (c) the bilayer model. The black and red lines indicate the spin up and spin down channels. Two flat bands represent the defect bands that arise due to the relatively small supercell and interaction of defect wave functions with their periodic images. The VBM of WS2 is referenced to zero.

SUPPLEMENTARY NOTE 3: CHARGE CORRECTION
It is well known that, in DFT simulation, the charged system with finite size supercell technique suffers from artificial electrostatic interaction between the periodic defect images and background charge. Therefore, correction scheme is needed to remove this interaction from both the total energy calculation and Kohn-Sham (KS) defect levels. The relative stability of charged defect could be evaluated by the formation energy. Here for the C S in WS 2 it is given by where E q d is the total energy of the system with defect at q charge state and E perfect is the total energy of the pristine system without defect. The µ c and µ s are the chemical potential of carbon and sulfur atoms, respectively, and can be derived from their pure bulk structure. The Fermi level Fermi represents the chemical potential of electron reservoir and should be aligned to the VBM energy of perfect WS 2 , perfect VBM . The E corr (q) is the correction term for the charged system due to the existence of electrostatic interactions with periodic boundary condition. This could be accomplished through supercell size scaling method and Freysoldt-Neugebauer-Van der Walle (FNV) correction as mentioned in the main text. [11] The FNV corrected formation energy demonstrates excellent agreement with the extrapolated value at infinite supercell size.
As a benchmark, we first perform the PBE calculation on various supercell sizes and all the atoms are fully relaxed. The models considered here range from 4 × 4 to 11 × 11. We fit the formation energy with a polynomial form as where  Table 1. We also did correction based on recently proposed self-consistent potential correction (SCPC) method, the correction energy is 0.42 eV and confirms the validation of our correction. Without charge correction, the HSE calculated occupied state is 0.51 eV above VBM while the unoccupied state is above the conduction band minimum (CBM). Supplementary Figure 5b shows the correction for defect KS levels with supercell size scaling method. For the negative charge state, both the occupied and unoccupied a 1 defect levels shift upward. The KS eigenvalue correction can be calculated with [12,13]: According to data in Supplementary Table 1, the correction energy is roughly 0.9 eV for both PBE and HSE functionals. With PBE functionals, the occupied state locates at 1.608 eV above VBM while the unoccupied state locates at 2.245 eV above VBM. Clearly the unoccupied state is already elevated above the CBM, which turns out that the observed unoccupied state does not belong to C − S . This is consistent with HSE result that the occupied state shifts to 1.406 eV above VBM while the unoccupied state resides at 4 eV above VBM. For the degenerate e state, the shift of defect levels is negligible.   Figure 7 that the second excited state (EX2) is much more complicated than EX1 is. Despite there is no occupied defect level in gap in the ground state electronic configuration, the unoccupied e level appears in the gap upon filling the state by photoexcitation which constitutes of a defect-to-defect transition. This could be due to the strong interaction between the defect states since both the e and a 1 states contain large d orbital components of tungsten atoms.

SUPPLEMENTARY NOTE 5: SOLUTION FOR THE JTE AND SOC
Previous studies on nitrogen-vacancy (NV) center in diamond pointed out that the 3 E excited state is unstable in C 3v symmetry due to the Jahn-Teller effect (JTE). [14][15][16] The JTE could reduce or even quench the SOC with damped factor p. This JTE is known as E⊗e dynamic JT system which was investigated by Ham and Bersuker. [17,18]  The symmetry breaking e phonons or local vibration modes reduce the high symmetry configuration that may couple the double degenerate E electron wave functions.
The SOC can be regarded as perturbation term if the JT energy is orders magnitude higher than SOC effect. Then λ Ham = | Ψ ± |Ĥ SOC | Ψ ± | = pλ. | Ψ ± is the vibronic wavefuntions caused by JT distortion. However, our calculated JT energy E JT from C 3v to C s is 89.9 meV which is very close to SOC and neither of them could be treated as perturbation. Then the SOC and electron-phonon interaction couple with each other and should be simultaneously solved.
We first discuss the system from pure electronic point of view. A is the usual anti-symmetrization operator: A|xy = 1 √ 2 (|xy − |yx ). So the four states in the main text could be expressed by: The single particle SOC operator is The two particle SOC operator is the sum of the two single particle SOC operators, A|e ↑ ± a ↓ and A|e ↓ ± a ↑ are eigenstates of H (2) Note that if there is no triplet-singlet splitting, |E xy and | 1 E would be degenerate. And if you switch on the SOC, it would lift the degeneracy, but the eigenstates would be still single determinant configurations.
So the total Hamiltonian isĤ Therefore, the 8 × 8 matrix that describes the electronic degrees of freedom is the following, In our present approximation, the Jahn-Teller effect originates in the real e x /e y character instead of the complex e ± = (e x ± ie y )/ √ 2 that is important towards the spin-orbit coupling. H whereσ z = |e x e x | − |e y e y |σ x = |e x e y | + |e y e x |.
We need to transform it to two-particle operator formalism aŝ JT .
Additionally, the Jahn-Teller interaction cannot connect the |E 12 , |A 12 with |E xy , | 1 E subspace. That is, all |E 12 , |A 12 's are maximal in spin, ↑↑ or ↓↓ (m S = ±1 ), but the singlets and m S = 0 states constitutes of antiparallel spins, ↑↓ or ↓↑. And since the Jahn-Teller operator cannot flip the spin there are no offdiagonal matrix elements appear between them. Also, the spin-orbit does not induce mixing, thus the m S = ±1 and m S = 0 subspace can be described fully independently. The two-particle Pauli matrices for Jahn-Teller effect acting inside the |E 12 , |A 12 states areσ After considering vibronic coupling with the additional Jahn-Teller effect, the SOC spliting between |E 12 and |A 12 levels will be reduced to pλ.
ForĤ JT , we solve it using series expansion where c nm and d nm are quantum oscillator coefficients. p could be calculated from p = nm c 2 nm − d 2 nm which is 0.017 indicating a giant reduction or almost quenched SOC due to JTE. Ham also numerically proposed fitting function p = exp −1.974 × χ 0.761 , where χ = E JT / ω e . Here χ = 2.99 and should be close to Huang-Rhys (HR) factor S between C 3v to C s configuration (2.55). Therefore, the deduced p = 0.011 from this approximation is close to our calculation with using the DJT Hamiltonian.

S
Based on Franck-Condon approximation, the phonon vibration could be evaluated through the overlap of phonon mode between electronic ground state and excited state. Huang-Rhys (HR) approximation further simplifies the mechanism that the interaction is equivalent between the defect and the lattice during excitation meaning the phonons are the same in the electronic ground and excited states.
The calculated HR factors are listed in Supplementary Table 3. We predict the Debye-Waller factor as e −S . With PBE functionals, the calculated HR factor for the first and second excited states are 3.00 and 4.82 in C 3v symmetry, respectively. These values agree well with previous findings. [19] We notice that the HSE functionals yield larger HR factors which might be attributed to the phonon calculation using the PBE functionals in the procedure.
The calculated phonon sideband of the PL spectrum of the two optical transitions are plotted in Supplementary  Figure 10. The partial HR factors S k represents the average number of phonons that emitted during excitation for specified phonon mode as illustrated in Supplementary Figure 11. For excitation 1 (EX1) with C 3v symmetry, we could screen out dominant localized phonon modes at 23 meV and 77 meV contributing to the neighbour sulphur atom. The 23-meV mode tends to shift the first neighbouring tungsten atoms towards the carbon impurity and then stretch the perpendicular distance between the carbon and underneath sulfur atoms, as shown in Supplementary  Figure 12. The 77-meV mode is an out-of-plane vibration of carbon solely. We also observe there are two delocalized modes at 14 meV and 18 meV. This might be attributed to the involvement of the band edge states during excitation. In the second excitation, the 23-meV mode still dominate while the 77-meV mode disappears because the carbon atom is relaxed out-of-plane during this excitation. One more delocalized mode appears at 51 meV with PBE level but not observable at HSE level.

SUPPLEMENTARY NOTE 7: HYPERFINE CONSTANT CALCULATION
The calculated hyperfine constants are shown in Supplementary Table 4. The electron spin coherence time (T 2 ) is particlular important for quantum applications and generally would be effected by proximate nuclear spins around the defect in the host material. In Supplementary Figure 13 we depict the distribution of the hyperfine constants of nuclear spins around the central carbon atom. The coherence time in WS 2 is about 11 ms therefore the farther tungsten atoms colored with purple could be used to store quantum information which do not significantly shorten the coherence time of the defect electron spin.
Supplementary Table 4. The calculated hyperfine constant with HSE functionals with contribution of the spin polarization of the core states (A1c) in the Fermi-contact term. The first neighbouring W atoms and the S atom below CS defect are considered. The natural isotopic abundance of 183 W (I = 1/2), 33 S (I = 3/2), and 13 C (I = 1/2) are 14.31%, 0.76%, and 1.07% respectively. The subnotations of W indicate the first, second and third nearest atoms. The average of the principal values of the hyperfine tensor (Average) is defined as (Axx + Ayy + Azz)/3.

Atom
Axx ( Figure 13. The hyperfine constants of tungsten atoms. The purple ones with hyperfine constants around 10 MHz and the blue ones indicate hyperfine constants less than 1 MHz.

SUPPLEMENTARY NOTE 8: DISCUSSION ON THE NATURE OF THE EXCITED STATE OF THE NEUTRAL C DEFECT
We assumed that the lowest energy excited state of the neutral C defect can be described as promoting an electron from the e orbital to the empty a 1 orbital in the gap. We calculated the excited state and the energy by ∆SCF method. The choice of the selection of the occupied and unoccupied (virtual) orbitals is based on the insight of the possible excitation routes but, in principle, it might be too restrictive and several pairs of occupied and virtual orbitals might be incorporated to describe the exciton wavefunction as the occupied orbital resides resonant in the valence band in the ground state electronic configuration.
We carried out further ∆SCF and many-body perturbation theory GW+BSE calculations [20,21] within Γ-point approximation. We here use a smaller supercell, 4 × 4, in order to reach close-to-converged calculation for the GW method which is very computationally demanding in the VASP implementation as more than 1500 bands were included in the single-shot G 0 W 0 calculation. The energy cutoff for the response function is set to be 150 eV. The Tamm-Dancoff approximation was used to solve BSE. The highest seven valence bands and seven lowest conduction bands are considered as basis for the excitonic state in the BSE procedure. We found that the results do not change by including more valence bands in this procedure. The calculations are based on the optimized HSE functional which resulted in the optimized geometry and the electronic structure of the neutral carbon defect in WS 2 . Here, we do not consider spin-orbit interaction. The ∆SCF procedure is applied at fixed geometry as obtained in the ground state, so vertical excitation energy is calculated. Although, this procedure is not entirely converged because of the small supercell but provides qualitatively good BSE results.
We find that the empty a1 defect level in the band gap by HSE agrees within 0.07 eV with G 0 W 0 quasi-particle level with respect to the calculated valence band maximum (see Supplementary Figure 14a). The ∆SCF method yields 1.43 eV vertical excitation energy in this supercell. The hole wavefunction is localized on the e x,y orbitals and the hole level pops up in the band gap, similarly to that in the large 6 × 6 supercell calculation. In the BSE absorption spectrum the first peak has a significant oscillator strength, and the weight of e x,y orbitals in the exciton state is about 86%. The dominant contribution of these orbitals in building up the hole part of the exciton justifies our procedure in ∆SCF. On the other hand, the first peak in the BSE absorption spectrum occurs at 1.77 eV (see Supplementary Figure 14d). The source of difference in the calculated ∆SCF and BSE excitation energies can be either the possibly missing 14% electron-hole pairs in building up the exciton in the ∆SCF procedure (with reducing the weight of the e x,y orbitals to 86%) or the missing back action on the quasi-particle wavefunction of the BSE method within Tamm-Dancoff approach as explained below.
According to our previous findings in other materials with similar electronic structure (defects in diamond, silicon carbide and silicon), the localization of the hole is highly critical which would decrease the excitation energy. This is analyzed by the inverse participation ratio method for the e and a 1 orbitals. The corresponding orbitals are plotted in Supplementary Figure 14b. One can see that the empty e x orbital in the ∆SCF orbital is more localized than the e x orbital in the ground state electronic configuration in the 4×4 calculation, and the empty e x level pops up in the gap in the ∆SCF procedure. We find that the e x orbital significantly localized when emptied whereas the a 1 orbital becomes less localized when occupied. As a consequence, the overlap between the e x and a 1 electron densities increases by ∆SCF with respect to that in the ground state DFT. This should result in a stronger attractive interaction between these particles with lowering its total energy. Obviously, the localization of hole part of the exciton wavefunction cannot occur by BSE method within Tamm-Dancoff approach (see Fig. 14b). We note that the localization of the e x orbital changes by +70% whereas it changes for a 1 orbital by −17% in the 6 × 6 supercell.
In summary, we conclude that the BSE method supports the dominant contribution of the selected orbitals in the ∆SCF method and the significant oscillator strength (optically allowed transition). The non-radiative transition rate between two electronic configurations under perturbation can be calculated with Here, W if is the electronic term and X if (T ) is temperature dependent phonon term. g is the equivalent energydegenerate atomic configurations and p in is the thermal population of the initial state (i), and f corresponds to the final state. The phonon matrix |χ i |Q − Q 0 |χ f | sums up the harmonic oscillator wave functions that enter the nonradiative recombination process. ∆E if is the calculated ZPL energy as described in the main text and ψ is the KS DFT wave function. Here two non-radiative processes are considered, one from 1 E singlet to 3 E triplet and the other one is from 3 E triplet to the singlet ground state, and the corrresponding W if are 0.017 and 0.001 eVamu −1/2Å −1 , respectively. We note that direct non-radiative process is possible (without any spin-flip) between 3 E and the singlet states because of the 2.1% singlet character in 3 E state. The dominant triplet character in 3 E may interact nonradiatively with the singlet states via spin-flip called intersystem crossing (ISC). The calculation method of the ISC rate is given in the main text. Then the ISC and non-radiative rates are calculated and shown in Supplementary  Figure 15. The 1 E singlet to 3 E triplet transition is very fast due to the small energy splitting and HR factor (0.12).
The total transition rate from singlet to triplet is 1 τ = 1 p 1 τ non + 1 p 3 τ ISC . p 1 and p 3 are the fraction of singlet and triplet characters in the 3 E triplet state (triplet character is dominant, p 1 p 3 ), respectively. The calculated lifetime of ISC is 0.68 ps and the direct non-radiative lifetime is about 38.2 ps (here, g = 3). The averaged lifetime is close to the ISC one going from 1 E to 3 E since it dominates the decay process. The non-zero spin sublevels 3 A 1 could be linked with singlets ground state with the non-axial part λ x,y of SOC. We find that either ISC or direct non-radiative decay process between the triplet 3 E and 1 A 1 ground state is relatively slow at low temperatures because of the large HR factor which greatly suppresses the |χ i |Q − Q 0 |χ f | term. Thus we conclude that the radiative decay dominates for 3 E state at low temperatures and the lifetime is still sufficiently long to carry out quantum information operations. At elevated temperatures, the ISC rate starts to dominate over the radiative rate (about 6.8 MHz, see main text) which may be disadvantageous for qubit operations.
Supplementary Figure 15. The intersystem crossing (ISC) and non-radiative (NRD) rates as a function of temperature for 1 E singlet to 3 E triplet (singlet to trplet) and 3 E triplet to 1 A1 ground state (triplet to ground state). The values at zero kelvin temperature are explicitly given.