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 devices1,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 properties3,4,5,6, or even to create exotic states such as those observed in twisted graphene7,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 manner10,11,12,13. 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 MoSi2N414 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 stability14. A number of interesting properties such as piezoelectricity and high thermal conductivity15, and spin-valley coupling16,17 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 excitonic 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.

Results

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 report14. The distance between two Si layers is 5.99 Å, which also agrees with experiment14. The strong polar covalent bonding and the relatively thick layer structure provide an energetically and mechanically stable structure for potential applications.

Fig. 1: Side view and top view of the crystal structure of monolayer MoSi2N4.
figure 1

The surface and interior N atoms are shown with red and blue colors, respectively.

Full size image

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 bilayer structures and the corresponding bulk structures with different stacking patterns. The structures are optimized using the PBE functional with the DFT-D318 correction to account for the interlayer vdW interaction. Details of the optimized structures are shown in the Supplementary Fig. 1 and Supplementary Table 1. The calculated inter-layer separation ranges from 3.270 to 2.825 Å for the bilayer structure, and 3.222 to 2.807 Å for the bulk phase. Not surprisingly, the lowest energy structures have the smallest interlayer separations. Table 1 compares the calculated interlayer separation d and binding energy of the most stable bilayer and bulk MoSi2N4 with a few selected layered materials. Experimental values for the bulk interlayer distances, if available, are also shown for comparison. 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.

Table 1 Comparison between the calculated interlayer binding energy (Ebind = (2Emonolayer − Ebilayer)/A, in meV/Å2) and interlayer separation d (in Å) of MoSi2N4 with a few selected layered materials using the PBE functional with the DFT-D3 correction.
Full size table

It is well known that interlayer chemical coupling can significantly modify the low energy electronic structure, especially the band gap, of many vdW layered materials. Such chemical coupling effects can be well described within DFT provided that the interlayer separations are properly determined. With the small interlayer separations and moderate interlayer binding energy, we would expect that MoSi2N4 experiences similar chemical coupling effects found in other layered materials.

Figure 2 compares the PBE band structure of monolayer, bilayer, and bulk MoSi2N4 calculated using the optimized structures. Monolayer MoSi2N4 has an indirect band gap of 1.79 eV calculated within PBE; the valence band maximum (VBM) locates at the Γ point and the conduction band minimum (CBM) locates at K. These results are consistent with other theoretical results14,16. 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 DFT-PBE band structures of the five metastable structures for the bilayer and bulk phases are shown in Supplementary Figure 2.

Fig. 2: DFT-PBE calculated band structures for MoSi2N4 from monolayer to bulk.
figure 2

The blue solid and red dotted lines in (a) represent the band structure without and with SOC effects for monolayer MoSi2N4, respectively.

Full size image

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 2, considering that these materials have comparable interlayer separations and binding energies (Table 1). For example, the DFT band gap of MoS2 changes from 1.72 eV (monolayer) to 0.88 eV (bulk), and C3N and C3B both show a band gap change of over 1.3 eV from monolayer to bulk.

Table 2 Interlayer chemical coupling effects on the band gap of MoSi2N4 compared with other layered materials.
Full size table
Table 3 Layer-dependent direct and indirect quasiparticle band gaps (in eV) of MoSi2N4 calculated within the GW approximation.
Full size table
Table 4 Layer-dependent GW band gap correction (ΔΣ) (to the direct gap at the K point), binding energy of the A exciton (\(E_{bind}^A\)), and optical gap (\(E_g^{opt}\)) of MoSi2N4. (in eV).
Full size table

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(a). The SOC effects result in a splitting of 130 meV of the top valence band at the K point; the SOC effects on the CBM (at the K point) and VBM (Γ) states are negligible. These results agree well with experiment14.

Robust band edge states protected from interfacial chemical coupling

The negligible layer dependence of the calculated DFT gap strongly hints at the presence of robust near-edge electronic states that are shielded from interfacial chemical coupling. This is only possible if the band edge states are derived primarily from interior (i.e., Mo and interior N atoms) atomic orbitals. In order to gain further understanding, we show in Fig. 3 (a) the decomposition of the Kohn-Sham wave functions near the band gap into contributions from atomic orbitals. It is obvious that the band edge states are largely derived from Mo d orbitals with a small admixture of other atomic orbitals. This is also clearly seen in Fig. 3 (b) and (c), in which we show the isosurface charge density plots of the band edge states at the K point: the charge density of the band edges states mainly locates inside the MoSi2N4 layer. Considering the relatively thick monolayer structure and the fact that Mo d orbitals are fairly localized, interlayer chemical coupling effects on the band edge states are greatly reduced.

Fig. 3: Atomic character of the near band edge states of MoSi2N4.
figure 3

