Intercalated architecture of MA2Z4 family layered van der Waals materials with emerging topological, magnetic and superconducting properties

1Monolayer two-dimensional (2D) materials have attracted tremendous interest for their unique electronic properties distinct from corresponding bulk phases, which show promising potential in a variety of fields such as energy storage and conversion,, nanoelectronics,, spintronics,, and superconductivity–. The electronic properties are strongly related to their geometric structures, which can be classified by the number of atomic layers (n).

2Up to now, most of the discovered monolayer 2D materials have a thickness with n ≤ 7. Graphene and hexagonal boron nitride (h-BN) monolayer are famous six-membered ring (SMR) materials of n = 1, which have been a central topic for over a decade–. In the silicon counterpart of graphene, the hexagonal honeycomb tends to be buckled into a double-atomic layer structure, that is, n = 2 (see Fig. 1a). As for n = 3, monolayer transition-metal dichalcogenides (TMDs) are the most studied SMR materials,. In particular, 1-WTe2 monolayer is predicted to be a quantum spin Hall insulator,. Another popular system of n = 3 is the 2D ferromagnetic semiconductor CrI3, of which the magnetism can be manipulated by external electric fields and electrostatic dopings–. Among the family of n = 4, the monolayer group III chalcogenides have accepted growing interest. For instance, the semiconducting monolayer InSe exhibits extraordinary photocatalytic properties. A strong electron–phonon coupling (EPC) can be obtained by an appropriate hole doping for this system, leading to interesting optical and transport properties,. The topological insulator Bi2Se3 and 2D ferromagnetic semiconductor CrGeTe3 are two representative systems of n = 5–. Recently, CaMg was found to have a sextuple layer (n = 6) structure.

3The exploration of monolayer 2D materials has been extended to septuple layer systems in recent years. For instance, a ferromagnetic semiconducting monolayer of MnBi2Te4 (Fig. 1a) can be obtained from its antiferromagnetic bulk phase, which has nontrivial topological properties–. More recently, we experimentally synthesized a new compound of septuple layer SMR material, that is, MoSi2N4 (see Fig. 1a), which has a bandgap of ~1.94 eV with excellent ambient stability. The geometric structures of these materials can be related to other layered systems. For instance, the structure of the MnBi2Te4 monolayer can be viewed as an intercalated structure that is formed by inserting a single layer of the (111) plane MnTe (n = 2) into a Bi2Te3 monolayer (n = 5). Likewise, the MoSi2N4 monolayer can be built by intercalating a 2H-MoS2-type MoN2 layer (n = 3) into an α-InSe-type Si2N2 (n = 4).

4As illustrated in Fig. 1a, we have summarized these known monolayer 2D SMR materials consisting of n = 1, 2, 3, 4, 5, 6, and 7 atomic layer thicknesses. With varying n number, compositions, and constituents, they will become richer in both structures and properties. However, the difficulties lie in how we effectively seek for more monolayer materials with unusual electronic properties. Here, we have proposed a general intercalated architecture approach to systemically construct the structures of MA2Z4 monolayers for the n = 7 family. We illustrate this scheme using density functional theory (DFT) calculations. In addition to reproducing the experimentally synthesized α1-MoSi2N4, α1-WSi2N4, and β5-MnBi2Te4 monolayer materials, we also predict 70, thermodynamically and dynamically, stable MA2Z4 monolayer materials, which exhibit diverse electronic properties including nontrivial topological properties, 2D ferromagnetism, Ising superconductivity, and robust electron valleys.

5We construct the structures of hexagonal MA2Z4 monolayers by intercalating a MoS2-type MZ2 layer into an InSe-type A2Z2 monolayer, that is, the MZ2 layer is sandwiched between two AZ layers. We consider two types of phases for each constituent, that is, the 2H and 1T phases for the hexagonal MZ2 and the α and β phases for the A2Z2 monolayer. Each type of the combinations of MZ2 and A2Z2 involves one pair of the phases. We note that on each side of the MZ2 layer, there are three high-symmetry sites for both A and Z atoms in the AZ layers. Therefore, there are 36 configurations for each combination of MZ2 and A2Z2. As a result, there are 144 possible configurations for each compound. We reduce the number of structures to 39 by removing the duplicated structures based on a symmetry analysis. In addition, we found that, among them, nine configurations are energetically unfavorable. Therefore, we have considered 30 structures (Supplementary Table 1) for our first-principles calculations for each compound, which can be classified into four types, that is, α, β, γ, and δ respectively (see Fig. 1b).

