Spin-defect qubits in two-dimensional transition metal dichalcogenides operating at telecom wavelengths

Abstract

Solid state quantum defects are promising candidates for scalable quantum information systems which can be seamlessly integrated with the conventional semiconductor electronic devices within the 3D monolithically integrated hybrid classical-quantum devices. Diamond nitrogen-vacancy (NV) center defects are the representative examples, but the controlled positioning of an NV center within bulk diamond is an outstanding challenge. Furthermore, quantum defect properties may not be easily tuned for bulk crystalline quantum defects. In comparison, 2D semiconductors, such as transition metal dichalcogenides (TMDs), are promising solid platform to host a quantum defect with tunable properties and a possibility of position control. Here, we computationally discover a promising defect family for spin qubit realization in 2D TMDs. The defects consist of transition metal atoms substituted at chalcogen sites with desirable spin-triplet ground state, zero-field splitting in the tens of GHz, and strong zero-phonon coupling to optical transitions in the highly desirable telecom band.

Introduction

An isolated point defect in a crystalline solid can be regarded as an artificial atom whose properties stem from the host material and bonding environment^1,2,3. The experimental demonstration of defects exhibiting long spin coherence times (T₂) and spin-selective optical transitions have made crystalline point defects one of the most promising platforms for the realization of long-distance quantum networks^1,4. However, finding a single-point defect that exhibits all the desirable traits for quantum entanglement network generation remains elusive. For example, the popular nitrogen-vacancy (NV) center (N_CV_C⁻¹ defect complex in diamond) does not operate at telecom wavelengths for low-loss fiber transmission (optical fiber telecom band: λ = 1260–1675 nm, or hν = 0.74–0.98 eV). On the other hand, Er-based qubits do but exhibit small optical oscillator strengths^5,6. All of the most promising point defects occur in three-dimensional (3D) bulk crystalline materials (diamond^7,8,9, SiC¹⁰, and oxides⁶), in which heterostructure fabrication, doping, and device fabrication remain challenging.

Here, we report on a family of point defects in 2D materials that combine moderate optical oscillator strengths, telecom operation, and low nuclear spin noise. Relative to 3D hosts, 2D hosts provide multiple advantages, including heterostructure engineering^11,12, reduced sensitivity to nuclear spin environment¹³, and ease of integration with photonic platforms. Furthermore, the placement of defects in a 2D layer (versus one buried in 3D) could be precisely controlled using a scanning tunneling microscope (STM)^14,15,16 or focused electron beam lithography¹⁷. In this regard, defects in monolayer hBN were theoretically investigated as qubit candidates in a 2D host^18,19,20. Long, ms-scale, longitudinal spin relaxation times have been demonstrated with a defect ensemble in hBN²¹, and optically detected magnetic resonance of single defects in hBN has been reported²². However, spin coherence times (also called quantum memory times) are limited to microseconds in hBN due to the nuclear spin environment (nuclear spins of all B and N isotopes). Ye et al. predicted that a nuclear-spin limited quantum memory time in MoS₂ can exceed milliseconds, even considering the natural abundance of nuclear spins before isotopic purification¹³. Moreover, the feasibility of isotopic purification of transition metal dichalcogenides (TMDs) further suppresses decoherence. Despite these promising properties, spin defect qubits in 2D TMDs remain uncharted territory even while defect-based single-photon emitters have been proposed and demonstrated^23,24,25,26.

The first-principles calculations based on density functional theory (DFT)^27,28 have extensively contributed to the characterization and identification of defect qubits in a wide range of solid hosts^{18,19,20,29,30,31,32,33}. In this work, we computationally search through defects for a spin defect qubit in 2D monolayer TMDs by means of hybrid DFT^34,35, known to be a quantitatively accurate method for solid-state defect calculations^29,36. As a result of the comprehensive characterization of fundamental qubit properties—electronic, magnetic, vibrational, optical properties, and thermodynamic stability—we report on a defect family of M_X in monolayer TMDs which turns out to be a promising candidate for quantum network applications.

Results

M_X defect family

