Electronic and Optical Properties of Vacancy Defects in Transition Metal Dichalcogenides

TL;DR

This study uses first-principle calculations to analyze how vacancy defects in single-layer transition metal dichalcogenides affect their electronic, optical, and magnetic properties, revealing defect-induced localized states and their optical signatures.

Contribution

It introduces a comprehensive analysis of vacancy defect effects in TMDCs, including localized defect states, optical responses, magnetic signatures, and a simplified tight-binding model capturing these phenomena.

Findings

Vacancy defects create localized defect states in TMDCs.

Spin-orbit coupling influences defect state splitting based on metal atomic number.

Magnetic moments arise in MoSe₂ with Mo vacancies, breaking time-reversal symmetry.

Abstract

A detailed first-principle study has been performed to evaluate the electronic and optical properties of single-layer (SL) transition metal dichalcogenides (TMDCs) (MX${}_2$; M= transition metal such as Mo, W and X= S, Se, Te), in the presence of vacancy defects (VDs). Defects usually play an important role in tailoring electronic, optical, and magnetic properties of semiconductors. We consider three types of VDs in SL TMDCs i) $X$-vacancy, $X_{2}$-vacancy, and iii) $M$-vacancy. We show that VDs lead to localized defect states (LDS) in the band structure, which in turn give rise to sharp transitions in in-plane and out-of-plane optical susceptibilities, $\chi_{\parallel}$ and $\chi_{\perp}$. The effects of spin orbit coupling (SOC) are also considered. We find that SOC splitting in LDS is directly related to the atomic number of the transition metal atoms. Apart from electronic and optical properties we also find magnetic signatures (local magnetic moment of $\sim\mu_{B}$) in MoSe$_{2}$ in the presence of Mo vacancy, which breaks the time reversal symmetry and therefore lifts the Kramers degeneracy. We show that a simple qualitative tight binding model (TBM), involving only the hopping between atoms surrounding the vacancy with an on-site SOC term, is sufficient to capture the essential features of LDS. In addition, the existence of the LDS can be understood from the solution of the 2D Dirac Hamiltonian by employing infinite mass boundary conditions. In order to provide a clear description of the optical absorption spectra, we use group theory to derive the optical selection rules between LDS for both $\chi_{\parallel}$ and $\chi_{\perp}$.

Tables

Table 1. Table 1: Calculated in-plane and out-of-plane band gaps E g ∥ subscript subscript 𝐸 𝑔 parallel-to {E_{g}}_{\parallel} and E g ⟂ subscript subscript 𝐸 𝑔 perpendicular-to {E_{g}}_{\perp} and spin-orbit splitting Δ S ​ O subscript Δ 𝑆 𝑂 \Delta_{SO} of the highest occupied valence band at K point.

System	MoS2	WS2	MoSe2	WSe2
$E_g^{}{}_{\parallel }{}^{}$[eV]	1.716	1.684	1.438	1.37
$E_g^{}{}_{⟂}{}^{}$[eV]	3.109	3.263	2.516	2.66
${\Delta }_{SO}$[meV]	150	438	195	482

Table 2. Table 2: SOC splitting Δ X , X 2 , M subscript Δ 𝑋 subscript 𝑋 2 𝑀 \Delta_{X,X_{2},M} in LDS appearing above the Fermi level for X, X 2 and M vacancies respectively, in different TMDCs

System	MoS2	WS2	MoSe2	WSe2
${\Delta }_X$[meV]	51	192	34	173
${\Delta }_{X_2}$[meV]	78	278	60	251
${\Delta }_M$[meV]	32.5	130	54	100

Table 3. Table 3: Character table of the group C 3 ​ v subscript 𝐶 3 𝑣 C_{3v} . E 𝐸 E , C 3 subscript 𝐶 3 C_{3} , σ v subscript 𝜎 𝑣 \sigma_{v} are the single group IRs and E 1 / 2 subscript 𝐸 1 2 E_{1/2} , E 3 / 2 subscript 𝐸 3 2 E_{3/2} are the corrresponding double group IRs.

$C_{3v}$	$E$	$C_3$	$3{\sigma }_v$
$A_1$	1	1	1
$A_2$	1	1	$-1$
$E$	2	$-1$	0
$E_{1/2}$	$2$ $-2$	$1$ $-1$	0 0
$E_{3/2}$	$2$ $-2$	$-2$ $2$	$0$ $0$

Table 4. Table 4: Character table of the group D 3 ​ h subscript 𝐷 3 ℎ D_{3h} . E 𝐸 E , σ h subscript 𝜎 ℎ \sigma_{h} , 2 ​ C 3 2 subscript 𝐶 3 2C_{3} , 2 ​ S 3 2 subscript 𝑆 3 2S_{3} , 3 ​ C 2 3 subscript 𝐶 2 3C_{2} , and σ v subscript 𝜎 𝑣 \sigma_{v} are the single group IRs and D 1 / 2 subscript 𝐷 1 2 D_{1/2} , 2 ​ S 1 2 subscript 𝑆 1 2S_{1} , 2 ​ S 2 2 subscript 𝑆 2 2S_{2} are the corrresponding double group IRs.

$D_{3h}$	$E$	${\sigma }_2$	$2C_3$	$2S_3$	$3C_2$	$3{\sigma }_v$
$A_1^{\prime }$	1	1	1	1	1	1
$A_2^{\prime }$	1	1	1	1	-1	-1
$A_1^{\prime \prime }$	1	-1	1	-1	1	-1
$A_2^{\prime \prime }$	1	-1	1	-1	-1	1
$E^{\prime }$	2	2	-1	-1	0	0
$E^{\prime \prime }$	2	-2	-1	1	0	0
$D_{1/2}$	2 -2	0	1 -1	$\sqrt{3}$ $-\sqrt{3}$	0	0
$2S_1$	2 -2	0	-2 2	0 0	0	0
$2S_2$	2 -2	0	1 $-1$	$-\sqrt{3}$ $\sqrt{3}$	0	0

Table 5. Table 5: Double group representations obtained from single group representaion for C 3 ​ v subscript 𝐶 3 𝑣 C_{3v} and D 3 ​ h subscript 𝐷 3 ℎ D_{3h} .

${\Gamma }_i\left(C_{3v}\right)$	$A_1$	$A_2$	$E$
${\Gamma }_i\times E_{1/2}$	$E_{1/2}$	$E_{1/2}$	$E_{3/2}+E_{1/2}$

Table 6. Table 6: Electric Dipole selection rules in C 3 ​ v subscript 𝐶 3 𝑣 C_{3v} and D 3 ​ h subscript 𝐷 3 ℎ D_{3h} symmetry. σ 𝜎 \sigma represents in plane transitions while π 𝜋 \pi represents out of plane transitions.

$C_{3v}$	$E_{1/2}$	$E_{3/2}$
$E_{1/2}$	$\sigma ,\pi$	$\sigma$
$E_{3/2}$	$\sigma$	$\sigma ,\pi$

Equations

$$\hat{H}^{TBM}_{X_{2}}=\left(\begin{array}[]{cc}\hat{H}_{e}^{X_{2}}&\hat{0}_{9\times 6}\\ \hat{0}_{6\times 9}&\hat{H}_{o}^{X_{2}}\end{array}\right),$$