6As a benchmark, we first apply the scheme to MnBi2Te4, MoSi2N4, and WSi2N4. Our calculations found that the β5 phase has the lowest energy for MnBi2Te4, which is exactly the structure of a single MnBi2Te4 layer in the bulk phase–,. For MoSi2N4 and WSi2N4, our results reveal that the α1 phase is the most stable structure, in agreement with our previous experiments. Moreover, we have applied a structure prediction method RG2 to search for possible structures for MoSi2N4 and WSi2N4. As a result, the predicted low-energy structures are within the 30 candidates. These results demonstrate the validity of our strategy in predicting monolayers with a septuple-atomic-layer structure.

7We now investigate the energetics and stability of MA2Z4 monolayers. We have performed calculations for MA2Z4 (M = elements of transition-metal groups IVB, VB, and VIB; A = Si, and Ge; and Z = N, P, and As). There are 54 compounds for this type of material, which have a number of valence electrons ranging from 32 to 34. Moreover, we have also carried out calculations for MA2Z4 monolayers (M = alkali-earth group elements Mg, Ca, and Sr and group IIB elements Zn, Cd, and Hg; A = Al and Ga; and Z = S, Se, and Te). There are 36 compounds in this family, which have a number of 32 valence electrons.

8In Fig. 2a, b, we show the enthalpies of formation of the five lowest-energy structures for each compound, for which the total energies of the α1 and β1 phases are used as the reference, respectively. More information about the enthalpies of formation and structural properties can be found in Supplementary Tables 2–7. We found that the β2 phase has the lowest energy for all the 32-electron systems, with M being the transition-metal group elements. The trends for the systems with 33 and 34 valence electrons are different. For MSi2N4 (M = V, Nb, Ta, Cr, Mo, W), MGe2N4 (M = Nb, Ta, Mo, W), and MSi2P4 (M = Nb, Ta), the α1 phase is the lowest-energy structure, whereas for MSi2P4 (M = Cr, Mo, W), MGe2P4 (M = V, Nb, Ta, Cr, Mo, W), and M(Si/Ge)2As4 (M = V, Nb, Ta, Mo, W), the α2 phase is found to be energetically lower than the others. We further investigate their stability by performing phonon calculations, for which the results are summarized in Supplementary Fig. 1. Our calculations suggest that among the studied systems, with M being the transition-metal group elements, there are 48 compounds dynamically stable. The other six compounds, that is, β2-TiSi2N4, β2-TiGe2N4, β2-VGe2N4, β2-CrGe2N4, β1-CrSi2As4, and CrGe2As4, may be dynamically unstable since they exhibit imaginary phonon branches.

9The monolayers based on the alkali-earth group and group IIB elements are in either the β1 or β2 phase. Twenty-four compounds are expected to be stable based on the phonon calculations (Supplementary Fig. 2). It should be mentioned that four of them, that is, β1-ZnAl2S4, β1-ZnGa2S4, β1-MgAl2S4, and β1-MgAl2Se4, can be found from the 2DMatPedia database, which was obtained by the high-throughput computational method. This agreement again indicates the validity of our intercalation method in predicting structures of the MA2Z4 monolayers. Moreover, we have performed ab initio molecular dynamics simulations for selected systems, that is, β2-SrGa2Te4, δ4-VSi2P4, α1-TaSi2N4, and α2-WSi2P4 (Supplementary Fig. 3), whose electronic properties will be discussed below. Our results indicate that they are thermally stable. In all, we have predicted 72 new monolayer materials with a septuple-atomic-layer structure.