To computationally discover spin defect qubits realized in monolayer TMDs, we first search through intrinsic (native) and dopant defects in H-MoS, where H is the notation for semiconducting monolayer TMDs³⁷. In this initial screening, we consider two criteria: (i) the spin-triplet ground state analogous to the NV center in diamond; (ii) spin-conserving intradefect optical transition without ionization of the defect²⁹. High spin states are desirable to decouple the spin from the S = 1/2 paramagnetic background and to allow spin control at zero magnetic field². Spin-conserving optical transitions are required for spin-state readout. First, we conducted DFT calculations based on the Perdew–Burke–Ernzerhof (PBE) functional³⁸ to quickly explore densities of states for various defect states (Supplementary Fig. 1). As a result of the initial screening in terms of the spin-triplet ground state, we identify three spin-triplet ground-state defect types: negatively charged donor-vacancy complexes (F_SV_S⁻¹ and Re_MoV_S⁻¹), Mo substitution for two S (Mo_S2), and Mo substitution for S (Mo_S). The donor-vacancy complexes have the spin-triplet ground state, but their occupied energy levels are not far enough away from the conduction band minimum to avoid ionization of the defects during intradefect optical excitation. Although Mo_S2 meets the two criteria, the defect is made of Mo_I and two V_S, so Mo_S2 is less likely to form than Mo_S, which can result in imprecise defect positioning, suffering from a random diffusion process during annealing; furthermore, a sulfur vacancy of Mo_S2 locates at the bottom sulfur layer of H-MoS₂ does not allow the STM tip manipulation. Out of the initial set of the spin-triplet donor-vacancy, substitution-type defects, we found Mo_S turns out to meet the screening criteria and was selected for a systematic study. We then further characterized the M_X defect family in the semiconducting H-MX₂ (M = Mo, W; X = S, Se, Te). WTe₂ is excluded because the most stable bulk phase of WTe₂ is the metallic T_d phase, not the semiconducting 2H phase³⁹, and thus unsuitable for hosting an optically active defect. In addition to the criteria above, it is practically desirable that the dopant M is different from the transition metal atoms constituting the host TMD so that we can optically distinguish the synthesized defect qubit from native anti-site defects and reach concentrations low enough for single qubit isolation.

Defect energy levels

We investigate the detailed electronic structures of selected defects in H-MX₂ using hybrid functional DFT calculations^34,35. The M_X defect family exhibits similar properties, and we focus the main text discussion on the W_Se defect in monolayer MoSe₂, which is found to have optical transitions in the telecom band along with W_S in MoS₂. The complete data for the family of defects investigated are listed in Table 1. We note that some defects in the family are not suitable spin qubit candidates. For example, W_Te in MoTe₂ does not have a spin-triplet ground state. For Mo_Se in WSe₂, the occupied defect levels calculated with spin–orbit coupling (SOC) are lower than the valance band maximum (VBM). Therefore, MoTe₂ and MoSe₂ can be excluded from the desirable host materials accommodating the M_X defect family. In addition to the M_X defects in monolayer TMDs, Table I includes our simulation results of the NV center in diamond²⁹ and the C_BV_N defect in monolayer hBN^18,19, which have been reported to meet the aforementioned criteria, although quantum chemistry approaches beyond the hybrid functional demonstrated that the ground state of the C_BV_N in hBN could be spin-singlet by taking into account multireference nature of the singlet state^26,40. These computational results are consistent with the previous reports for the NV center in diamond and C_BV_N in hBN, confirming that our simulation approaches are well-founded for defect qubit predictions. Comparison to these known centers also highlights the distinct features of the M_X defect family.

Table 1 Summary of calculated defect properties in diamond, hBN, and TMDs

The structural geometry and spin density of the N_CV_C⁻¹, C_BV_N, and W_Se defects are shown in Fig. 1a–c. All three defects possess the spin-triplet ground states with optical excitation pathways of spin-conserving intradefect transitions. The optical transitions lie within the bandgap E_g, prohibiting single-photon ionization of the defect [Fig. 1d–f]; since its estimation based on Kohn–Sham eigenvalues can be erroneous owing to the ambiguous interpretation of the Kohn–Sham eigenvalues, we further confirmed this from the comparison of the zero-phonon line energy and the ionization energy determined by the charge transition level, more precisely (e.g., the zero phonon line energy of W_Se in MoSe₂ is 0.79 eV, and the ionization energy of that is 1.2 eV). Similar to the NV center in diamond, the W_Se defect belongs to the C_3v point group, and the electron configuration of the 2D defects is identical to the hole configuration of the NV center. Two majority-spin electrons occupy doubly degenerate e_x and e_y orbitals, and the optical transition takes place between e_x,y and a₁ orbitals. The quantities between parentheses in Table 1 are given to estimate the SOC effects with heavy elements, where SOC reduces E_g and lifts the degeneracy of the e_x and e_y orbitals. More detailed calculations are required to determine SOC effects on spin coherence times, coherent spin-light interactions, and inter-system crossing² and will be addressed in future research work.

Fig. 1: Defect geometries and calculated electronic structures of (a,d) N_CV_C⁻¹ in diamond, (b, e) C_BV_N in hBN, and (c,f) W_Se in MoSe₂.