$$\hat{H}_{e}^{X_{2}}=\left(\begin{array}[]{ccc}\hat{H}_{e}^{1,1}&\hat{H}_{e}^{1,2}&\hat{H}_{e}^{1,3}\\ \hat{H}_{e}^{2,1}&\hat{H}_{e}^{2,2}&\hat{H}_{e}^{2,3}\\ \hat{H}_{e}^{3,1}&\hat{H}_{e}^{3,2}&\hat{H}_{e}^{3,3}\end{array}\right),\quad\hat{H}_{o}^{X_{2}}=\left(\begin{array}[]{cc}\hat{H}_{o}^{4,4}&\hat{H}_{o}^{4,5}\\ \hat{H}_{o}^{5,4}&\hat{H}_{o}^{5,5}\end{array}\right),$$

$$\hat{H}_{e(o)}^{\alpha,\beta}=\left(\begin{array}[]{ccc}\epsilon_{e(o)}^{\alpha,\beta}&t_{e(o)}^{\alpha,\beta}e^{i\theta}&t_{e(o)}^{\alpha,\beta}e^{-i\theta}\\ t_{e(o)}^{\alpha,\beta}e^{-i\theta}&\epsilon_{e(o)}^{\alpha,\beta}&t_{e(o)}^{\alpha,\beta}e^{i\theta}\\ t_{e(o)}^{\alpha,\beta}e^{i\theta}&t_{e(o)}^{\alpha,\beta}e^{-i\theta}&\epsilon_{e(o)}^{\alpha,\beta}\end{array}\right).$$

$$\hat{H}^{SOC}_{X_{2}}=\frac{\Delta}{2}\textbf{L}\cdot\textbf{S}=\frac{\Delta}{2}\left(\begin{array}[]{cc}\hat{L}_{z}&\hat{L}_{-}\\ \hat{L}_{+}&-\hat{L}_{z}\end{array}\right),$$

$$\hat{L}_{z}=\left(\begin{array}[]{ccccc}\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}\\ \hat{0}_{3}&\hat{0}_{3}&2i\times\hat{I}_{3}&\hat{0}_{3}&\hat{0}_{3}\\ \hat{0}_{3}&-2i\times\hat{I}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}\\ \hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&-i\times\hat{I}_{3}\\ \hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&i\times\hat{I}_{3}&\hat{0}_{3}\\ \end{array}\right),$$

$$\hat{L}_{+}=\left(\begin{array}[]{ccccc}\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&\sqrt{3}\times\hat{I}_{3}&i\sqrt{3}\hat{I}_{3}\\ \hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&-i\times\hat{I}_{3}&-1\times\hat{I}_{3}\\ \hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&-1\times\hat{I}_{3}&i\times\hat{I}_{3}\\ -\sqrt{3}\times\hat{I}_{3}&i\times\hat{I}_{3}&\hat{I}_{3}&\hat{0}_{3}&\hat{0}_{3}\\ -i\sqrt{3}\times\hat{I}_{3}&\hat{I}_{3}&-i\times\hat{I}_{3}&\hat{0}_{3}&\hat{0}_{3}\\ \end{array}\right)$$

$$|\Psi ⟩=\alpha |\zeta ⟩|\uparrow ⟩+\beta |\xi ⟩|\downarrow ⟩.$$

$$\begin{split}\hat{H}_{X_{2}}&=I_{2}\otimes\hat{H}_{X_{2}}^{TBM}+\hat{H}_{X_{2}}^{SOC}\\ &=\left(\begin{array}[]{cc}\hat{H}_{X_{2}}^{TBM}+\frac{\Delta}{2}\hat{L}_{z}&\hat{0}_{9\times 9}\\ \hat{0}_{9\times 9}&\hat{H}_{X_{2}}^{TBM}-\frac{\Delta}{2}\hat{L}_{z}\end{array}\right).\end{split}$$

$${\hat{H}}_M^{TBM}$$

$${\hat{H}}_e^M$$

$${\hat{H}}_o^M$$

$$\hat{L}_{z}=\left(\begin{array}[]{cccccc}\hat{0}_{3}&-i\times\hat{I}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}\\ i\times\hat{I}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}\\ \hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}\\ \hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&-i\times\hat{I}_{3}&\hat{0}_{3}\\ \hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&i\times\hat{I}_{3}&\hat{0}_{3}&\hat{0}_{3}\\ \hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}\\ \end{array}\right),$$

$$\hat{L}_{+}=\left(\begin{array}[]{cccccc}\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&-1\times\hat{I}_{3}\\ \hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&i\times\hat{I}_{3}\\ \hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{I}_{3}&-i\times\hat{I}_{3}&\hat{0}_{3}\\ \hat{0}_{3}&\hat{0}_{3}&-1\times\hat{I}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}\\ \hat{0}_{3}&\hat{0}_{3}&i\times\hat{I}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}\\ \hat{I}_{3}&-i\times\hat{I}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}&\hat{0}_{3}\\ \end{array}\right),$$

$${\mathcal{H}}_{\tau }=\tau {\sigma }_z\Delta +v\mathbf{\sigma }\cdot \mathbf{p}.$$

$$\Phi ={\Phi }_{+}\oplus {\Phi }_{-}$$

$$\mathcal{H}={\tau }_z\otimes {\sigma }_z\Delta +v{\tau }_0\otimes \mathbf{\sigma }\cdot \mathbf{p},$$

$$\sigma_{r}=\left(\begin{array}[]{cc}0&e^{-i\phi}\\ e^{i\phi}&0\end{array}\right),\quad\sigma_{\phi}=\left(\begin{array}[]{cc}0&-ie^{-i\phi}\\ ie^{i\phi}&0\end{array}\right).$$

$${\tilde{\mathcal{H}}}_{\tau }{\tilde{\Phi }}_{\tau }=ϵ{\tilde{\Phi }}_{\tau },$$

$$\tilde{\mathcal{H}}=\tau {\sigma }_z\Delta -iv\left[{\sigma }_x\left(\frac{\partial }{\partial r}+\frac{1}{2r}\right)+{\sigma }_y\frac{1}{r}\frac{\partial }{\partial \varphi }\right].$$

$$\begin{split}\widetilde{\Phi}_{\tau;m}&(r)=\sqrt{Qr}\\ \times&\left(\begin{array}[]{c}\frac{1}{\sqrt{\epsilon-\tau\Delta}}\left(a_{\tau;m}h_{m-1}^{(1)}(Qr)+b_{\tau;m}h_{m-1}^{(2)}(Qr)\right)\\ \frac{i}{\sqrt{\epsilon+\tau\Delta}}\left(a_{\tau;m}h_{m}^{(1)}(Qr)+b_{\tau;m}h_{m}^{(2)}(Qr)\right)\end{array}\right),\end{split}$$

$$\widetilde{\Phi}_{+,m}=N_{+,m}\left(\begin{array}[]{c}\frac{1}{\sqrt{\Delta-\epsilon}}g_{m-1}(\kappa r)\\ \frac{i}{\sqrt{\Delta+\epsilon}}g_{m}(\kappa r)\end{array}\right),$$

