Efficient and accurate defect level modelling in monolayer MoS$_2$ via GW+DFT with open boundary conditions

Within the framework of many-body perturbation theory (MBPT) integrated with density functional theory (DFT), a novel defect-subspace projection GW method, the so-called p-GW, is proposed. By avoiding the periodic defect interference through open boundary self-energies, we show that the p-GW can efficiently and accurately describe quasi-particle correlated defect levels in two-dimensional (2D) monolayer MoS$_2$. By comparing two different defect states originating from sulfur vacancy and adatom to existing theoretical and experimental works, we show that our GW correction to the DFT defect levels is precisely modelled. Based on these findings, we expect that our method can provide genuine trap states for various 2D transition-metal dichalcogenide (TMD) monolayers, thus enabling the study of defect-induced effects on the device characteristics of these materials via realistic simulations.


Introduction
The physical dimension of Si logic transistors is approaching the atomic limit, thus requiring novel architectures and/or high-mobility channel materials for future technology nodes. Logic switches based on two-dimensional (2D) transitionmetal dichalcogenide (TMD) monolayers have thus been proposed to continue Moore's scaling law, thanks to their remarkable electronic properties. However, several works [1,2] reported that various defects inside these monolayers may limit their performance as logic devices, mainly through charged impurity scattering and defect-induced trap levels. In particular, the "mid-gap" states introduced by those impurities are presumably at the origin of large Schottky barriers (SB) and high contact resistances. Therefore, in order to understand the physics related to defects in 2D TMD monolayers and to guide device design, ab initio simulations are required. In this work, we propose an efficient GW algorithm combined with density functional theory (DFT) to accurately describe defect levels in 2D TMD monolayers. In conventional GW calculations, environmental effects from substrates are included to obtain the realistic bandgap of 2-D monolayers, which requires huge computational resources [3]. Our method, so-called projected GW (p-GW), overcomes this issue by projections onto a defect subspace while removing spurious interactions between periodic images by means of open boundary conditions. This algorithm can correctly predict the position of defect levels in the bandgap and ensure efficiency by resorting to the DFT-level bandgap. We then apply this method to the most common defects in MoS 2 monolayers: S vacancy and adatom.

Algorithm
The p-GW algorithm is based on Green's function theory and aims at describing isolated defects. The starting point is the definition of a device region Ω D containing a defect and consisting of integer repetitions of a unit cell called "principal layer" (PL), as illustrated in Fig. 1. We want to compute a device Green's function G D which includes r e g i o n defect Figure 1: Schematic view of a device region Ω D containing a defect and consisting of integer repetitions of a PL. The coupling to the bulk states of the surrounding material is described by Σ B . Electron-electron interactions beyond the mean-field level of theory are included only for a narrow region surrounding the defect and denoted as "GW region". the correlation of the electrons localized around the defect and couples to the Bloch states of pristine MoS 2 at the boundaries. This is achieved through the following procedure: 1. A DFT calculation of Ω D builds the Hamiltonian of the defect+MoS 2 system at a mean-field level. 2. A boundary self-energy replaces the periodic boundary conditions (PBCs) of DFT. 3. Projection onto an orthogonal subspace surrounding the defect defines the GW region.
Our p-GW algorithm has been developed for the general case of a non-orthogonal basis set employed in the DFT calculation. The flowchart in Fig. 2 shows how the above three steps can be used to obtain G D . We describe each step in more details below.

Boundary self-energy
The boundary self-energy Σ B describes the coupling of Ω D with the electrons in the Bloch states of MoS 2 . Ω D must be choosen large enough to safely assume that the potential at the boundaries is converged to the one of pristine MoS 2 . Σ B can then be efficiently obtained from a k-point calculation of the PL as summarized in the following (see also Refs. [4,5] for more details): 1. An iterative method such as Sancho-Rubio [6] calculates the PL surface Green's function along one direction. The last step requires the evaluation of a disconnected Green's function given by where z is a complex energy with an infinitesimal shift along the imaginary axis. H and S are the Hamiltonian and overlap matrices of the defect-free Ω D region, which can be conveniently obtained from the corresponding PL matrices exploiting the translational symmetry of the supercell structure [5]. This efficient and precise algorithm allows us to treat the system as "open" and effectively simulates the defect as isolated. Indeed, this avoids undesired interferences or bound state patterns related to the PBCs.