a–c Top (top) and side (bottom) views of defect geometries and spin densities of the defect qubits in the ground state (isosurface level = 0.003 Å⁻³). d–f Energy levels of the defect qubits. The green arrows indicate spin-conserving intradefect optical transition. Detailed physical quantities of possible combinations of M_X defects and MX₂ hosts are summarized in Table 1.

Defect formation energy

Defect formation energy is a crucial quantity to determine whether a proposed defect can be physically realized in a host solid. The defect formation energy of a defect X^q in a charge state q is given by^41,42,43

$${E}^{f}[{X}^{q}]={E}_{{{{{{\rm{tot}}}}}}}[{X}^{q}]+{E}_{{{{{{\rm{corr}}}}}}}^{q}-{E}_{{{{{{\rm{tot}}}}}}}[{{{{{\rm{pristine}}}}}}]-\mathop{\sum}\limits_{i}{n}_{i}{\mu }_{i}+q({{\epsilon }}_{F}+{{\epsilon }}_{{{{{{\rm{VBM}}}}}}}^{{{{{{\rm{pristine}}}}}}}-\Delta {V}_{0/{{{{{\rm{p}}}}}}})$$

(1)

where E_tot[X^q] and E_tot[pristine] are the total energies of a supercell with and without the defect X^q, respectively. n_i is the number of atoms of a species i added (positive) or removed (negative) from the pristine supercell, μ_i is the chemical potential of a species i. The chemical potential range was determined by considering competing phases (Supplementary Fig. 2) given in phase stability diagrams provided by Materials Project⁴⁴; based on the phase stability diagrams, we further computed the chemical potentials within the HSE06 hybrid functional to plot the formation energy diagrams at extreme conditions, such as the M-rich condition. The chemical potentials of C and N are obtained in the diamond crystal and the N₂ molecule, respectively. ${E}_{{{{{{\rm{corr}}}}}}}^{q}$ is the electrostatic correction, ${\epsilon }_{F}$ is the Fermi level, ${\epsilon }_{{{{{{\rm{VBM}}}}}}}^{{{{{{\rm{pristine}}}}}}}$ is the VBM energy level in the pristine supercell, and ΔV_0/p is the potential alignment term. The electrostatic correction is employed to take into account spurious image charge due to periodic cells and uniform background charge, where the Freysoldt–Neugebauer–Van de Walle (FNV) correction scheme^41,42,43,45 enables us to handle defects in anisotropic medium, such as 2D materials.

Figure 2 and Supplementary Fig. 4 show that the formation energy of an M_X defect is lower than the sum of formation energies of the two independent defects of M_I and V_X (i.e., E^f[M_x] < E^f[M_I] + E^f[V_x]), indicating that the formation of M_X defects is favorable. Compared with the NV center in diamond and C_BV_N in hBN, the formation energy of M_X in TMDs is small, so that the M_X defect family is expected to be readily created. Based on this formation energy, we can create the M_X defects by annealing a system with preexisting M_I and V_X defects. The formation of an antisite defect Mo_S in a MoS₂, which is among the M_X defect family, has been confirmed experimentally^46,47. Along with the experimental observation of Mo_S, the similar formation energy diagrams for the M_X defects in the family (Supplementary Fig. 4) support the feasible creation of the M_X defect family. Note that the dopant M needs to be different from transition metal atoms constituting the h ost TMD to distinguish the intentionally created defect. Since V_X is prevalent in TMDs⁴⁸, the M_X defect would be formed near the additional M_I after annealing. Supplementary Figure 6 shows defect formation energies of possible competing defects, where V_Se is much easier to be formed than V_Mo; thus, once we introduce W_I in the presence of abundant V_Se, the W_Se complex can be readily formed. The M atom could be incorporated via ion implantation or STM lithography.

Fig. 2: Defect formation energy diagrams.

Defect formation energies of the defect qubits for a N_CV_C⁻¹ in diamond, b C_BV_N in hBN in N-rich condition, and c W_Se in MoSe₂ in Mo-rich condition. The M-rich condition for MX₂ (N-rich condition for hBN) provides lower defect formation energies than the X-rich condition (B-rich condition) (Supplementary Fig. 3). c W_I + V_Se indicates the sum of formation energies of the two independent defects, W_I and V_Se. The orange-shaded area shows the range of stability of N_CV_C⁻¹, C_BV_N⁰, and W_Se⁰. Defect formation energy diagrams for other M_X in TMDs are displayed in Supplementary Fig. 4.

Zero-phonon line emission