$$g_m(z)=(-z{)}^m{\left(\frac{d}{zdz}\right)}^m\frac{e^{-z}}{z}.$$

$$g_{-m}(z)=g_{m-1}(z),$$

$${\tilde{\Phi }}_{-,-m}=c{\sigma }_x{\tilde{\Phi }}_{+,m},$$

$$\widetilde{\Phi}_{+,m}=\widetilde{N}_{+,m}\left(\begin{array}[]{c}\cos(\chi_{m}/2)\\ i\sin(\chi_{m}/2)\end{array}\right),$$

$$F_m=\frac{g_m(\kappa r)\sqrt{\Delta -ϵ}}{g_{m-1}(\kappa r)\sqrt{\Delta +ϵ}},$$

$$F_m(ϵ)=1/F_{-m}(-ϵ).$$

$$F_m(ϵ,r_0)=1.$$

$${ϵ}_0=0.$$

$$r_0>R_m=\frac{v}{\Delta }\left(|m|-\frac{1}{2}\right).$$

$$\tan \left(\frac{\beta }{2}\right)=\tau \frac{\sqrt{\Delta +ϵ}-\sqrt{\Delta -ϵ}}{\sqrt{\Delta +ϵ}+\sqrt{\Delta -ϵ}}.$$

Full text

I Introduction

Single-layer (SL) transition metal dichalcogenides (TMDCs) have attracted a lot of attention due to their intriguing electronic and optical properties, with a wide range of promising applications.Wang et al. (2012); Mak and AU Shan (2016) SL TMDCs are direct band gap semiconductorsMak et al. (2010); Splendiani et al. (2010), which can be used to produce smaller and more energy efficient devices, such as transistors and integrated circuits. Moreover, the band gap lies in the visible region which makes them highly responsive when exposed to visible light, a property with potential applications in optical detection. In contrast to graphene, SL TMDCs exhibit large intrinsic spin-orbit coupling (SOC), originating from the d orbitals of transition metal atoms. The presence of considerably high SOC (up to few hundred meV)Xiao et al. (2012); Liu et al. (2013); Song and Dery (2013) makes them a candidate material for exploring spin physics and spintronics applications.

Wafer-scale production of SL TMDCs is required to fully appreciate their technological potential. The most common experimental techniques used to produce large chunks of SL MoS2 are i) mechanical exfoliation, ii) chemical vapor deposition, and iii) physical vapor deposition. It has been observed that samples produced by all of these techniques have considerably lower carrier mobility than the theoretically predicted values.Radisavljevic et al. (2011); Fuhrer and Hone (2013) It has recently been suggested that this discrepancy between the predicted and observed values of carrier mobility is due to the presence of impurities created during the growth process.Yu et al. (2014); Komsa et al. (2012a) The most common and energetically favorable types of impurities are vacancy defects (VDs).Hong et al. (2015)

Defects usually play an important role in tailoring various electronic and optical properties of two dimensional materials and have been the subject of intense research over the last few decades. VDs in semiconductors act as trapping centers for charge carriers and their interaction with charge carriers becomes stronger at reduced dimensionalities. Point defects in SL TMDCs have been explored both theoretically and experimentally.Kunstmann et al. (2016); Pulkin and Yazyev (2016); Li et al. (2016); Srivastava et al. (2015); He et al. (2015); Koperski et al. (2015) Recent photoluminescence (PL) experiments Srivastava et al. (2015); He et al. (2015); Koperski et al. (2015) reveal that localized excitonic states related to VDs can serve as single-photon emitters in WSe2. Magnetism in low dimensional systems is another area of interest.Yazyev (2010); Yazyev and Helm (2007) It has been shown thatMa et al. (2011) certain LDS related to VDs can induce ferromagnetism in SL TMDCs, suggesting that they could be good candidates for spin channels in spintronic devices. In addition, LDS can be used to open and tune a band gap in graphenePedersen et al. (2008) and SL MoS2.Huang et al. (2013); Shao et al. (2014) Various atomic defects can be realized artificially by using different experimental techniques. It has been shown that hexagonal pits (3$\times$MoS2) can be removed through etching of MoS2 crystals by using XeF2 as a gaseous reactant.Huang et al. (2013) Point defects can be induced by irradiating the SL TMDCs with $\alpha$-particles or by thermal annealing.Tongay et al. (2013) Several experimental studies have been reported regarding the effects of point defects or of grain boundaries on SL TMDCs.Nan et al. (2014); Zhou et al. (2013); Komsa et al. (2012b) Strong PL enhancement has been observed as a result of oxygen adsorption at sulfur vacancy sites.Zhou et al. (2013) Also, sulfur vacancies are observed in MoS2 through transmission electron microscopy experiments.Komsa et al. (2012b)

Pristine TMDCs are invariant with respect to the reflection $\sigma_{h}$ about the Mo or W plane of atoms ($z=0$ plane). Therefore, electron states can be classified into two catagories: even and odd or symmetric and antisymmetric with respect to the $z=0$ plane. We found that the even and odd bands in TMDCs have two different band gaps $E_{g\parallel}$ and $E_{g\perp}$, respectively.Erementchouk et al. (2015) $E_{g\perp}$ has been usually neglected for pristine TMDCs because of its substantially larger value (Table 1) and weak optical response (Fig: 2 b) as compared with $E_{g\parallel}$. Earlier studiesErementchouk et al. (2015); Shao et al. (2014) show that the presence of VDs gives rise to LDS in addition to the normal extended states present in conduction or valence bands in SL MoS2. These LDS appear within the band gap region and they can also be present deep inside the valence band depending on the type of VD. Optical transitions between LDS across Fermi level appear as resonance peaks, both in $\chi_{\parallel}$ and $\chi_{\perp}$, which shows that odd states are necessary for understanding the properties of VDs in SL MoS2.Erementchouk et al. (2015)

In this paper, our aim is fourfold. First, we provide a comprehensive study of VDs in 4 types of SL TMDC materials, MoS2, MoSe2, WS2, and WSe2. Second, we provide detailed analytical models about the description of LDS within the Dirac equation formulation and by using the tight binding model. Third, we include the effects of SOC on VDs, which has not been considered so far. As mentioned earlier, SOC in these materials is large and therefore needs to be taken into account in order to obtain a better understanding of the electronic and optical properties of TMDCs. Fourth, we briefly discuss defect induced magnetism in some cases. Throughout this work, we consider 3 types of VDs: i) Single $X$-vacancy ii) $X_{2}$-vacancy, and iii) $M$-vacancy.

This paper is organized as follows. Section II describes the numerical results obtained for band structures. Sections III and IV describes qualitative models for the existence of defected states. Section V deals with the optical response of defected SL TMDCs.

II Band Structure