Projection onto a GW region
GW corrections are computed only for a narrow region Ω C surrounding the defect, as shown in Fig. 1. We start by selecting a region Ω P ⊂ Ω D where the formation of electron-hole pairs is expected to screen the electrons in the defect level states. Because of the strong atomic orbital character of these states, this region can only comprise the defect and a few atoms nearby. Taking into account the non-orthogonality between the DFT basis set, we write the Green's function projected onto Ω P as [7] where S DP is the overlap matrix between orbitals in Ω D and in Ω P , respectively, and S P is the overlap matrix within Ω P . Starting from G P and the bare Coulomb matrix calculated for the basis functions in Ω P , we obtain the screened Coulomb interaction W P within the random phase approximation [8,9]. Finally, we multiply the part of the screened interaction W C in Ω C by the Green's function G C projected from G P in analogy with Eq. 2 to obtain the GW self-energy Σ GW .

Device Green's function
The device Green's function allows to access the electronic and transport properties of the defected monolayers. It is obtained as [10] where H and S are the device Hamiltonian and overlap matrices with removed PBCs and where V xc is the DFT exchange-correlation (XC) potential that needs to be subtracted to avoid double counting of the correlations included in Σ GW . δV H is the deviation from the DFT Hartree potential and is calculated from the change in the density matrix D C in the Ω C region. Because δV H and Σ GW depend on G D themselves, Eq. 3 is solved self-consistently until convergence of D C .

Results
We study the effect of S vacancies (S-) and adatoms (S+) in 2D MoS 2 monolayers. The device region is composed of 4 × 6 repetitions of a PL composed of 6 Mo and 12 S atoms, as shown in Fig. 1. The electronic structure calculation of the PL is over-sampled with a 11 × 6 × 1 k-mesh to obtain a Σ B that precisely describes the bulk MoS 2 states. The GW region is shown in Figs. 3a and 3b together with the wavefunction of the states created by the defect. Ω P includes up to the 2nd nearest neighbor to the defect, i.e. 12 Mo and 13 S atoms for S-and 3 Mo and 15 S atoms for S+. The defect states have a strong atomic orbital character: they are essentially a superposition of the 3d Mo orbitals closest to the vacancy for S-and an unpaired electron in the in-plane p orbital of the S adatom for S+. This allows us to define Ω C as the 3 Mo and their surrounding S atoms for S-and the single S adatom for S+. Σ C is then computed in this region only.  We calculated the correspondingdensity-of-states (DOS), the projected DOS (PDOS), and the electron transmission and report these results in Fig. 4. It is apparent from the DOS and the PDOS that the effect of the many-body correction is to shift the energy levels of the defect while preserving the DFT properties, i.e. the bandgap, as also corroborated by the conservation of the bulk-like electronic transmission. Previous k-point GW studies of full defect+MoS 2 S-structures found similar positions of the defect level with respect to the corresponding band edges [3]. This indicates that our p-GW algorithm can accurately predict trap-levels with minimal computational burden. The DFT study for S+ predicts a shallow state close to the valence band. The GW correction pulls the defect-level position down into the valence band, as suggested by experimental studies that show a strong p-type behaviour in the presence of S adatoms [11]. This indicates that such defects may act as doping centre, a behaviour that is not captured by pure DFT.

Conclusions
We proposed a novel algorithm to locally and efficiently apply many-body corrections using GW to a region surrounding a defect. Periodic self-interactions are removed by virtue of an efficient boundary self-energy calculation. The presented algorithm is then applied to S vacancy and adatom defects in a MoS 2 monolayer. Our method is a first step toward the inclusion of many-body methods beyond DFT in large scale simulations of realistic devices.    The onset of the electron transmission around the fundamental gap is also preserved by the p-GW correction.