Photon emission of defects plays a key role in qubit operation. Spin-conserving cycling transitions are utilized to read out the spin-qubit state. Zero-phonon-line (ZPL) transitions are utilized to realize spin-photon entanglement, which is required for generating spin-entangled quantum networks via photon measurement^49,50. The ZPL emission is also utilized as a spectroscopic fingerprint to identify the defect qubit^2,31,32. A photoluminescence line shape is composed of the ZPL and phonon sidebands. The contribution of the ZPL emission to the total emission is estimated by the Debye-Waller (DW) factor^31,32. Only the ZPL emission is useful for photon-spin entanglement schemes, and thus a high DW factor is desirable. The configuration coordinate diagram (adiabatic potential energy against configuration coordinate) is often utilized to investigate the ZPL emission (Fig. 3)^51,52. Here, the configuration coordinate displacement ΔQ is calculated as^51,52

$$\Delta Q=\sqrt{\mathop{\sum}\limits_{\alpha,i}{m}_{\alpha }{({R}_{e;\alpha i}-{R}_{g;\alpha i})}^{2}},$$

(2)

where i = {x, y, z}, m_α is the mass of atom α, ${R}_{g\left(e\right);{\alpha }i}$ is the equilibrium position in the ground (excited) state. The electronic excited state is calculated by using the constrained DFT⁵³. The number of phonons emitted during the optical transition can be quantified by the Huang–Rhys factor S. In the one-dimensional (1D) effective phonon approximation⁵¹, $S=\frac{\Delta E}{\hslash \omega }$, where ΔE and ℏω are described in Fig. 3a, and the effective phonon frequency ω is obtained using the harmonic oscillator approximation $E=\frac{1}{2}{\omega }^{2}{Q}^{2}$. Here, ΔE is the difference between the ZPL energy E_ZPL and the vertical emission energy. The DW factor is given by DW = e^{−S 31,32,52}. The ZPL energies, the Huang–Rhys factors, and the DW factors for the M_X defect family are shown in Table 1. The M_X defect family exhibits larger DW factors than the NV center in diamond and the C_BV_N in hBN except for in the MoTe₂ host (which is not promising from the energy level point of view, as discussed earlier). The large DW factors stem from the small curvature of the M_X defect family configuration coordinate diagram compared with the NV center in diamond and the C_BV_N in hBN (Fig. 3). In Table 1, the SOC-corrected ZPL energies between parentheses are approximated by estimating shifts in the defect energy levels shown in the same table. The ZPL energies of the M_X defect family typically lie around 1 eV, close or in the telecom band, with the SOC-corrected ZPL energies at 0.74 eV and 0.94 eV of the W_Se in MoSe₂ and W_S in MoS₂, respectively. As we show further below, the 2D host environment enables fine-tuning of the ZPL energy by applying strain to further minimize optical fiber transmission loss. While photoluminescence measurements in 2D TMDs have identified localized excitons and chalcogen vacancies^54,55,56, a ZPL that can be attributed to the M_X defect family has not been experimentally identified yet.

Fig. 3: Configuration coordinate diagrams responsible for ZPL emission.

Configuration coordinate diagrams of the defect qubits for a N_CV_C⁻¹ in diamond, b C_BV_N in hBN, and c W_Se in MoSe₂. The solid horizontal lines correspond to the phonon energy levels in the harmonic approximation. The ground state and the excited state are labeled as ³A₂ and ³E for N_CV_C⁻¹ in diamond³¹ and M_X in TMDs; the states are labeled as (1)³B₁ and (2)³B₁ for C_BV_N in hBN⁴⁰.

Zero-field splitting and hyperfine tensors

Magnetic properties of a spin qubit are of paramount importance for realizing quantum information applications^2,31,32. For instance, the zero-field splitting (ZFS) between the m_s = 0 and m_s = ±1 spin sublevels corresponds to the microwave energy to manipulate the qubit state at zero applied magnetic field and enables spin-selective resonant optical excitation. Furthermore, utilizing ZFS TMD-based defects may also be attractive for quantum sensing applications. In these applications, changes in the ZFS can be used to sense electric fields, strain, and temperature^57,58, while splitting of the m_s = ±1 state due to the Zeeman effect is used for magnetic sensing applications^59,60. Also important is the magnetic coupling of the spin-qubit to the crystal spin bath. A ZFS allows one to decouple the qubit spin from a paramagnetic electron spin-1/2 bath. In addition to a qubit spin state coupling to paramagnetic electron spins, there will also be hyperfine couplings to the crystal host nuclear bath¹³. In spin Hamiltonian, the ZFS and the hyperfine interaction are described as ${\sum}_{n}{\hat{{{{{{\bf{S}}}}}}}}^{{{{{{\bf{T}}}}}}}\cdot{{{{{{\bf{A}}}}}}}^{\left(n\right)}\cdot{\hat{{{{{{\bf{I}}}}}}}}^{(n)}$ and ${\hat{{{{{\bf{S}}}}}}}^{{{{{{\rm{T}}}}}}}\cdot{{{{{\bf{D}}}}}}\cdot{\hat{{{{{{\bf{S}}}}}}}}$, respectively, where $\hat{{{{{{\bf{S}}}}}}}$ is the electron spin, ${\hat{{{{{{\bf{I}}}}}}}}^{(n)}$ is the nuclear spin of nucleus n, D is the ZFS tensor, and A⁽ⁿ⁾ is the hyperfine tensor.