The model system consists of a periodic 2D superlattice of TMDCs (Fig: 1 a, b, c). All numerical calculations are carried out using density functional theory (DFT). The local density approximation (LDA) is used with the Perdew-Zunger (PZ) parametrizationPerdew and Zunger (1981) of the correlation energy of a homogeneous electron gas calculated by Ceperley-Alder Ceperley and Alder (1980). The calculations are implemented within Atomistix Toolkit 2015.1QW_ in order to be able to perform DFT calculations on large supercells in a reasonable amount of time. The periodic structure of the superlattice allows one to characterize the electron states by the bandstructure $\epsilon_{n}(\mathbf{k})$, where $\mathbf{k}$ is the vector in the first Brillouin zone of the superlattice and $n$ enumerates different bands. We consider a $7\times 7\times 1$ (Fig: 1) supercell having $147$ number of atoms with an edge length of $21.354$ Å. The Brillouin zone of the supercell is sampled by a $3\times 3\times 1$ $k$-mesh. All the structures are geometrically optimized with a force tolerance of $0.05$ eV/Å. SOC is taken into account via the norm conserving pseudo potentials. Hamann et al. (1979); Bachelet et al. (1982) Band structures are calculated along the $\Gamma$ $-$ $M$ $-$ $K$ $-$ $\Gamma$ path. Band structures of SL TMDCs for the pristine cases are plotted in Fig. 2 and calculated values are given in Table 1. The results are in good agreement with previously reported values both for band gap and SOC energy. Xiao et al. (2012); Liu et al. (2013); Zhu et al. (2011); Zibouche et al. (2014) We consider LDA because it is computationally less expensive and therefore allows us to perform DFT calculations on large supercells. A drawback of the generalized gradient approximation (GGA) is that the Atomistix Toolkit 2015.1 gives rise to an indirect band gap for SL TMDCs, which is in contradiction to the already established results for TMDCs. Nonetheless, we obtain approximately the same values for both band gap and SOC using either LDA or GGA. Figures 4 and 5 show the band structure of various SL TMDCs in the presence of vacancies. Black lines denote regular electronic states within the valence or conduction bands while colored lines denote the LDS. Vertical arrows show some of the allowed optical transitions observed in the optical spectra (see Fig. 10).

III Tight Binding Model (TBM) and Symmetries

III.1 General considerations

The simplest qualitative model that can explain the existence of LDS in the band structure due to VDs is the TBM. Within the TBM aproximation the electron wavefunction can be presented as $\ket{\psi}=\sum_{j,\mu\in O_{j}}\psi_{\mu}^{(j)}\varphi_{\mu}^{(j)}(\mathbf{r}-\mathbf{R}^{(j)})$, where $j$ enumerates atomic positions surrounding the vacancy and $\mu$ runs over the atomic orbitals $O_{j}$. In our tight binding analysis only the atoms surrounding the VD are considered in order to make the calculations simple enough for capturing the essential physical properties of the problem. The three VDs can be classified into two groups on the basis of symmetries. The $X$-vacancy lacks spatial inversion symmetry with respect to the $M$-plane of atoms, i.e. the $\sigma_{h}$ symmetry is broken, and is therefore described by the group $C_{3v}$. In contrast, the $X_{2}$ and $M$-vacancy preserve the $\sigma_{h}$ symmetry of the crystal and thus can be described by the group $D_{3h}$.Cheiwchanchamnangij et al. (2013); Song and Dery (2013) For the latter the electronic states break down into even and odd parity with respect to the $\sigma_{h}:z\mapsto-z$ symmetry. $d$-orbitals of the transition metal and $p^{(t,b)}$- orbitals ($t$ and $b$ denoting the top and bottom layers) of the chalcogen atoms give the largest contribution to the conduction and valence band structure of TMDCs. Zhu et al. (2011); Cappelluti et al. (2013) Based on the $\sigma_{h}$ symmetry, the even and odd atomic orbitals are spanned by the bases $\{d_{x^{2}-y^{2}},d_{xy},d_{z^{2}},{~{}}p_{x,y}^{e}=(p_{x,y}^{(t)}+p_{x,y}^{(b)})/\sqrt{2},{~{}}p_{z}^{e}=(p_{z}^{(t)}-p_{z}^{(b)})/\sqrt{2}\}$ and $\{d_{xz},d_{yz},{~{}}p_{x,y}^{o}=(p_{x,y}^{(t)}-p_{x,y}^{(b)})/\sqrt{2},{~{}}p_{z}^{o}=(p_{z}^{(t)}+p_{z}^{(b)})/\sqrt{2}\}$, respectively. We also include the effects of intrinsic SOC of the form $\sim\textbf{L}\cdot\textbf{S}$. The resulting spin-orbit states transform according to irreducible representations (IRs) of the double groups $C^{D}_{3v}$ and $D^{D}_{3h}$. Group representation theory is an efficient tool for determining the allowed optical transitions across the Fermi level in solids. This will be discussed in detail in the last section. The aim of this section is to present a qualitative description of LDS appearing in the band structure (Fig: 2 b, c, d). Here following Refs. Liu et al., 2013; Cappelluti et al., 2013; Erementchouk et al., 2015 we first develop the TBM Hamiltonian by allowing the hopping between atomic orbitals of the atoms surrounding the VD only. Also we consider a large supercell in order to suppress the intervacancy couplings. Consequently, the effects of SOC are considered as VD onsite couplings.

III.2 X- and X2-vacancy

Both X- and X2-vacancies are surrounded by 3 M atoms, as shown in Fig. 3. As mentioned earlier for M atoms, $d$-orbitals provide the main contribution. Considering 3 atomic sites $A,B,C$ with 5 $d$-orbitals on each site, we have 15 species of $d$-electrons. We will suppress the spin indices and denote electron operators collectively as a vector by $\psi=(\psi_{1},\psi_{2},\psi_{3},\psi_{4},\psi_{5})$, with $\psi_{\tau}=(d_{\tau}^{A},d_{\tau}^{B},d_{\tau}^{C})$, where $d_{\tau}^{P}$ denotes the annihilation operator of electrons for orbital $\tau$ at site $P$ with $\tau=1,2,3,4,5$ standing for $d_{z^{2}},d_{xy},d_{x^{2}-y^{2}},d_{xz},d_{yz}$, respectively. The spinless representation of the Hamiltonian can be expressed in block form as

[TABLE]

where $\hat{H}_{e}^{X_{2}}$ and $\hat{H}_{o}^{X_{2}}$ are $9\times 9$ and $6\times 6$ blocks with even ($e$) and odd ($o$) parity, respectively, with respect to $\sigma_{h}$, and $\hat{0}_{m\times n}$ denotes a zero matrix of dimension $m\times n$. Here we take advantage of the inversion symmetry $\sigma_{h}$ by decoupling the orbitals with opposite parities.

Also orthogonality between different orbitals on the same atomic site is enforced. The submatrices in Eq. (1) are given by

[TABLE]

where each of $\hat{H}_{e(o)}^{i,j}$ is a $3\times 3$ matrix $i,j={1,2,3,4,5}$. The blocks in the above Hamiltonians read

[TABLE]