a Decomposition of Bloch wave function of MoSi2N4 into atomic contributions, b the isosurface charge density plot \(( {\rho = 0.02e/{{{\mathrm{{\AA}}}}}^3})\) for the top of the valence band at the K point, and c the isosurface charge density plot for the bottom of the conduction band at the K point.

Full size image

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 2. 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.

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 excitonic 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 materials10,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.

Figure 4 compares the GW band structures of monolayer, bilayer, and bulk MoSi2N4. The band gaps of all three structures remain indirect (Γ to K) in nature within the GW approximation. Although the DFT band gap of MoSi2N4 shows a very small layer dependence, changing from 1.79 to 1.70 eV from the monolayer to the bulk phase as we have discussed earlier, the quasiparticle band gap shows more a substantial layer dependence as summarized in Table 3. The calculated GW band gaps (indirect) are 2.82 eV for monolayer, 2.67 eV for bilayer, and 2.41 eV for bulk MoSi2N4, representing a self-energy correction ranging from 1.03 (monolayer) to 0.71 (bulk) eV. Note that even calculated at the GW level, the layer-dependence of the quasiparticle band gap of MoSi2N4 is still significantly weaker than other well-known 2D materials such as MoS219,20 and black phosphorus5. For instance, the GW band gap correction for monolayer MoS2 is 0.96 eV19, and that for bulk is 0.41 eV20. Those for black phosphorus are 1.17 and 0.26 eV5. 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, excitonic properties can be modulated through engineering the dielectric screening, as we will discuss in the next section.

Fig. 4: GW band structures for MoSi2N4 from monolayer to bulk.
figure 4

a–c represent the GW band structure of monolayer, bilayer, and bulk MoSi2N4, respectively. Important electron and hole states involved in the optical transitions that give rise to the A (red), B (green), and C (blue) low energy absorption peaks in monolayer MoSi2N4 are also shown.

Full size image

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 developments21,22 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 tests21. In addition, the combined subsampling and analytical BZ integration approach22 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 publications3,22,23,24. We have also tested the convergence of the calculated results with respect to the thickness of the vacuum layer. More detail of the convergence test of our GW and BSE calculations can be found in Supplementary Figs. 3 and 4.

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)25, which is transformed into a simplified eigenvalue problem after decoupling the excitations and de-excitations (also known as the Tamm-Dancoff approximation)25:

$$(E_{c\vec k} - E_{v\vec k})A_{vc\mathop{k}\limits^{\rightharpoonup} }^S + \mathop {\sum}\limits_{v^\prime c^\prime \mathop{k}\limits^{\rightharpoonup}{\,}^{\!\!\prime }} { \,<\, vc\mathop{k}\limits^{\rightharpoonup} |\hat K^{eh}|{v^\prime c^\prime \mathop{k}\limits^{\rightharpoonup}{\,}^{\!\!\prime}}> A_{v^\prime c^\prime \mathop{k}\limits^{\rightharpoonup}{\,}^{\!\!\prime }}^S} = \Omega^{S}A_{vc\mathop{k}\limits^{\rightharpoonup} }^S,$$
(1)

where \(\hat K^{eh}\)is the e-h interaction kernel, \(E_{c\vec k}\) \((E_{v\vec k})\) is the quasiparticle energy of conduction (valence) states, S labels the excitonic state \(|S >\) with eigenvalue \(\Omega ^S\)and eigenstate function constructed using the electron and hole wave functions with the e-h coefficient \(A_{vc\mathop{k}\limits^{\rightharpoonup} }^S\):

$$\psi ^S(\vec r_e,\vec r_h) = \mathop {\sum}\limits_{vc\mathop{k}\limits^{\rightharpoonup} } {A_{vc\mathop{k}\limits^{\rightharpoonup} }^S\varphi _{c\vec k}(\vec r_e)\varphi _{v\vec k}(\vec r_h)}$$
(2)

In our calculations, four (eight) valence and four (eight) conduction states are included in the expansion of the excitonic wave functions for the monolayer (bulk or bilayer) MoSi2N4, which should be adequate to cover excitations of up to at least 6 eV. We have also carried out calculations using more valence and conduction bands; the results for low energy excitations are practically unchanged as shown in Supplementary Table 3. 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 \(|A_{\mathop{k}\limits^{\rightharpoonup} }^S|^2 = \mathop {\sum}\nolimits_{vc} {|A_{vc\mathop{k}\limits^{\rightharpoonup} }^S|^2}\). The imaginary part of the frequency-dependent macroscopic dielectric function including the excitonic effects is then given by (using atomic units)