The ZFS tensor determines the dipolar spin–spin interaction between electrons and is given by^32,61

$${D}_{{ab}}=\frac{1}{2}\frac{{\mu }_{0}}{4\pi }\frac{{g}_{e}^{2}{\mu }_{B}^{2}}{S\left(2S-1\right)}\mathop{\sum }\limits_{i > j}^{{{{{{\rm{occupied}}}}}}}{\chi }_{{ij}}\left\langle {\Psi }_{{ij}}\left({{{{{{\bf{r}}}}}}}_{1},\,{{{{{{\bf{r}}}}}}}_{2}\right) \left|\frac{{r}^{2}{\delta }_{{ab}}-3{r}_{a}{r}_{b}}{{r}^{5}\,}\, \right|{\Psi }_{{ij}}\left({{{{{{\bf{r}}}}}}}_{1},\,{{{{{{\bf{r}}}}}}}_{2}\right)\right\rangle,$$

(3)

where μ₀ is the magnetic permeability of vacuum, g_e is the electron g-factor, μ_B is the Bohr magneton, χ_ij is +1 for parallel spins and −1 for antiparallel spins, r_a,b = (r₁ − r₂)_a,b, and |Ψ_ij(r₁, r₂)〉 is the slater determinant of ith and jth Kohn-Sham orbitals. After diagonalizing the ZFS tensor, one can obtain the ZFS value $D=\frac{3}{2}{D}_{{zz}}$, presented in Table 1 (see Supplementary Table 1 for the tensor elements, D_xx, D_yy, and D_zz). For the NV center in diamond, D calculated in this work is 2.86 GHz, which is close to the reported one^2,32. D of the M_X defect family are 10–20 GHz, about an order of magnitude larger than that of the NV center, which is within the experimentally accessible range of microwave control^62,63 and could enable higher-temperature resonant spin readout as well as the compatibility of higher Purcell factors⁶⁴ with resonant optical spin selectivity. The large D of the M_X defect family is attributed to stronger dipolar spin-spin interaction due to the more localized electron wavefunctions than the NV center. Note that because of the additional contribution of SOC⁶⁵, the ZFS could be even greater than the value presented in Table 1, especially with a heavy element, such as W.

The hyperfine tensor of nucleus n at r = 0 is calculated by using^32,66

$${A}_{{ab}}^{\left(n\right)}=\frac{{\mu }_{0}}{4\pi }\frac{{g}_{e}{\mu }_{B}{g}_{n}{\mu }_{n}}{S}\int {d}^{3}r{n}_{s}\left({{{{{\bf{r}}}}}}\right)\left[\left(\frac{8\pi }{3}\delta \left(r\right)\right)+\left(\frac{3{r}_{a}{r}_{b}}{{r}^{5}}-\frac{{\delta }_{{ab}}}{{r}^{3}}\right)\right],$$

(4)

where n_s(r) is the electron spin density [Fig. 1(a–c)], g_n is the nuclear g-factor⁶⁷, and μ_n is the nuclear magneton. In Eq. (4), the first parenthesis is the non-dipolar Fermi contact term, and the second parenthesis is the dipole–dipole interaction term. Table 2 displays the calculated and diagonalized hyperfine tensors of the NV center in diamond, C_BV_N in hBN, and W_Se in MoSe₂ (see Supplementary Table 2 for a full list of hyperfine tensors, including other M_X in the family) at the defect and nearest neighbor sites. The ¹⁸³W and ⁷⁷Se nuclear spins of the W_Se defect exhibits large hyperfine tensor elements, similar to the on-site interaction C_BV_N defect in hBN and the nearest neighbor ¹³C in the NV. Considering the number of equivalent sites, the total hyperfine coupling between the electron spin and nearby nuclear spins is not necessarily stronger than the NV center. Furthermore, the advantageous dimensionality¹³ and the isotopic purification for 2D TMDs are expected to provide an exceptionally coherent time, whereas 2D hBN is incapable of excluding spinful nuclear isotopes. One intriguing possibility with 2D TMDs is to completely engineer a nuclear spin quantum memory register⁶⁸ by STM lithography¹⁵. In this case, one would begin with an isotope-purified spin-0 host and incorporate a handful of nonzero spin nuclei in proximity to the defect.