The diagonal elements ${\epsilon_{e(o)}^{\alpha,\beta}}$ (0 for $\alpha\neq\beta$) and the off-diagonal elements $t_{e(o)}^{\alpha,\beta}$ are phenomenological parameters describing the onsite energy and hopping between $d$-orbitals at different atomic sites, respectively. It can be easily shown that Hamiltonian (3) is invariant under $\sigma_{h}$ and C3 symmetry operations, for $\theta=0,\pm 2\pi/3$. But, in addition to $\sigma_{h}$ and C3 symmetry operations, D3h group also contains $\sigma_{v}$ symmetry operations, i.e. reflection by the planes perpendicular to the xy-plane and passing through the lines OA, OB and OC (Fig. 3). $\sigma_{v}$ demands all the complex factors appearing in Equation (3) to be 1 or equivalently $\theta=0$. Equations (2, 3) provides an initial insight into the nature of LDS. One can easily show that e.g. $\hat{H}_{o}^{X_{2}}$ has a pair of 3 eigenvalues, i.e. $\bar{\epsilon}+t-\sqrt{(\delta\epsilon+h)^{2}+4t_{h}^{2}}$, $\bar{\epsilon}-t/2-\sqrt{(\delta\epsilon-h/2)^{2}+t_{h}^{2}}$, $\bar{\epsilon}-t/2-\sqrt{(\delta\epsilon-h/2)^{2}+t_{h}^{2}}$ and $\bar{\epsilon}+t+\sqrt{(\delta\epsilon+h)^{2}+4t_{h}^{2}}$, $\bar{\epsilon}-t/2+\sqrt{(\delta\epsilon-h/2)^{2}+t_{h}^{2}}$, $\bar{\epsilon}+t/2-\sqrt{(\delta\epsilon-h/2)^{2}+t_{h}^{2}}$, where $\bar{\epsilon}$, $\delta\epsilon$ are related to addition and subtraction of onsite energies for orbitals dxz, dyz; $t$, $h$ are related to addition and subtraction of hopping parameters of the same orbitals at different sites, and $t_{h}$ is the hopping parameter of different orbitals at different atomic sites. Each pair contains a two-fold doublet, which explains the existence of triplets within the band structure.Erementchouk et al. (2015) However, the apparent two-fold degeneracy, which arises from the overlap of neighboring atomic orbitals, is lifted in the presence of SOC. Here we emphasize that each $d$-orbital appears in the form of triplets in the band structure. Thus, there is a total of 15 LDS (in the absence of SOC) for the case of the X2-vacancy. It may appear that the simplest TBM may contradict the numerical results in Fig. 4, where calculations show a lower number of LDS. A closer inspection of the numerical results, however, resolves this contradiction in favor of the TBM. In fact, in addition to the LDS appearing within the band gap region, there are also LDS deep inside the valence bands, with possibility to mix with the extended states in the bulk.

SOC in the Hamiltonian is included by a pure atomic termLiu et al. (2013) and for simplicity we consider only the onsite contribution arising from the M atoms surrounding the vacancy. Using the basis $\ket{d_{z^{2}},\uparrow},\ket{d_{xy},\uparrow},\ket{d_{x^{2}-y^{2}},\uparrow},\ket{d_{xz},\uparrow},\ket{d_{yz},\uparrow}$ and $\ket{d_{z^{2}},\downarrow},\ket{d_{xy},\downarrow},\ket{d_{x^{2}-y^{2}},\downarrow},\ket{d_{xz},\downarrow},\ket{d_{yz},\downarrow}$, we can write the SOC Hamiltonian as

[TABLE]

where

[TABLE]

[TABLE]

and $\hat{L}_{-}=\hat{L}_{+}^{\dagger}$. The off-diagonal elements $\hat{L}_{\pm}$ in Eq. (4) couple the even to the odd blocks of the Hamiltoian matrix shown in Eq. (1) and are related to the spin flip processes due to the SOC, which give rise to virtual transitions.Roldán et al. (2014) Because of the large spatial anisotropy of an atomically thin layer of TMDC, for the pristine case these off-diagonal terms can be neglected, which is substantiated by our DFT calculations (see below). A generalized SOC states has the form

[TABLE]

Here, $\ket{\zeta}$ and $\ket{\xi}$ are orbital states of the spin-up and spin-down states $\ket{\uparrow}$ and $\ket{\downarrow}$, respectively, and $\alpha$, $\beta$ are probablility amplitudes for the up and down spinors. DFT calculations reveal that for SOC Bloch states corresponding to LDS, either $\alpha\ll\beta$ or $\alpha\gg\beta$ in the majority of cases (Fig: 6), corresponding to strong polarizations of the LDS in z direction (b), which is due to the large spatial anisotropy. It can be calculated that LDS are spin polarized for X2-vacancy or the Bloch states for X2-vacancy preserve the $\sigma_{h}$ symmetry. Therefore we anticipate that the effects of $\hat{L}_{\pm}$ can safely be neglected for the X2-vacancy. The full tight binding Hamiltonian can be written as

[TABLE]

The Hamiltonian appears to be block diagonal, which indicates that spin states in $z$ direction are not mixed by spin flip processes and therefore the spin in $z$ direction is still a good quantum number due to $\sigma_{h}$ symmetry. As mentioned above, the abscence of spin flip processes can be attributed to the 2D character of TMDCs or due to the large anisotropy between $xy$-plane and $z$-axis. In case of the $X$-vacancy, due to lack of $\sigma_{h}$ symmetry, defect states appear with no definite parity (Fig: 6) (a). Therefore, here we argue that for the $X$-vacancy the off-diagonal terms $\hat{L}_{\pm}$ in Eq. (4) need to be taken into account.

In the absence of SOC each energy band is doubly degenerate (spin-up and spin-down states at each k point). SOC lowers the symmetry and can break the spin degeneracy at k points away from high symmetry points. However, time reversal symmetry leads to the condition that $\varepsilon(\textbf{k},\uparrow)=\varepsilon(-\textbf{k},\downarrow)$, commonly known as Kramers degeneracy. This degeneracy is reflected in the band structure (Fig. 4), where each energy level is doubly degenerate for both types of vacancies. In solids or 2D surfaces spin splitting depends both on the size of atomic SOC and of the gradient of electric potential.Petersen and Hedegård (2000) This difference in the gradient of electric potential leads to the different spin splittings for same types of defects in different TMDCs as shown in Fig. 4 and in Table 2.

III.3 M-vacancy

There are 6 chalcogen X atoms in the top and bottom layers, surrounding the tranisition metal M-vacancy. Thus, there are 18 species of electrons corresponding to six possible combinations of $p$-orbitals (3 of them even and 3 odd with respect to $\sigma_{h}$) at 3 in-plane atomic positions. Proceeding as before the TBM for the M-vacancy can be written as

[TABLE]

where each of $\hat{H}_{e(o)}^{j}$ is a $3\times 3$ matrix corresponding to even (odd) combinations of $p$ orbitals, with $j={6,7,8,9,10,11}$ being indices reserved for the ${p_{x}^{e},p_{y}^{e},p_{z}^{e},p_{x}^{o},p_{y}^{o},p_{z}^{o}}$ orbitals, respectively. Each $\hat{H}_{e(o)}^{\beta}$ in Eq. (15) has the same form as in Eq. (3).