10The intercalation method leads us to find a series of septuple layer MA2Z4 materials with diverse electronic properties, summarized in Supplementary Figs. 4–7, which are distinctly different from those of the constituents MZ2 and A2Z2. For instance, the 2H-NbN2 monolayer is a metallic ferromagnet, while α-Si2N2 is a semiconductor. However, the intercalated system α1-NbSi2N4 is surprisingly a ferromagnetic semiconductor (Supplementary Fig. 5). We found that the electronic properties of MA2Z4 can be classified according to the total number of valence electrons of these systems. The MA2Z4 nitrides with 32 valence electrons are semiconductors. However, their phosphides and arsenides with such a number of valence electrons are metallic, which is due to a much weaker band hybridization between the M and Z atoms in these systems than in the nitrides. Most of the MA2Z4 systems with 34 valence electrons are semiconductors, except for CrGe2N4, CrSi2As4, and CrGe2As4, which are ferromagnetic metals. In all, we found 21 semiconductors among the families of 32 and 34 valence electrons. We list the sizes of the bandgaps of all the semiconducting monolayers in Supplementary Table 2, which are from both the standard DFT and HSE hybrid functional (HSE06) calculations.

11Besides, we have obtained many unexpected electronic properties for MA2Z4 monolayers, such as nontrivial topological properties, 2D ferromagnetism, Ising superconductivity, and robust electron valleys with spin-momentum locking. In the following, we will discuss them in more detail.

12We first demonstrate that some of these systems exhibit nontrivial topological properties. We pay our attention to the family of MA2Z4 (M = Mg, Ca, Sr, Zn, Cd, and Hg; A = Al and Ga; and Z = S, Se, and Te). Most of them are semiconductors. The semiconducting nature can be understood based on the structural properties and the chemical bonding in this type of system. In the InSe-like structure of the A2Z2 monolayer, there is a covalent bond between the two A atoms. However, this bond is broken upon intercalating an MZ2 monolayer into the A2Z2 unit to form the intercalated system MA2Z4. Consequently, each atom tends to be charge compensated, which makes this system semiconducting. We have predicted a series of topological insulators in this family on the Perdew–Burke–Ernzerhof (PBE) level (Supplementary Table 3). For instance, our calculations show that β2-SrGa2Se4 is a zero-gap semiconductor (Supplementary Fig. 8). Inclusion of spin–orbit coupling (SOC) leads to a gap opening of ~68 meV and band inversion at Γ. Our calculations see gapless states in this system and found that it has a value of Z2 = 1. These results suggest that β2-SrGa2Se4 is a topological insulator. We further investigate its topological properties by using the advanced HSE hybrid functional method. As a result, it becomes topologically trivial.

13In contrast, we found that β2-SrGa2Te4 is topologically nontrivial by both the PBE and HSE06 methods. In Fig. 3a, we show the electronic bands from PBE calculations without SOC, which indicates that this system has a metallic band structure. The inclusion of SOC induces a gap opening near Γ, which results in a separation of the valence and conduction bands (see Fig. 3b). Similar to the case of β2-SrGa2Se4, we have Z2 = 1 and gapless edge states. We then investigate the topological property by performing HSE06 calculations without and with SOC, of which the results are shown in Fig. 3c, d, respectively. The calculation without SOC finds that there is a small direct bandgap at Γ. The inclusion of SOC induces a downward bending of the valence band at Γ, which signals the band inversion seen in the topological insulator Bi2Se3. We further confirm this band inversion by the orbital-projected band structure shown in Fig. 3d. We have also calculated the evolution of the Wannier charge centers, which also gives Z2 = 1 (see Fig. 3e). Moreover, we have observed gapless edge states for β2-SrGa2Te4 (see Fig. 3f) based on the Wannier functions obtained from the HSE06 calculations, which again confirms that this system is topologically nontrivial.

14These systems are in either the ferromagnetic state or the nonmagnetic state in this family. Among them, there are nine ferromagnetic systems, for which the magnetic configurations and the energetics are shown in Supplementary Fig. 9 and Supplementary Table 9, respectively. Our standard PBE calculations found that among them, NbSi2N4 and MGe2N4 (M = V, Nb, and Ta) monolayers show a half-metallic behavior (Supplementary Fig. 5), whereas VSi2Z4 (Z = N, P, and As) and VGe2Z4 (Z = P and As) exhibit a gapless semiconducting behavior,.