$$\varepsilon _2(\omega ) = \frac{{16\pi ^2}}{{\omega ^2}}\mathop {\sum}\limits_S {|\vec \lambda \cdot \,<\, 0|\vec v|S \,>\, |^2\delta (\omega - \Omega ^S)}$$
(3)

where \(\vec \lambda\) is the unit vector of the polarization direction of light and \(\vec v\) is the velocity operator. A Gaussian smearing of width 25 meV is used in the calculation of \(\varepsilon _2(\omega )\).

Figure. 5 compares the calculated \(\varepsilon _2(\omega )\)for monolayer (a), bilayer (b), and bulk (c) 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 that the dielectric function is volumetric, therefore it scales with the thickness of the vacuum layer for 2D systems26,27,28. 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.50 eV for the monolayer MoSi2N4 arises from two degenerate excitonic states, to be compared with the reported first optical absorption peak at 2.21 eV14. In addition to the underlying accuracy of the theory, factors that may contribute to the difference between experiment and theory include temperature and SOC effects; both effects turn to shift the absorption energy downwards. It is interesting that two excitonic states can give rise to such a strong absorption. The k-resolved e-h pair amplitudes \(|A_{\mathop{k}\limits^{\rightharpoonup} }^S|^2\)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 MoS229. Therefore, for the A excitons, we can define a non-interacting electron-hole pair excitation gap \(E^{{{{\mathrm{non - int}}}}} = E_{c{{{\mathrm{K}}}}}^{QP} - E_{v{{{\mathrm{K}}}}}^{QP}\), which is the minimum quasiparticle band gap at the K point. Measuring from this non-interacting electron-hole pair excitation gap, we obtain the exciton binding energy of about 0.63 eV for the first two excitonic states. A more rigorous definition for the exciton binding energy is the difference between the mean quasiparticle excitation gap (for a given excitonic state) and the excitonic eigenvalue, i.e.,

$$E_{bind}^S = \mathop {\sum}\limits_{vc\mathop{k}\limits^{\rightharpoonup} } {|A_{vc\mathop{k}\limits^{\rightharpoonup} }^S|^2(E_{c\mathop{k}\limits^{\rightharpoonup} }^{QP} - } E_{v\mathop{k}\limits^{\rightharpoonup} }^{QP}) - \Omega ^{S}.$$
(4)

Using this definition, we obtain an exciton binding energy of 0.78 eV for the A exciton for the monolayer MoSi2N4.

Fig. 5: Layer-dependent optical properties of MoSi2N4.
figure 5

a–c show the absorption spectra of MoSi2N4 from monolayer to bulk, and (d–f) e-h pair amplitudes of the first three absorption peaks of monolayer MoSi2N4.

Full size image

Peak B (at around 2.79 eV) consists of eight energetically close lying (i.e., nearly degenerate) excitonic states; these excitonic states are derived from e-h pairs around the Γ point and near the K (K’) point with close quasiparticle band gaps. This can be clearly seen in Fig. 5 (e), which shows \(\mathop {\sum}\nolimits_{S \in \{ B\} } {|A_{\mathop{k}\limits^{\rightharpoonup} }^S|^2}\), summation of the k-resolved e-h amplitudes of the eight excitonic states that contribute to the B absorption peak. Another pronounced adsorption peak (peak C) appears at about 2.97 eV. This absorption peak comprises ten bright excitonic states. The summation of the e-h amplitudes, \(\mathop {\sum}\nolimits_{S \in \{ C\} } {|A_{\mathop{k}\limits^{\rightharpoonup} }^S|^2}\), of these states is shown in Fig. 5 (f), suggesting that this absorption peak primarily originates from e-h pairs states with wave vectors near the Γ point along the Γ-M directions, forming a ring-like structure around the Γ point with a small mixture of e-h pairs near the K (K’) point. These features are consistent with the GW band structure shown in Fig. 4. Note that SOC effect is not included in these calculations. Including SOC will result in a double-peak structure of the first excitonic peak with a splitting of about 130 meV. The SOC splitting for the B and C peaks, however, is substantially smaller (about 50 meV). This is because the e-h pairs involved in the formation of the excitons that are responsible for the B and C absorption peaks spread out near the Γ and K (K’) points [see Fig. 5(e) and (f)] where the SOC splitting is significantly reduced at the one-particle level as we have discussed earlier.