SOC is included by considering the Hamiltonian described in Eq. (4). The M-vacancy also preserves the $\sigma_{h}$ symmetry. Again, the $\hat{L}_{z}S_{z}$ term gives the largest contribution due to the large anisotropy. The Bloch states are shown in Fig. 6, from which it can be concluded that LDS are spin polarized also in the case of M-vacancy. The matrices for operators $\hat{L}_{z}$ and $\hat{L}_{\pm}$ in the case of the M-vacancy can be written as

[TABLE]

,

[TABLE]

and $\hat{L}_{-}=\hat{L}_{+}^{\dagger}$. The DFT calculations show that the Kramers degeneracy is preserved for MoS2 and WS2 while it is broken for MoSe2 and WSe2. In Ref. Ma et al., 2011, it has been shown that presence of Mo vacancies in MoSe2 can induce spin polarization and results in long range antiferromagnetic coupling between local magnetic moments, even at a distance above 13 Å due to the large spatial extensions of spin density. The local magnetic moment on each M-vacancy breaks the time reversal symmetry and therefore lifts Kramers degeneracy in MoSe2 in the presence of Mo vacancies. In Fig. 5 (b) we show plots for the density of states (DOS) obtained by the local density spin approximation (LSDA) method for the M-vacancy in different TMDCs, in order to confirm that indeed the Mo-vacancy in MoSe2 exhibits magnetic signatures. Our spin-polarized DFT calculations show that the exchange correlation potential leads to a spin splitting only for the MoSe2 system. In Fig. 7(a) the isosurface plot for the spin density is shown along with the magnetic moment $\mu$ calculated by means of the Mulliken population analysis at all the nearest-neighboring atomic sites (Se atoms) and the next-nearest neighboring sites (Mo atoms) surrounding the vacancy. Our results show that the main contribution to the magnetism is due to the $p$-orbitals localized at the Se atoms and the $d$-orbitals localized at the next-nearest Mo atoms surrounding the Mo vacancy. The electronic band structure calculated by using LSDA is shown in Fig. 7(b). The spin states of the LDS are split whereas the bulk states do not show any magnetic moment. The magnetic moment of the LDS is governed by the unpaired electron spins according to Hund’s rules, which is 1 $\mu_{B}$ (inset of Fig. 7(b)). The calculations yield a slightly smaller magnetic moment of 0.6$\mu_{B}$, which is acceptable within standard DFT limits. Finally, we plot the magnetic moment $\mu$ vs defect density $\rho$ in Fig. 7(c), where the localization of the magnetic moment is demonstrated for densities $\rho<25\times 10^{12}$ cm*-2*. We also notice splittings in the case of WSe2, but no magnetic signatures are found. We attribute this splitting as a result of interaction between adjacent vacancies due to large spatial extensions of W and Se orbitals.

IV Defect states within the Dirac equation formalism

The main analytical tool for describing properties of electron states in transition metals monolayers is the Dirac equation, which emerges within the $k\cdot p$-approximation as the two-band model.Xiao et al. (2012); Kormányos et al. (2015, 2014) Such description is valid when the main role is played by low lying excitations near the band edges. This assumption, however, is ill-justified for the case of vacancies and, indeed, as will be demonstrated below the Dirac equation fails to reproduce many important features of the defect states. At the same time, the Dirac equation allows one to establish several important features and, first of all, the sole existence of the states bound to vacancies. For example, within the framework provided by the Schrodinger equation the vacancy is naturally represented by a repulsive potential, which cannot support bound states. The Dirac equation, in turn, provides more options for describing defects and, as will be shown below, when the special boundary conditions are enforced at the boundary of the defect, the states localized near the boundary appear.

The formalism of Dirac equation can be introduced as follows. Within the $k\cdot p$-approximation the electron states are described by $\Psi(\mathbf{r})e^{i\mathbf{K}\cdot\mathbf{r}}$, where $\Psi(\mathbf{r})$ is a smooth function of coordinates. Adopting the two-band approximation, $\Psi(\mathbf{r})$ is presented as a two-component spinor, which satisfies a 2D Dirac-like equation. For example, for MoS2 the two-band approximation is often implemented retaining only the dominating contribution of Mo’s $d$ orbitals,Xiao et al. (2012); Liu et al. (2013) so that near the inequivalent $K$-points of the Brillouin zone, $\mathbf{K}_{\tau}=\tau\mathbf{K}$ with $\tau=\pm 1$, spinors $\Psi_{\tau}(\mathbf{r})$ are spanned by $\ket{d_{x^{2}-y^{2}}-i\tau d_{xy}}$ and $\ket{d_{z^{2}}}$, the states representing the top of the valence band and the bottom of the conduction band, respectively. Spinors $\Psi_{\tau}$, thereby, satisfy $\left[\sigma_{z}\Delta+v\left(\tau\sigma_{x}p_{x}+\sigma_{y}p_{y}\right)\right]\Psi_{\tau}(\mathbf{r})=\epsilon\Psi(\mathbf{r})$, where $2\Delta$ is the width of the gap and $\epsilon$ is the energy counted from the center of the gap $\epsilon_{c}$.

In order to eliminate the valley dependence of the Hamiltonian governing the spatial distribution of $\Psi_{\tau}$, it is convenient to introduce $\Phi_{+}=\Psi_{+}$ and $\Phi_{-}=\sigma_{y}\Psi_{-}$, which satisfy $\mathcal{H}_{\tau}\Phi_{\tau}=\epsilon\Phi_{\tau}$, where

[TABLE]

Thus solutions for electrons in different valleys are related by simple reverting the sign of $\Delta$. Combining $\Phi_{\pm}$ into a single 4-spinor

[TABLE]

the equations of motion for different valleys can be presented in a unified form $\mathcal{H}\Phi=\epsilon\Phi$ with

[TABLE]

where $\tau_{i}$ with $i=x,y,z$ and $\tau_{0}$ are the Pauli matrices and the identity matrix, respectively, acting on the valley space.

Hamiltonian (24) possesses the cylindrical symmetry, which can be employed by presenting $\boldsymbol{\sigma}\cdot\mathbf{p}=-i\sigma_{r}\partial/\partial r-ir^{-1}\sigma_{\phi}\partial/\partial\phi$, where

[TABLE]

The explicit angular dependence is eliminated by introducing $\widetilde{\Phi}_{\tau}=e^{i\sigma_{z}\phi/2}\Phi_{\tau}$, which accounts for rotation of the spinors $\Phi_{\tau}$, while encircling the origin. It should be noted that due to the relation $\Phi_{-}=\sigma_{y}\Psi_{-}$ the rotation directions of $\Psi_{+}$ and $\Psi_{-}$ are different: $\widetilde{\Psi}_{\tau}=e^{i\tau\sigma_{z}\phi/2}\Psi_{\tau}$. Thus the winding numbers of spinors corresponding to electrons belonging to different valleys have opposite signs.

Spinors $\widetilde{\Phi}_{\tau}$ satisfy

[TABLE]

where $\widetilde{\mathcal{H}}_{\tau}=e^{i\sigma_{z}\phi/2}\mathcal{H}_{\tau}e^{-i\sigma_{z}\phi/2}$ has the form

[TABLE]