15A close inspection of the PBE-derived band structure finds that the Fermi level of VSi2P4 touches both the valence band maximum (VBM) and conduction band minimum (CBM), which are the spin majority and minority states, respectively (Fig. 4a). In addition, we have calculated the magnetocrystalline anisotropy energy of VSi2P4, for which it is 56 μeV. This result indicates that its easy axis is out of plane. Furthermore, our HSE06 calculation reveals that it is nearly a direct bandgap semiconductor (Fig. 4b). We realized that the electron correlations related to the partially filled 3d orbitals of the transition-metal atoms may be important for the band structure of such a kind of system. We then estimate the effective Coulomb interaction (U) by using the linear-response method, which gives U = 4 eV for the V-3d orbitals in VSi2P4. In Fig. 4c, we show the band structure of VSi2P4 from our PBE calculation with U (PBE + U). One can see that inclusion of the Hubbard U opens a gap of ~0.3 eV, which is then increased up to ~1.0 eV by the HSE06 calculation (Fig. 4d). Our calculations show that the HSE06 results are consistent with the PBE + U results that both predict a semiconducting behavior for the VSi2P4 monolayer. Moreover, we make a comparison of the band structures for the HSE06 calculations with and without Hubbard U and found that the sizes of the predicted bandgaps are comparable, although there are visible differences in details of the electronic bands near K. Based on these observations, we then performed HSE06 calculations for all the ferromagnetic systems. As a result, they all show a sizable gap ranging from 0.2 to ~1.0 eV (Supplementary Fig. 5).

16We use the Heisenberg spin Hamiltonian to understand the magnetic interactions in the VSi2P4 monolayer. We extract the exchange interaction parameters between spins by fitting the total energies from our DFT calculations of various spin configurations to the Hamiltonian with the first nearest neighbors only, that is, H = − J∑<ij>Si⋅Sj. We found that J is ~8 meV for S = 1. In fact, J2 is negligibly small (<0.1 meV) if we fit the Hamiltonian up to the second nearest neighbors. Based on the magnetic exchange parameters, we obtain a Curie temperature of ~90 K by performing Metropolis Monte Carlo simulations of the Heisenberg Hamiltonian (Supplementary Fig. 11).

17The nonmagnetic systems of this family are metallic since they have one unpaired electron, which results in a half-filled electronic band (Supplementary Fig. 5). We note that among the systems, α1-TaSi2N4 has a special band structure, that is, it has a disentangled band contributed by the Ta-5d orbitals. Such a feature leads to the speculation of strong electron correlation in this system. Like in the case of VSi2P4, we have also carried out a linear-response calculation, which gives a value of ~2.5 eV for the U parameter of the Ta-5d orbitals. Then, we investigate the effects of Hubbard U on the band structure of TaSi2N4 and found that the half-filled band remains almost unchanged (Supplementary Fig. 12). Note that, in this system, the strong SOC along with the inversion symmetry breaking induces a large valley-contrasting spin splitting at K and K′. Such a band structure favors the type-I Ising superconductivity, as already observed in NbSe2. Following this inspiration, we have calculated the phonon spectrum and Eliashberg function (α2F(ω)) as well as derived the EPC strength (λ) for α1-TaSi2N4 (Fig. 5). Using λ = 0.66 and the calculated logarithmic average phonon frequency of 305.58 cm−1, we have derived a superconducting transition temperature of Tc = 5.42–13.61 K via the Dynes modified McMillan formula with the effectively screened Coulomb repulsion constant of μ from 0.15 to 0.05. For α2-TaGe2P4, Tc is ~3.75 K for μ = 0.1 (Supplementary Fig. 13). Moreover, we have examined the effects of SOC and strain on λ and Tc of α1-TaSi2N4. We found that SOC has a minor effect on both λ and Tc (see Fig. 5g, e and Supplementary Fig. 13). We summarize our results of the strain effect (−3 to 6%) in Fig. 5e, which show that Tc varies nonlinearly with strain. Tc is only reduced by ~32% under a large strain up to 6%. While Tc is increased by ~24% under a compressive strain −3%. This trend indicates that the superconductivity in α1-TaSi2N4 is robust against strain, which favors the observation of Ising superconductivity in the epitaxially grown monolayer.