Figure 5 (b) and (c) show the imaginary part of the dielectric function of the bilayer and bulk phases. Similar to the case of the monolayer, three strong low energy absorption peaks can be clearly identified for the bilayer system. Due to the increasingly stronger dielectric screening, therefore weakened e-h interaction, the calculated exciton binding energies decrease to 0.40 eV and 0.12 eV for the bilayer and bulk phases, respectively. The position of the absorption peak A shifts slightly upwards with the increasing layer thickness: 2.50 eV for monolayer, 2.53 eV for bilayer, and 2.59 eV for bulk. These shifts, however, are likely within the numerical accuracy of the calculations. The fairly stable optical band gap (with respect to the layer thickness) can be traced to the fact that the band edge states of MoSi2N4 are protected from the interlayer coupling, as we have discussed in sections II A and B. Additional results comparing the percentage absorbance between the monolayer and bilayer MoSi2N4, and the imaginary part of the dielectric function calculated using a larger broadening parameter of 40 meV are shown in Supplementary Figs. 6 and 7.

Compared with other 2D materials, the DFT band gap of MoSi2N4 does not experience a significant reduction with the increasing number of layers, as shown in Table 3. The calculated DFT direct gap at the K point actually increases slightly from 2.03 eV (monolayer) to 2.10 eV (bulk). The GW quasiparticle correction and the excitonic binding energy partially cancel out, as shown in Table 4. This leads to a net correction to the optical excitation gap \((\Delta \Sigma - E_{bind}^A)\) (from the DFT value) that is nearly unchanged going from monolayer to bulk. As a result, the optical gap of MoSi2N4 is also rather insensitive to the number of layers or environmental perturbations, which may be advantageous for applications that require a stable optical gap. We would like to emphasize that the partial cancellation between the quasiparticle self-energy correction and the excite binding energy alone cannot fully explain a stable optical gap in this system. There are plenty of example of 2D materials with significant layer dependent optical gaps; examples include GeSe30, black phosphorus5, C3N and C3B.

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 the 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 that 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 excitonic properties through dielectric (or Coulomb) engineering.

Fig. 6: Layer-dependent band gap results of MoSi2N4.
figure 6

Evolution of the indirect minimum band gap (left panel), the direct band gap at the K point and the optical gap (right panel) as a function of the number of layers calculated at different theory levels.

Full size image

Before we conclude, we would like to mention that calculations of e-h excitations at the GW-BSE level are extremely difficult to converge19,29,31, 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 24 × 24 × 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 an 18 × 18 × 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 (Supplementary Table 2, Supplementary Figs. 4 and 5). We believe the results should converge to within about 50 meV. We have also tested the adequacy of the number of valence and conduction bands included in the expansion of the excitonic wave functions (Supplementary Table 3).

Discussion

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 band edge states of MoSi2N4 are derived primarily from interior atomic orbitals (i.e., Mo and interior N atoms) and thus are shielded from surface or interfacial chemical perturbations. The quasiparticle and excitonic properties of MoSi2N4 can still be modulated through the nonlocal dielectric screening effects, as demonstrated by our GW-BSE calculations. The quasiparticle band gap of MoSi2N4 varies from 2.82 eV for the monolayer to 2.41 eV for the bulk phase. Interestingly, the optical band gap actually increases slightly with the increasing number of layers, changing from 2.50 eV (monolayer) to 2.60 eV (bulk). Our findings suggest that MoSi2N4 offers robust band edge states that are largely protected from undesirable environmental chemical interactions. The quasiparticle and excitonic properties can then be solely tuned by nonlocal screening effects, which can be achieved through the choice of appropriate substrate and/or controlling the layer thickness. This unique property, together with the moderate band gap and the thermodynamic and mechanical stability, may pave the way for applications of MoSi2N4 in areas including energy, 2D electronics, and optoelectronics.

Methods

First-principles GW and GW-BSE calculations

Structural optimizations for the monolayer, bilayer, and bulk phase MoSi2N4 with different layer stacking patterns are carried out using the Vienna ab initio simulation package (VASP)32,33. The Perdew-Burke-Ernzerhof (PBE)34,35 functional is used for the basic electronic band structure calculations, and the DFT-D3 correction of Grimme et al.18 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 cell of monolayer and bilayer systems.

The layer-dependent quasiparticle and optical properties of MoSi2N4 are calculated using the BerkeleyGW package36 within the GW37 and Bethe-Salpeter equation (BSE)25 approach. Norm-conversing pseudopotentials38 are used for GW calculations, with semicore states (i.e., 4s and 4p) of Mo included as valence electrons. Cutoff for the plane wave expansion of the Kohn-Sham wave functions is set at 125 Ry, and that for the dielectric function is 50 Ry. We use the Hybertsen–Louie generalized plasmon-pole model (HL-GPP)37 to extend the static dielectric function to finite frequencies. For GW calculations of monolayer and bilayer systems, the truncated Coulomb potential39 is used to eliminate the artificial image interactions. We employed two recently implemented accelerated methods21,22, 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 publications21,22. Other computational details have been discussed in Results and Discussion, and also in Supplementary Information.