Table 2 Calculated hyperfine tensors for N_CV_C⁻¹ in diamond, C_BV_N in hBN, and W_Se in MoSe₂

Radiative decay

In addition to the DW factor, the radiative recombination rate is an important optical property. For quantum information protocols, recombination rates should be fast enough to realize efficient spin initialization and readout^2,31,32. Practically, radiative rates should also exceed the rates of any nonradiative recombination processes. The radiative recombination rate, which is the inverse of the radiative recombination lifetime τ_R, is calculated using^32,69

$$\frac{1}{{\tau }_{{{{{{\rm{R}}}}}}}}=\frac{n{E}_{{ij}}^{3}{\left|{\mu }_{{ij}}\right|}^{2}}{3{\epsilon }_{0}\pi {c}^{3}{{{\hslash }}}^{4}},$$

(5)

where n is the refractive index, ϵ₀ is the vacuum permittivity, E_ij is the excitation energy that is substituted with E_ZPL, and ${\mu }_{{ij}}=\left\langle {\psi }_{j}|e{{{{{\bf{r}}}}}}\,|{\psi }_{i}\right\rangle$ is the transition dipole moment between the initial state |ψ_i〉 and the final state |ψ_j〉. Under the Frank-Condon approximation, we consider only the electronic component of the initial and final wavefunctions, which are occupied and empty Kohn-Sham orbitals of the spin-triplet ground state. Table 1 shows the calculated τ_R for the systems that we have examined so far. The W_Se in MoSe₂ exhibits a 4.2 µs decay time, which is four times shorter than the 20.5 µs decay time of the W_S in MoS₂. Overall, τ_R of the M_X defect family is 100–1000 times larger compared with the NV center in diamond and C_BV_N in hBN. In the M_X defect family, the optical transition between e_x,y, and a₁ is smaller due to the orbital selection rule (Laporte rule)⁷⁰ associated with distinct d orbital characters of their defect states, e_x, e_y, and a₁ (Supplementary Figure 6). While (slightly) shorter τ_R may be desirable, we note τ_R is already 5 orders of magnitude shorter than the current most promising defect telecom qubit, Er:³⁺Y₂SiO₅ where the intra-f-shell transitions are utilized, unlike the transition metal defects with d-orbital physics⁶. Moreover, for efficient photon collection, cavity integration is required, which can reduce τ_R by 4 orders of magnitude via the Purcell effect⁶. Due to the large ZFS, the system should still retain frequency-selective spin excitation for spin-photon entanglement and spin read-out even with the 4 orders of magnitude frequency broadening. TMDs can also provide multiple advantages in sensing. Due to the proximity to the surface, the exposed defect qubit on the surface of monolayer TMDs can compensate for the low radiative decay rate by suppressing internal reflection. Together with the radiative process, nonradiative recombination is a vital process determining quantum yield. The absence of crossing between the potential energy curves of ³E and ³A₂ shown in Fig. 3c indicates that the nonradiative transition between the triplet states is less likely to occur; however, further investigation is necessary to make sure the rare nonradiative transition because the transition could depend on many critical factors.

Intersystem crossing (ISC)

The transition between a triplet state and a singlet state can play an important role in a nonradiative process and can enable the low-fidelity room-temperature optical initialization and readout of the qubit-based sensors. The M_X defect family symmetrically resembles the antisite defect in monolayer TMDs and is expected to exhibit symmetry-allowed ISCs as in the antisite defect⁷¹. ISC is mediated by a combination of SOC and electron-phonon interaction. The crossing rate was calculated by the application of Fermi’s golden rule according to the formula^72,73:

$${\Gamma }_{{{{{{\rm{ISC}}}}}}}=4\pi {\lambda }_{\perp }^{2}{\widetilde{X}}_{{if}},$$

(6)

$${\widetilde{X}}_{{if}}=\mathop{\sum}\limits_{m}{w}_{m}\mathop{\sum}\limits_{n}{\left |\left\langle {\phi }_{{im}} \bigg|{\phi }_{{fn}}\right\rangle \right| }^{2}\delta \left(\Delta {E}_{{if}}+m{{\hslash }}{\omega }_{i}-n{{\hslash }}{\omega }_{f}\right),$$

(7)