18Among this family, we focus on Mo- and W-based MA2Z4 (A = Si and Ge and Z = N, P, and As) monolayers. The reason is that some of them even show better electronic properties than the dichalcogenide counterparts of Mo and W, that is, the TMDs monolayers. For the nitrides, that is, MSi2N4 and MGe2N4 (M = Mo and W), the lowest-energy structures are in the α1 phase, whereas for the phosphides and arsenides, the α2 structure becomes favorable. Despite this difference, the inversion symmetry is absent in both structures. Among these systems, MoSi2P4, MoSi2As4, WSi2P4, and WSi2As4 are predicted to be direct bandgap semiconductors by both the PBE and HSE06 calculations. One prominent feature of the band structure is that like the TMDs monolayers, there are electron valleys at K and. The strong atomic SOC in Mo and W induces large spin splittings at K and. For α2-WSi2P4 (Fig. 6a), the spin splitting at K is ~0.4 eV, which is comparable to that of 2H-WSe2 monolayer,. In addition, the strong SOC together with the inversion symmetry breaking gives rise to spin-momentum locking at the valleys and a Berry curvature contrasting at the two valleys (see Fig. 6c). These features provide opportunities for exploring the spin-valley physics and manipulating their electronic properties via doping, interfacing, and layer stacking including Morie patterns,.

19The unique structural and electronic structures of our systems show advantages over the TMDs layers in several aspects. First, the electron valleys in them are robust against interlayer interactions. In the case of the TMDs bilayers, by the layer hybridization, the conduction band at Γ is pushed up to a higher energy and the valence band at kc is pushed down to a lower energy than that at K and (see Fig. 6e and Supplementary Fig. 15), resulting in a direct–indirect bandgap transition. Note that there is a large energy difference between the VBM at K and Γ for α2-WSi2P4. Moreover, the valence band and conduction band are dominated by the d-orbitals of W (Mo) (Supplementary Fig. 16), which is sandwiched by double-atomic AZ layers. Therefore, the VBM and CBM can be less affected by interfacing and layer stacking than those of the TMDs layers. In Fig. 6b, we show the band structure of the α2-WSi2P4 bilayer. One can see that the nature of the direct bandgap is maintained well upon the layer stacking. In addition, we found that α2-WSi2P4 exhibits a large hole mobility of up to 460 cm2 V−1 s−1 and an electron mobility of about 150 cm2 V−1 s−1 (see Supplementary Fig. 17 and Supplementary Table 8). These values are ~150% of those for the 2H-WSe2 monolayer,, indicating that α2-WSi2P4 may have better electronic transport properties than the TMDs monolayers.

20In summary, we have presented an intercalation approach to construct septuple-atomic-layer MA2Z4 monolayers. We have illustrated this strategy by performing first-principles calculations for systems with 32–34 valence electrons. Our calculations predict that among 90 candidates, 72 compounds are both thermodynamically and dynamically stable. The systems with 32 and 34 valence electrons are mostly semiconductors, whereas those with 33 valence electrons are either nonmagnetic metals or ferromagnetic semiconductors. In addition, we found that these systems exhibit a number of novel electronic properties. In particular, β2-SrGa2Te4 is found to be a topological insulator by both the PBE and HSE06 calculations. In addition, our study finds that among the family of the systems with 33 valence electrons, the V-based MA2Z4 monolayers are basically ferromagnetic semiconductors as revealed by our HSE06 calculations. As for the nonmagnetic metals like α1-TaSi2N4 in this family, the inversion symmetry breaking and the strong SOC lead to a valley-contrast spin splitting in the half-filled band with out-of-plane spin polarizations. We further predict that such features favor the type-I Ising superconductivity in α1-TaSi2N4. Moreover, our study reveals that α1-WSi2P4 not only has a direct bandgap with robust electron valleys against layer interactions but also shows superior electronic transport properties compared to the TMDs monolayer WSe2. Finally, we would like to emphasize that our currently proposed intercalated architecture approach can be indeed extended to MA2Z4 monolayer materials with M for late transition metal elements, such as MnBi2Te4 for which our current calculations also correctly capture the agreements to experiments. Furthermore, it can be generalized to a wider way. For instance, n=7 MA2Z4 monolayer materials can be constructed by intercalating n = 2 silicene-like monolayer into n = 5 Bi2Se3-like monolayer and, we can also even combine n = 3 MoS2-like monolayer and n = 5 Bi2Se3-like monolaye to form new type n = 8 monolayer materials, and so on. Our study thus provides an effective way of designing septuple layer systems with unusual electronic properties.