Equation (26) is solved by separating variables $\widetilde{\Phi}_{\tau}(r,\phi)=\sum_{m=-\infty}^{\infty}\widetilde{\Phi}_{\tau;m}(r)e^{im\phi}$. For amplitudes $\widetilde{\Phi}_{\tau;m}(r)$, we find the general solution

[TABLE]

where $h_{m}^{(1,2)}(Qr)$ are the spherical Hankel functions, $Q=v^{-1}\sqrt{\Delta^{2}-\epsilon^{2}}$, $a_{\tau;m}$ and $b_{\tau;m}$ are arbitrary constants.

We are interested in bound states and, therefore, in solutions of Eq. (26) corresponding to energies inside the gap. For such energies, we have $Q=i\kappa$ with non-negative $\kappa=\sqrt{\Delta^{2}-\epsilon^{2}}$. From the regularity condition at infinity, it follows that $b_{\tau}^{(m)}=0$, while $a_{\tau}^{(m)}$ are determined from the normalization condition. The solution can be written as

[TABLE]

where we have denoted the normalization constant by $N_{+,m}$. The functions $g_{m}(z)$ are related to the modified spherical Hankel functions $g_{m}(z)=2k_{m}(z)/\pi$ and for $m>0$ can be presented as

[TABLE]

Taking into account the relation

[TABLE]

we can use Eq. (30) for finding $g_{m}(z)$ with $m<0$ as well. With the help of this relation, one can show, starting from Eq. (28), that

[TABLE]

where $c$ is a phase factor, $|c|=1$. Such connection between solutions corresponding to electrons from different valleys allows us to limit our attention to $\tau=+1$.

The functions $g_{m}(z)$ can be shown to be non-negative. Thus, we can rewrite

[TABLE]

with $0\leq\chi_{m}\leq\pi$. This representation shows that at any chosen distance from the center of the vacancy, the defect states have the form of a spin coherent state Aravind (1999); Combescure and Robert (2012) lying in the plane perpendicular to $\mathbf{n}_{B}$, the vector normal to the boundary of the anti-dot and directed outward. The angle $\chi_{m}=2\arctan(F_{m})$, where

[TABLE]

has the meaning of the polar angle of the vector characterizing the direction of the spin coherent state. Its dependence on $m$ is illustrated by Fig. 8, which shows that $\chi_{m}$ monotonously increases from $\chi_{-\infty}=0$ to $\chi_{\infty}=\pi$. It is also a monotonous function of $r$ (increasing for $m<0$ and decreasing for $m>0$) and monotonously increasing function of energy. Taking into account Eq. (31) one can see the important symmetry

[TABLE]

Thus, the states with $\tau=1$ lie in the half plane corresponding to the positive projection on the vector tangent to the boundary $\mathbf{e}_{\phi}=\mathbf{e}_{z}\times\mathbf{n}_{B}$, while the states with $\tau=-1$ have the opposite orientation, due to $\innerproduct{\widetilde{\Phi}_{-,m}}{\widetilde{\Phi}_{+,m}}=0$, which can be easily checked.

The energies of the defects states (and their sole existence) are determined by the boundary condition on the boundary of the anti-dot. The general form of the condition is found requiring that the radial component of the probability current must vanish at the boundary,Berry and Mondragon (1987); Tworzydło et al. (2006) $\matrixelement{\Phi(r_{0})}{\mathbf{n}_{B}\cdot\mathbf{J}}{\Phi(r_{0})}=0$, where $\mathbf{J}=v\tau_{0}\otimes\boldsymbol{\sigma}$. This condition is equivalent to $M\Phi=\Phi$, where $\Phi$ is the 4-spinor defined by Eq. (23) and the Hermitian matrix $M$ has the eigenvalues $\pm 1$ and anticommutes with the radial component of the current operator $\{\mathbf{J}\cdot\mathbf{n}_{B},M\}=0$. Within the infinite mass model,Berry and Mondragon (1987); Tworzydło et al. (2006) the anti-dot is represented as a region with renormalized width of the gap $\Delta\to\Delta(1+d(r))$ with $d(r)=0$ for $r>r_{0}$ and $d(r)\to\infty$ when $r<r_{0}$, so that in this case $M=(\boldsymbol{\tau}\cdot\mathbf{e}_{z})\otimes(\boldsymbol{\sigma}\cdot\mathbf{e}_{\phi})$. In other words, within this model in order to have decaying electron distribution inside the anti-dot $\widetilde{\Phi}_{\tau}(r)$ must be proportional to $\ket{\tau\mathbf{e}_{\phi}}$ as $r$ approaches $r_{0}$.

The condition $\widetilde{\Phi}_{+,m}(r_{0})\propto\ket{\mathbf{e}_{\phi}}$, or $\chi_{m}=\pi/2$, constitutes the condition imposed on the energy of the bound state

[TABLE]

In virtue of Eq. (35), if for some $m$ there exists a bound solution with the energy $\epsilon$, then there is the solution corresponding to $m^{\prime}=-m$ with the energy $-\epsilon$. Thus within the infinite mass model the spectrum of the defect states is symmetric with respect to the center of the gap.

For $m=0$, Eq. (36) has the simplest form $F_{0}=\sqrt{\Delta-\epsilon}/\sqrt{\Delta+\epsilon}=1$ with the solution

[TABLE]

Thus the anti-dot independently of its size supports a bound state with the energy at the center of the gap.

States with $|m|>0$, in turn, appear only when the defect is sufficiently large. In order to find the condition of supporting the state with some $m$ we notice that $F_{m}(\epsilon,r_{0})$ is monotonously decreasing function of energy while $\epsilon$ changes from $-\Delta$ to $\Delta$. Since $g_{m}(z\to 0)\sim(2m-1)!!/z^{m+1}$, we find that the energy of the $m$-th state is inside the gap, if

[TABLE]

Conversely, for the given radius $r_{0}$ the number of bound defect states is given by $N=4+8\lceil r_{0}\Delta/v+1/2\rceil$, where $4$ accounts for states from different valleys and with different spins at $m=0$ and the second term accounts for states with $m>1$, here $\lceil\ldots\rceil$ denotes taking the integer part and $8$ in addition to the spin and valley degeneracies accounts for the symmetry $m\to-m$. The dependence of energies of the defect states on the radius of the anti-dot is shown in Fig. 9.

The spatial electron distribution corresponding to the defect state is conveniently characterized by the probability density $\rho_{\tau,m}(r)=\innerproduct{{\Phi}_{\tau,m}}{{\Phi}_{\tau,m}}$ and the vector of orientation of the (pseudo)spin coherent state $\mathbf{S}_{\tau,m}=\matrixelement{{\Phi}_{\tau,m}}{\boldsymbol{\sigma}}{{\Phi}_{\tau,m}}/\rho_{\tau,m}$. As follows from Eq. (33), the pseudospin state is transversal, $S_{r}=0$, with the spatial variation of the projection of $\mathbf{S}_{\tau,m}$ onto the $(\mathbf{e}_{\phi},\mathbf{e}_{z})$-plane depending on $m$.

In the simplest case $m=0$ the pseudospin remains in the plane of the layer, $S_{\tau,z}=0$. States with nonzero $m$ are characterized by out of the plane distribution of the pseudospin (for $r>r_{0}$). The angle of maximum deviation from the plane is