where λ_⊥ is the transverse SOC constant between spin-singlet and spin-triplet states, ${\widetilde{X}}_{{if}}$ is the phonon wavefunction overlap between initial state i with phonon quantum number m and final state f with phonon quantum number n, ϕ_im and ϕ_fn are the phonon wavefunctions, ω_i and ω_f are the phonon frequencies, w_m is the occupation number of phonon according to Bose-Einstein distribution, and ΔE_if is the energy difference between the initial state and final state (See Methods for further details of phonon wavefunction overlap and SOC strength calculations). The ISC from the triplet excited states ³E to the singlet shelving state ¹A₁ can be symmetrically allowed when m_s = ±1⁷¹. The simulated transition rate of ISC from the triplet excited state to the singlet shelving state is 0.031 μs, which is shorter than the radiative lifetime 4.2 μs of the triplet excited state, which tells us that the proposed quantum defect can exhibit the initialization and readout operation via the spin-selective decay pathways (Fig. 4). We note, however, that for the high-fidelity initialization and readout required for computation and network, resonant, spin selective excitation is required along with avoided or minimized ISCs⁷⁴. Since SOC underlies the ISC transition⁷², we will be able to engineer ISC by utilizing various transition metal dopants with different SOCs.

Fig. 4: Sublevel structure of W_Se in MoSe₂.

The radiative processes are shown in the orange vertical line. The blue dashed lines show the symmetry-allowed ISC transitions from the triplet excited state ³E to the singlet state ¹A₁ and the transition from ¹A₁ to ³A₂, which are responsible for spin-selective decay, enabling the initialization and readout operations. The purple circular arrows within ZFS indicate the manipulation of qubit states by microwave.

Strain engineering

Strain can be effective in altering the dominant d orbital character by reducing the defect-crystal symmetry, which significantly modulates defect qubit properties, including the optical transition properties. Therefore, we can modify the radiative recombination rate under applied strain. As shown in Fig. 5, uniaxial strain along x or y breaks the C_3v symmetry and lifts the e_x and e_y degeneracy. Technically, the notation e_x,y is not valid when uniaxial strain is applied, and the notation is associated with their original orbital without strain. Although the biaxial strain does not break the C_3v symmetry, the biaxial strain affects orbital mixing, resulting in the modulation of τ_R and E_ZPL. If we consider the lowest excitation for qubit operation (e_x → a₁ for uniaxial strain along x, e_y → a₁ for uniaxial strain along y), uniaxial strain is always beneficial to achieve a shorter lifetime. The tensile biaxial strain would also be helpful. In addition to τ_R, the strain technique can be used to engineer the ZPL energy. As shown in Fig. 5d–f, strain shifts energy levels and changes the energy gaps between a₁ state and e_x,y states by a few hundred meV, which would provide a useful way to tune defects to a single operational frequency in a targeted communication band. 2D host materials are beneficial for the strain engineering of defect qubit properties because a single atomic sheet can accommodate a more significant mechanical strain (up to a few %) than bulk materials (typically less than 0.1%), and the strain will depend on the 2D material-substrate interfaces⁷⁵. Interestingly, δE pertaining to the second-lowest energy excitation, changes abruptly with small uniaxial strains as a consequence of the lifted degeneracy due to the symmetry breaking. The drastic response to external strain could be promising for highly susceptible quantum strain sensors.

Fig. 5: Strain effects on radiative recombination lifetime and ZPL energy.

a, b Uniaxial and c biaxial strain effects on radiative recombination lifetime τ_R. The inset in a shows the directions x and y. d, e Uniaxial and f biaxial strain effects on δE, modulation of the gap between corresponding eigenvalues (a₁ to e_x,y) in the ground state. The δE approximates the change of the ZPL energy, assuming the vertical shift of adiabatic potential energy curves occurs in the configuration coordinate diagram.

Discussion

We proposed the M_X defect family in monolayer TMDs as a promising solid-state defect qubit through systematic computational investigation of essential criteria: defect energy levels, defect formation energy, ZPL emission, ZFS, hyperfine tensor, radiative recombination rate, and ISC transition rate. Compared with the NV center in diamond and the C_BV_N defect in hBN, the proposed defects exhibited desirable qubit properties, operating at telecom wavelengths. Finally, we demonstrated strain effects on radiative recombination lifetime and defect energy levels, which provides a technique that we can exploit for further engineering qubit properties and applications to sensitive quantum strain sensors.

Among the various combinations of M and X, the W_Se defect in MoSe₂ and the W_S defect in MoS₂ are particularly promising candidates for quantum network applications with a ZPL transition in the telecom band. However, many of the family’s defects are promising candidates for the first demonstration of experimentally detected spin defect qubits in a 2D TMD host. Computationally, there is also further work to be performed. In particular, the role of spin-orbit coupling on spin T₁ lifetimes and coherent spin-light interactions should be investigated. The Debye temperatures of TMDs⁷⁶ are an order of magnitude smaller than those of the NV center in diamond and the C_BV_N defect in hBN; thus, it is reasonable to expect that spin relaxation time T₁ of the M_X defect family could be shorter than those of the counterparts due to a strong spin-phonon interaction^3,77. If T₂, such that T₂ ≤ 2T₁, is limited by T₁, one can explore different combinations of defects and hosts in the family to mitigate the spin-phonon interaction by reducing SOC. Other transition metal atoms in adjacent columns of the periodic table can also be explored to substitute for X along with a nonzero change state, implying expansive room for further exploration and qubit property engineering of 2D quantum defect systems. Having theoretically discovered and characterized the promising spin-defect qubits in monolayer TMDs, we opened a new door to the 2D world of research on spin-defect qubits.