21First-principles calculations were performed by applying the Vienna ab initio simulation package (VASP),, and the PBE exchange-correlation functional was used to calculate the enthalpy of formation and band structure. To further get a more accurate bandgap, a hybrid function (HSE06) was applied. The plane-wave energy cutoff was adopted to be ranging from 250 eV (MgGa2Te4) to 500 eV (MoSi2N4). The 2D Brillouin zone was sampled by a 15 × 15 k-mesh for the self-consistent calculations. Structural relaxations were done with a threshold of 10−3 eV Å−1 for the residual force on each atom. The energy convergence criteria were set as 10−6 eV. To minimize the interactions between the monolayer and its periodic images, a vacuum of 20 Å was used for all the calculations.

22Phonon dispersions were obtained using density functional perturbation theory and Phonopy package. We used a large supercell (either a 4 × 4 supercell or a 5 × 5 supercell) for the force calculations.

23For metallic materials, the EPC constant λ is given bywhere α2F(ω) is the Eliashberg function, which is defined aswhere N(ϵF) is the density of states at the Fermi level, ωqν is the phonon frequency of the mode ν at wavevector q and γqν is the phonon linewidth or lifetime.

24The Eliashberg function α2F(ω) is used to calculate logarithmic average phonon frequencies by. We use the Allen–Dynes-modified McMillan formula to estimate the superconducting transition temperature, that is,.

25The above calculations were performed using the Quantum-ESPRESSO package with ultrasoft pseudopotentials and local-density approximation exchange-correlation functional,. For α1-TaSi2N4, the kinetic energy cutoff and charge density cutoff of the plane-wave basis are 60 and 480 Ry, respectively. A 32 × 32 k-mesh with Marzari–Vanderbilt cold smearing of 0.02 Ry is used for self-consistent calculations. A 4 × 4 q-mesh is used to obtain dynamic matrix and EPC constant, respectively. For α2-TaGe2P4, a kinetic energy cutoff of 80 Ry and a charge density cutoff of 640 Ry are used. A 36 × 36 k-mesh and 6 × 6 q-mesh are used to calculate EPC constant and the superconducting Tc.

26The intrinsic carrier mobility μ of 2D materials was derived based on the deformation potential (DP) approximationwhere C is the elastic modulus defined as, m* is the effective mass at the CBM or VBM, and T is the temperature. Here, room temperature T = 300 K was used. E1 is the DP constant defined as, where ΔE is the change of the eigenvalue at CBM or VBM and Δl is the lattice dilation along deformation direction.

27The enthalpy of formation of a MA2Z4 monolayer (per atom) can be expressed as:where Etot is the total energy of the system, and EM, EA, and EZ are the ground state total energies of the elementary crystals of M, A, and Z, respectively.

28The Berry curvature of a 2D material with n bands can be defined as–:where un,k is the lattice periodic part of the Bloch wave functions. Our calculations were performed using the Wannier90 package, which constructs un,k via an ab initio tight-binding method on the basis of maximally localized Wannier functions.

29The values of the topological invariant Z2 were obtained by calculating the Wannier charge centers. The edge states were obtained using an iterative Green functions method. The above calculations were based on the ab initio tight-binding parameters from the Wannier90 calculations.