[TABLE]

Thus for $\tau=1$ the pseudospin “sticks out” of the plane for $\epsilon>0$ (that is for $m>0$) and has the negative projection on the $z$-axis for $\epsilon<0$. For $\tau=-1$ the direction of the pseudospin is reversed.

V Optical Response

The presence of LDS in the band structure gives rise to sharp peaks in the optical spectrum. In Ref. Li et al., 2014 the relative dielectric functions $\epsilon_{r}$ of various TMDCs have been measured experimentally, which in turn are related to the electric susceptibilities by the standard formula $\epsilon_{r}=1+\chi$. In Fig. 10 results for both in plane $\chi_{\parallel}$ and out of plane $\chi_{\perp}$ components of the electric susceptibility are presented for different VDs in TMDCs. The electric susceptibility provides valuable insight into the optical selection rules for transitions between states across the Fermi level. We are interested in transitions involving states with energy near the gap edges or inside the gap. The electric susceptibility tensor is evaluated using the Kubo-Greenwood formula

[TABLE]

where $\pi_{pq}^{j}=\matrixelement{{\psi}_{p}}{x_{j}}{\psi_{q}}$ is the $j$th component of dipole matrix element between states $p$ and $q$, $V$ the volume, $f$ the Fermi function and $\Gamma$ is the broadening, which is set to be 0.01 eV. Appearance of states inside the band gap $E_{g}$ or close to the band edges leads to the resonances at single energy $E_{pq}=|E_{p}-E_{q}|$. The dipole matrix element $\pi_{pq}^{j}$ determines the strength of an optical transition and whether it is allowed or prohibited by symmetries.

When considering defects in a crystal, the LDS transform according to the IRs of the symmetry group of the crystal site in which the defect resides. While the translational symmetry of the crystal is broken, point group symmetries are partially or completely preserved. M and X2 vacancies keep the $D_{3h}$ symmetry whereas the X-vacancy exhibits the lower $C_{3v}$ symmetry. The character tables for $C_{3v}$ and $D_{3h}$ with single and double group IRs are shown in Table 3 and Table 4, respectively. Table 5 contains the decomposition of the direct products of the single group representations with the representation according to which the spin matrices transform, i.e. $E_{1/2}$ and $D_{1/2}$ respectively. In Fig: 11 band structures for WS2 with S and W vacancies are shown for with and without SOC. Following Ref. Dresselhaus, 1957; Birman, 1959, in Fig: 11 we show how the single group IRs in the absence of SOC can be mapped to the corresponding double group IRs in the presence of SOC. Note that due to the nature of the DFT calculation, the superlattice defined by the supercell introduces an artificial translational symmetry, which in some cases leads to artificial splittings. These can be typically recognized by systematically changing the size of the supercell.

The appearance of LDS inside the band gap leads to sharp resonances in $\chi_{\parallel}$ and $\chi_{\perp}$ at frequencies corresponding to the energy differences between LDS. However, not all transitions are allowed. Instead, several transitions are prohibited due to symmetry, i.e. when $\pi_{pq}^{j}$ does not transform according to the symmetric representation of the symmetry group of the superlattice. In the matrix element $\chi_{pq}^{j}$, the initial state $\psi_{p}$, the final state $\psi_{q}$, and the position operator $x_{j}$ transform according to the IRs $\Gamma(\psi_{p})$, $\Gamma(\psi_{q})$ and $\Gamma(x_{j})$, respectively. An electric dipole transition between two states is allowed if the direct product $\Gamma(\psi_{p})\otimes\Gamma(x_{j})\otimes\Gamma(\psi_{q})$ includes $\Gamma(I)$ in its decomposition in terms of a direct sum. $\Gamma(I)$ denotes the IR for the identity i.e., $A_{1}$ and $A_{1}^{\prime}$ for $C_{3v}$ and $D_{3h}$, respectively. This is strictly related to the polarization of the radiation. One needs to consider separately the in plane and out of plane components of $\pi_{pq}^{j}$ because they transform according to different IRs of $C_{3v}$ and $D_{3h}$. The selection rules for electric dipole transitions for the double group IRs are summarized in Table 6.

The presence of SOC couples the spin and orbital angular momenta, thereby requiring the consideration of the double group IRs. In our case, we need to consider $D_{1/2}\otimes C_{3v}$ and $D_{1/2}\otimes D_{3h}$, where $D_{1/2}$ is the 2-dimensional spin representation. The electromagnetic field couples to the orbital part of the state, i.e. either to $\left|\zeta\right>$ or to $\left|\xi\right>$ [see Eq. (7)], polarizations induced by electromagnetic waves will not be changed in the presence of SOC. The role of the SOC is to lift some degeneracies, which gives rise to extra absorption peaks compared with the case without SOC.Erementchouk et al. (2015) In the susceptibility, these extra peaks lie close to the energies predicted by the susceptibility without SOC. It is important to take care when dealing with selection rules described by double groups because double groups may allow some transitions that are prohibited by the single groups. Then such selection rules must be discarded. One such example is the $\pi$-transition for $X$-vacancy. In the abscence of SOC the $\pi$-transition is allowed only between states with the IRs $E$, i.e. $E\otimes A_{1}\otimes E=A_{1}\oplus A_{2}\oplus E$; this transition is not allowed since orbitals of the IR $E$ exist above the Fermi level only. In the presence of SOC the $\pi$-transition is allowed by the double group, i.e. $E_{1/2}\otimes A_{1}\otimes E_{3/2}=A_{1}\oplus A_{2}\oplus E$, but is not seen in the susceptibility including SOC. This can be understood as the artefact of double groups since $\pi$-transitions are prohibited in the abscence of SOC.

As a final note, as mentioned earlier, all samples are geometrically optimized before performing electronic calculations. Geometrical optimization may break certain symmetries and can affect certain selection rules or can result in concealing of some of the resonances.

VI Conclusion

In this paper we have provided numerical and analytical descriptions of electronic and optical properties of SL TMDCs in the presence of VDs. We have shown that the presence of LDS gives rise to sharp transitions both in $\chi{\parallel}$ and $\chi_{\perp}$. In order to understand these transitions, odd states need to be considered in addition to even states. A central result of our paper is that group theory can be used to derive strict selection rules for the optical transitions, which are in excellent agreement with the susceptibility calculated using the Kubo-Greenwood formula using the DFT orbitals. SOC induced splitting is observed in LDS and is seen to be larger for VDs in WX2 than in MoX2. Interestingly, our findings suggest magnetic properties of MoSe2 in the presence of Mo vacancy, which may be enhanced by increasing the density of defects. In order to provide a qualitative explanation of the existence of LDS, we performed analytical calculations based on the TBM and 2D Dirac formulation. All these results considerably improve the understanding of VDs in SL TMDCs and should benefit their potential applications in optoelectronic and nanoelectronic devices.

VII Acknowledgments

We acknowledge support provided by NSF CCF-1514089, the Airforce Summer Faculty Fellowship 2016, and AFOSR LRIR 15RY159COR. We thank Kavir Dass, Shin Mou, and Xiaodong Xu for useful discussions.