Prediction of protected band edge states and dielectric tunable quasiparticle and excitonic properties of monolayer MoSi$_2$N$_4$

The electronic structure of two-dimensional (2D) materials are inherently prone to environmental perturbations, which may pose significant challenges to their applications in electronic or optoelectronic devices. A 2D material couples with its environment through two mechanisms: local chemical coupling and nonlocal dielectric screening effects. The local chemical coupling is often difficult to predict or control experimentally. Nonlocal dielectric screening, on the other hand, can be tuned by choosing the substrates or layer thickness in a controllable manner. Therefore, a compelling 2D electronic material should offer band edge states that are robust against local chemical coupling effects. Here it is demonstrated that the recently synthesized MoSi$_2$N$_4$ is an ideal 2D semiconductor with robust band edge states protected from capricious environmental chemical coupling effects. Detailed many-body perturbation theory calculations are carried out to illustrate how the band edge states of MoSi$_2$N$_4$ are shielded from the direct chemical coupling effects, but its quasiparticle and excitonic properties can be modulated through the nonlocal dielectric screening effects. This unique property, together with the moderate band gap and the thermodynamic and mechanical stability of this material, paves the way for a range of applications of MoSi$_2$N$_4$ in areas including energy, 2D electronics, and optoelectronics.


I. Introduction
It is difficult to overstate the research interest in two-dimensional (2D) materials due to their rich physics and potential applications in next-generation electronic devices [1,2]. While significant progress has been made in our understanding of the fundamental physics and properties of 2D materials, their practical applications have certainly lagged behind.
Much recent effort has been put into exploiting the interlayer coupling effects of 2D materials to tune their electronic, optical, or magnetic properties [3][4][5][6], or even to create exotic states such as those observed in twisted graphene [7][8][9].
However, the fact that the low energy electronic structure, especially the band edge states, of 2D materials are prone to environmental perturbations may become a major issue facing their practical applications in electronic or optoelectronic devices. There are two fundamental mechanisms through which a 2D material may couple with its environment: local chemical interaction and nonlocal dielectric screening. Nonlocal screening effects can be engineered with the choice of the substrate and/or the layer thickness, which can actually serve as an advantageous mechanism to tailor the quasiparticle and optical properties of 2D semiconductors in a controllable manner [10][11][12]. Local chemical coupling, which may arise from unintentional surface adsorption and/or interfacial chemical coupling, on the other hand, is particularly difficult to predict and control in experiment but can significantly affect the band edge states of 2D materials. Note that even very weak chemical coupling at typical van der Waals (vdW) distances, i.e., without the formation of conventional chemical bonds, can significantly affect the band edge states. These interactions may strongly modify the electronic structure of 2D semiconductors, thus their device performance. Therefore, a compelling 2D electronic material should offer band edge states that are robust against weak surface or interfacial chemical interactions.
Here we demonstrate that the recently synthesized MoSi2N4 [13] may provide such an ideal 2D semiconductor with band edge states being protected from capricious chemical couplings. MoSi2N4 belongs to a family of emerging 2D materials with the chemical formula MSi2N4, where M is a transition metal and Si and N may be replaced by other group IV and V elements, respectively. This material has demonstrated excellent ambient stability [13]. A number of interesting properties such as piezoelectricity and high thermal conductivity [14], and spin-valley coupling [15,16] have also been predicted. In this work, we have carried out thorough density functional theory (DFT) and many-body perturbation theory calculations to illustrate how the band edge states of MoSi2N4 are essentially shielded from direct chemical coupling to its environment, and at the same time, the quasiparticle and optical properties of MoSi2N4 can be modulated through the nonlocal dielectric screening effects which renormalize the electron-electron (e-e) and electron-hole (e-h) interactions.