Methods

First-principles calculations

We used Vienna Ab initio Simulation Package (VASP)^78,79 to perform the first-principles calculations based on density functional theory (DFT)^27,28. The Heyd–Scuseria–Ernzerhof (HSE06) hybrid functional^34,35, partially incorporating the Hartree-Fock exchange interaction, is used to overcome the bandgap problem with local exchange-correlation functionals. The pseudopotential is given by the projector-augmented wave method^80,81. The energy cutoff for the plane-wave basis set is 250 eV for monolayer TMDs (350 eV for diamond and monolayer hBN). We prepared a supercell of 6 × 6 × 1 primitive cells for pristine monolayer TMDs and hBN (3 × 3 × 3 cubic unit cells for diamond), including a 15-Å-thick vacuum region. The single Γ-centered k-point is adopted for the Brillouin zone sampling. A pristine cell geometry is optimized until the maximum atomic force is smaller than 0.02 eV/Å; then, a defective cell geometry is relaxed within a fixed cell shape and volume based on the optimized pristine cell. The SOC is not considered unless otherwise stated. We utilized subroutines implemented in VASP to compute the magnetic properties—the ZFS tensors and the hyperfine tensors. We used the Corrections For Formation Energy and Eigenvalues (CoFFEE) code⁴³ to calculate defect formation energies with the FNV charge correction scheme⁴².

Phonon wavefunction overlap and SOC strength

ISC is attributed to a combination of SOC and electron-phonon interaction. To obtain the phonon wavefunction overlap between the initial and final state, a one-dimensional harmonic oscillation approximation was used, which introduces the general configuration coordinate diagram. The potential surfaces of spin-triplet excited state ³E and spin-singlet state ¹A₁ were obtained by linearly interpolating between initial ³E and final ¹A₁ structures involved in the ISC. Energies of the interpolated structure were calculated using constrained-occupation DFT⁷³. Since Kohn-Sham DFT theory cannot describe states composed of multiple Slater determinates, approximate electron occupations—|a₁ e_x〉 for ³E and $|{e}_{x}{\bar{e}}_{y}\rangle$ for ¹A₁—were adopted, where ${\bar{e}}_{y}$ indicate the different spin channel of e_y orbital, and we made an approximation to access the energy of the ¹A₁ at the equilibrium geometry following Mackoit-Sinkeviciene et al.⁸². All constrained DFT computations were performed using VASP, facilitated by modified Nonrand⁸³ preprocessing and postprocessing for interpolated structure energy calculation. The calculated configuration coordinate diagram for ³E and ¹A₁ is shown in Supplementary Fig. 7.

SOC strength was computed with the ORCA code⁸⁴ using time-dependent density functional theory⁸⁵. Different from VASP, ORCA does not have the feature of periodic boundary conditions. We thus constructed cluster models for both N_CV_C⁻¹ and W_Se defects by cutting relaxed structures from bulk and saturating dangling bonds to reproduce the electronic structures of bulk structures. The dangling bonds in the diamond cluster are easily saturated by H, while TMD is well-known for complicated edge states and charge transfer between edges and defects for over 10 Å⁸⁶. After testing with different sizes, boundaries, and termination groups, a cluster with hybrid zigzag and arm-chair boundary and termination groups of H, OH, and NH was found using B3LYP functional to have both the same spin density as a periodic result [Supplementary Fig. 8a, b] and HOMO-LUMO gap of 1.22 eV to get reasonably excited states [Supplementary Figure 8(c)]. We obtained SOC values of 4.71 GHz for λ_∥ and 44.6 GHz for λ_⊥ for N_CV_C⁻¹ defect using PBE functionals with def2-TZVP basis, which agrees well with previously computed values and experimentally measured values^72,73,87. With the calculated λ_⊥, we obtained the ³E → ¹A₁ ISC rate for the NV center in diamond at 30.6 MHz which is in fair agreement with the literature-reported value of 60.7 MHz⁸⁸. We then computed the SOC strength for the axial λ_∥ and non-axial λ_⊥ components of the W_Se defect in MoSe₂ using B3LYP functionals to be 69 and 109 GHz, respectively.