A. Electronic and structural properties of MoSi2N4: from monolayer to bulk
The crystal structure of the newly synthesized monolayer MoSi2N4 is shown in Fig. 1, which can be considered as an MoN2 layer sandwiched between two SiN layers. The optimized in-plane lattice constant for monolayer MoSi2N4 is 2.896 Å, which is consistent with an earlier report [13]. The distance between two Si layers is 5.99 Å, which also agrees with experiment [13]. The strong polar covalent bonding and the relatively thick layer structure provide an energetically and mechanically stable structure for potential applications.
Such a multi-atomic layer structure also offers the possibility of having protected band edge states from weak interfacial chemical coupling if these states are primarily derived from orbitals of the interior atoms. To investigate such possibility, we carry out DFT calculations of the layer-dependent electronic structure of MoSi2N4. We first construct six  Compared with other layered materials, MoSi2N4 has slightly smaller interlayer distances d, and the interlayer binding energy also falls within that of a typical layered material.  other theoretical results [13,15]. Surprisingly, the calculated PBE band gap of bilayer MoSi2N4 is 1.77 eV, and that of the bulk phase is 1.70 eV, which are only 20 meV and 90 meV smaller than that of the monolayer structure, respectively.
The dispersion of the low energy electronic structures also shows negligible changes going from monolayer to bulk. The fact that interlayer coupling has negligible effects on the calculated band gap and the low energy band dispersion (at the DFT level) of MoSi2N4 is in stark contrast with other layered materials, as shown in Table II, considering that these materials have comparable interlayer separations and binding energies (Table I) Before we proceed, we briefly discuss the spin-orbit coupling (SOC) effect in this system.
The PBE band structure of monolayer MoSi2N4 calculated with SOC included is shown with red dotted lines in Fig. 2     The presence of small admixtures of atomic orbitals derived from Si and surface N atoms in the band edge states explains the slight layer-dependent DFT band structure and band gap as shown in Fig. 2 and Table II. Interlayer chemical coupling (hybridization) will occur for states which have significant contributions from outer atomic orbitals. Interestingly, these states are either significantly above or below the band gap. Therefore, MoSi2N4 offers an exciting material system with band edge states protected from interfacial chemical coupling. It should be pointed out that although we only discuss interlayer chemical coupling effects here, we expect that the band edge states of MoSi2N4 are similarly protected from other surface or interfacial perturbations such as physical adsorptions, or weak chemical coupling with substrates (i.e., without the formation of strong chemical bonds). As we have mentioned earlier, these perturbations are difficult to control or predict but may significantly affect the performance of 2D electronic devices.

C. Layer-dependent quasiparticle properties: Tailoring the electronic structure with nonlocal dielectric screening effects
Our results have unequivocally demonstrated that MoSi2N4 is a unique layered material with band edge states protected from surface or interfacial perturbations. This does not mean, however, that the quasiparticle or optical properties of MoSi2N4 are not affected by nonlocal interlayer or substrate coupling effects. In fact, substrate or layer-dependent dielectric screening effects can strongly renormalize the e-e or e-h interactions, thus the quasiparticle and excitonic properties of 2D materials. Fortunately, unlike local chemical coupling, which is challenging to control in experiment, substrate or interlayer screening effects can serve as an effective means to modulate the quasiparticle and optical properties of 2D materials [10][11][12]. To this end, we first carry out fully converged GW calculations for monolayer, bilayer, and bulk MoSi2N4, aiming at illustrating the dielectric screening effects on the quasiparticle properties, in particular, the quasiparticle band gap, of MoSi2N4.   [5]. For instance, the GW band gap correction for monolayer MoS2 is 0.96 eV [23], and that for bulk is 0.41 eV [24]. Those for black phosphorus are 1.17 and 0.26 eV [5]. These results again suggest that the band edge states of MoSi2N4 are not as sensitive to environmental perturbation as other 2D materials. On the other hand, controlled fine tuning of the quasiparticle band gap (from 2.82 to 2.41 eV) is still possible with the choice of substrate and/or layer thickness. Similarly, optical properties can be modulated through engineering the dielectric screening, as we will discuss in the next section.
It should be pointed out that fully converged GW calculations for 2D materials remain a challenge due to the large cell size and the unusual analytical behaviors of 2D dielectric functions and electron self-energies, which make conventional GW calculations using the band-by-band summation and the uniform Brillouin zone (BZ) integration approach rather inefficient. Our work takes advantage of the recent code developments [25,26] that can drastically speed up GW calculations for 2D materials. Using these new developments, we can effectively include all conduction bands in the GW calculations, thus eliminating the need for tedious band convergence tests [25]. In addition, the combined subsampling and analytical BZ integration approach [26] greatly improves the efficiency of the BZ integration of the quasiparticle self-energy. In this work, the self-energy integration in the BZ is carried out using a 6×6×1 k-grid with four subsampling points; these parameters are sufficient to converge the GW band gap to within 0.03 eV as discussed in our previous publications [3,22,26,27]. We have also tested the convergence of the calculated results with respect to the thickness of vacuum layer. More detail of the convergence test of our GW and BSE calculations can be found in Supplementary Material.

D. Layer-dependent electron-hole excitations and optical properties
Now we investigate the layer-dependent e-h excitations and optical properties of MoSi2N4.
The quasiparticle properties are calculated within the GW approximation as described in the previous section, and the e-h excitation spectra are obtained by solving the Bethe-Salpeter equation (BSE) [28], which is transformed into a simplified eigenvalue problem after decoupling the excitations and de-excitations (also known as the Tamm-Dancoff approximation) [28]: In our calculations, four valence and four conduction states are included in the expansion of the excitonic wave functions, which should be adequate to cover excitations of up to at least 6 eV. We have also carried calculations using more valence and conduction bands; the results for low energy excitations are practically unchanged as shown in Supplemental Material (Table S3). The contribution to an excitonic state > S | from e-h pairs in the BZ can be conveniently visualized using the k-dependent excitonic amplitude function . The imaginary part of the frequency-dependent macroscopic dielectric function including the excitonic effects is then given by (using atomic units) where λ  MoSi2N4. We include results calculated with (blue) and without (black) e-h interaction.
Including the e-h interaction not only significantly shifts the absorption edge but also produces prominent excitonic absorption peaks. Note that it is well-known the dielectric function is volumetric, therefore it scales with the thickness of the vacuum layer for 2D systems [29][30][31]. However, the excitonic features (i.e., the position and relative strength of the excitonic peaks) can be compared directly with experiment as long as the results are adequately converged with respect to the interlayer separation. The first absorption peak (i.e., peak A) at 2.18 eV for the monolayer MoSi2N4 arises from two degenerate excitonic states. This excitation energy agrees well with the measured first optical absorption peak at 2.21 eV [13]. It is interesting that two excitonic states can give rise to such a strong absorption. The k-resolved e-h pair amplitudes 2 | | S k A  for these two excitonic states are shown in Fig. 5 (d). It is clear that the dominant contribution to the first two excitonic states comes from the e-h pairs near the minimum direct gap at the K and K' points [see Fig. 4 (a) for the quasiparticle band structure of monolayer MoSi2N4]. This is very similar to that in monolayer MoS2 [32]. Therefore, for the A excitons, we can define a non-interacting electron-hole pair excitation gap   Finally, we summarize in Fig. 6 the evolution of indirect minimum band gap, the direct band gap at the K point, and the optical gap as a function of number of layers. As we have discussed earlier, the band gap of MoSi2N4 calculated at the PBE level is surprisingly stable thanks to the protected band edge states which are not susceptible to interlayer chemical coupling. Including the nonlocal self-energy correction within the GW approximation, however, leads to a significant layer-dependent quasiparticle band gap. Therefore, MoSi2N4 provides an interesting material system with robust band edge states that are spared from undesirable chemical coupling, but at the same time, tunable quasiparticle and optical properties through dielectric (or Coulomb) engineering.
Before we conclude, we would like to mention that calculations of e-h excitations at the GW-BSE level are extremely difficult to converge [23,32,34], especially with respect to the k-point sampling density. To alleviate the computational burden of BSE calculations, one often carries out the calculations on a relative coarse k-grid, the e-h kernel matrices are then interpolated to finer k-grids to obtain the final results. For the monolayer and bilayer structures, the e-h interaction kernel matrices are first calculated using a coarse 12×12×1 k-grid, the results are then interpolated onto a fine 72×72×1 k-grid from which the e-h excitations and optical absorption are obtained. For the bulk system, we use a 9×9×2 coarse k-grid and a 36×36×8 fine k-grid in our calculations. We have carefully tested the convergence of our results with respect to the k-point sampling density (Table S2 and (Table S3 in Supplementary Material).

III. Conclusion
We have performed detailed first-principles calculations for the newly discovered layered material MoSi2N4. We find that the interlayer binding energy between MoSi2N4 layers is comparable to those in MoS2 and other well-studied layered materials, but the calculated DFT-PBE band gap of MoSi2N4 is barely affected by the interlayer interaction, in stark contrast to other layered materials. Detailed analyses reveal that the wave functions of the

IV. Computational Methods
Structural optimizations for the monolayer, bilayer, and bulk phase MoSi2N4 with different layer stacking patterns are carried using the Vienna ab initio simulation package (VASP) [35,36]. The Perdew-Burke-Ernzerhof (PBE) [37,38] functional is used for the basic electronic band structure calculations, and the DFT-D3 correction of Grimme et al. [17] is employed in structural optimizations to account for the van der Waals (vdW) interactions.
To minimize the fictitious interaction between periodic image layers, a vacuum layer of over 20 Å is included in the unit cells of monolayer and bilayer systems.  [25,26] which greatly improve the efficiency of fully converged GW calculations for 2D systems. These two methods address two fundamental bottlenecks of GW calculations of 2D materials, leading to an overall speed-up factor of over three orders of magnitude compared with the conventional approaches as we have discussed in earlier publications [25,26]. Other computational details have been discussed in Results and Discussion, and also in Supplemental Material.
[39] Deslippe J, Samsonidze G, Strubbe D A, Jain M, Cohen M L and Louie S G  Figure S1 shows the six different stacking patterns of bilayer and bulk MoSi2N4 studied in this work. Table S1 shows the corresponding structural parameters. Configuration D is the most stable stacking pattern for both bilayer and bulk phases. The DFT-PBE band structures of the five metastable structures for the bilayer and bulk phases are shown in Fig. S2.    Figure S3 shows the convergence behavior of the calculated GW band gap of monolayer MoSi2N4 with respect to the kinetic energy cutoff for the dielectric matrices (the upper panel) and the number of empty states (the bottom panel) included in the GW calculations. The lower horizontal axis of the bottom panel shows the number of bands effectively included in our calculation, whereas the upper horizontal axis shows the number of integration points for evaluating the band summation using our new energy integration method [1], suggesting a speed-up factor of more than 20. The final results we reported in the manuscript are calculated using a 50 Ry cutoff for the dielectric matrix and effectively including all (more than 30,000) conduction bands, which should give well-converged results with respect to these two parameters.

Convergence tests of GW-BSE results
Fully converged GW calculations for 2D materials are more challenging than those for their 3D counterparts. GW (and BSE) results should also be checked against other convergence parameters such as k-point density for the Brillouin zone integration and the thickness of the vacuum layer included in the calculation (to reduce the spurious interlayer coupling). Both issues have been well-recognized and discussed in literature. To reduce the dependence of the calculated results on the thickness of the vacuum layer, we use the Coulomb truncation scheme proposed by Ismail-Beigi [2]. However, even with this truncation scheme, one still needs to include a sizeable vacuum layer in the calculation. In our original manuscript, the c lattice constant (perpendicular to the layer direction) for the GW/BSE calculations of monolayer MoSi2N4 is set at 30 Å, which, from our experience, is large enough to converge the results to within 0.1 eV. We then use our newly developed mini-BZ subsampling and analytical integration method [3] to evaluate the BZ integration of the self-energy.
The upper panel of Fig. S4 shows the calculated direct band gap of monolayer MoSi2N4 as a function of lattice c using different approaches. The blue curve shows the results calculated with a truncated Coulomb potential using a 6×6×1 uniform k-point sampling, whereas the red curve shows the result calculated using the 1/r long-range Coulomb potential with 6×6×1 uniform k-grid. Interestingly, using a truncated Coulomb potential would seem to give worse results if the BZ integration is not converged with respect to the k-point sampling density (i.e., using a uniform 6×6×1 grid). In other words, with the use of a truncated Coulomb potential, the GW results converge slower with respect to the BZ kpoint sampling density than using a long-range Coulomb potential. Obviously, neither results are converged. This subtle interplay among the convergence of the BZ integration, the use (or not) of truncated Coulomb potential, and the thickness of the vacuum layer has been discussed in literature, and our results are consistent with the results of Rasmussen et al. [4]. The black curve in the upper panel shows GW results as a function of the c lattice constant using our new BZ integration technique [3]. We can confidently say that the results converge to about 0.1 eV with respect to the c-lattice. Finally, the bottom panel shows the calculated GW band gap as a function of BZ k-point sampling density, with and without the use of the mini-BZ integration technique we developed [3], using a lattice constant c of 30 Å. Using our newly developed method, the calculated band gap converges even with a very coarse k-grid 3×3×1. In comparison, conventional BZ integration using a uniform k-grid would likely require an 18×18×1 (or denser) grid to achieve the same level of convergence. The results reported in our original manuscript were calculated using a 6×6×1 k-grid (with our new mini-BZ integration) and a c lattice constant of 30 Å. Table S2 shows the positions and exciton binding energies of the first excitonic peaks of monolayer, bilayer, and bulk MoSi2N4 with respect to the k-point density. The final results that we present in the main text should converge to within 50 meV with respect to the kgrid density for the monolayer and bilayer systems. For the bulk phase, the results should convergence to within 10 meV. Figure S5 shows the real and imaginary parts of the dielectric function for monolayer MoSi2N4 based on GW-BSE calculations using different fine k-grids. Finally, Table S3 shows the variation of the calculated positions of the first three excitonic peaks of monolayer MoSi2N4 with different number of valence and conduction bands included in the BSE calculations.