Stacking order and driving forces in the layered charge density wave phase of 1T-MX2 (M = Nb, Ta and X = S, Se)

The group-V transition metal dichalcogenides (TMD) have attracted a lot of research due to their unique structures and rich physical properties. In these materials, charge density waves (CDW) are still the subject worthy of in-depth research despite being a popular issue. Based on first principles, the stacking effect of the T- MX2 (M = Nb, Ta and X = S, Se) CDW phase are comprehensively explored, with the interplay of correlation effect and magnetic order. Without correlation effect, T- MX2 with one specific stacking order (AA_AC_AA) is most structurally stable and leads to a natural band insulator due to interlayer dimerization. In contrast, same materials with the other stacking orders (AA_AB_AA, AA, AB, AC) are metallic phase. In the presence of correlation effect, whether the systems are insulators or metal highly depends on the magnetic order. AA_AC_AA stacking T- MX2 with antiferromagnetic order end up with band insulator. T- MX2 with AA, AB, AC stacking are metal even in the presence of correlation effect, but with magnetic order, they become Mott insulator. The complication is that AA_AB_AA stacking is located at the intersection of the band insulator and Mott insulator. From this article, we can see T-MX2 with different stacking structure with/without correlation effect and magnetic order show different phases. We outline a rich landscape and systematically explain the causes of the insulating characteristic of the CDW phase and emphasizes the critical role of correlation effect and magnetic order, extending the underlying mechanism of metal-insulation transitions that previously relied only on Mott localization as a driving force.


Introduction
Two-dimensional (2D) materials have opened new horizons for exploring new types of selective electronic states [1][2][3]. Particularly intriguing is the charge density wave (CDW) emergence in the recent research hotspot, which cause a redistribution of the charge density and a periodic spatial modulation, below a certain critical temperature [4][5][6]. Layers held together by weak van der Waals (vdW) forces interact delicately in 2D layered materials, leading to the promising application [7][8][9]. In consequence, the rich physical properties can be modulated by interlayer interactions such as doping, strain, twist or stacking order [10][11][12][13]. Especially in recent years, the unexpected superconductivity, correlated insulating behavior, ferromagnetism, and topological phase found in the magic-angle twisted bilayer graphene have inspired researchers' interests [11,12,[14][15][16].
The group-V transition metal dichalcogenides (TMDs) denoted by MX 2 (M = Nb, Ta and X = S, Se) are among the most studied CDW systems. They all share similar crystal structures with isoelectronic chemical properties [17][18][19][20], and also exhibit the similar commensurate CDW (CCDW) phase distinguished by√by × √bytstar-of-David (SD) clusters [21][22][23]. However, the reason why the CDW phase of 1T-MX 2 is Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
insulating is a subject of intense discussion. The conventional view holds that an unpaired electron exists in the CCDW phase, forming a half-full band, and that the insulating state is driven by the Mott-Hubbard mechanism [24][25][26][27]. Nonetheless, a combination of density functional theory (DFT) calculations [28][29][30], and experimental results [31][32][33][34][35] provides successively strong demonstration that the insulating ground state in the CCDW phase of 1T-TaS 2 is dominated by SD cluster out-of-plane stacking rather than electron-electron correlation [29][30][31][32][33][34][35]. Interestingly, this does not completely exclude the idea that a Mott phase exists in the bulk phase 1T-TaS 2 , because extremely small temperature windows or unpaired SD clusters are where Mott insulating state can be found [36]. Additional research has supported the existence of some stacking randomness along the c-direction [37][38][39]. Furthermore, similar in the CCDW phase of T -MX 2 , they all have a single unpaired electron with S = 1/2 spin moment localizing in the centers of each SD cluster [16,[40][41][42].There is no denying that the stacking pattern is a key factor in determining the origin of the insulating phase. Recently, the behavior of band insulator caused by interlayer dimerization have been also corroborated in the stacking studies on T-phase NbS 2 [16] and T-phase TaSe 2 [43]. However, the comprehensive fundamentals of the specialized stacking order are still debatable, which provides crucial research space left in the whole MX 2 system. Therefore, based on first-principles DFT calculations, in order to better understand the roles that stacking effect, electron-electron correlation, and magnetic order play in the stability and electrical properties of the 1T-MX 2 family (M = Nb, Ta and X = S, Se), this research aims to further investigate the out-of-plane stacking order of the SD clusters in the CCDW phase. It is hoped that these findings will inspire additional theoretical and experimental investigation into the electrical characteristics of the CCDW phase and fill a void in MX 2 research in this area.

Calculation methods
Density functional theory (DFT) calculations were performed by the VASP package [44,45] within the DFT framework [46]. Plane wave basis with an energy cutoff of 520 eV was used to describe the valance electrons, and core electrons were described by the projector-augmented waves (PAWs) potentials. With a grid density of 0.01 Å −1 , the Γ-centered Monkhorst-Pack k-point sampling [47] was used. For the exchange correlation potential, the PBE function of the generalized gradient approximation (GGA) was used [48]. The van der Waals interactions were considered by using the optB88-vdW-DF [49,50]. The geometric positions of atoms were relaxed until the force on each atom is smaller than 0.001 eV Å. The criterion for electronic self-consistent calculation was 10 −6 eV. To deal with the electronic correlation, the GGA + U was applied to the electronic structure calculations, with a Coulomb repulsion of 2.95 eV [51] for the Nb 4d orbital and 2.27eV [52] for the Ta 5d orbital. For monolayer, a vacuum region of ∼20 Å was adopted to minimize the interaction between adjacent slabs.
One SD cluster has a total of 13 atoms, which can be simplified to five sites: A, B, C, B', C', due to the triple rotational symmetry. In order to make the study process more simplified and based on previous studies [28], the stacking sites were concentrated to three: the central atom (site A, red), the second surrounding atom (site B, green) and the outermost atom (site C, blue). The stacking pattern is described as a lateral offset of the central atom with regard to the adjacent layer A site, denoted as AA, AB, and AC, where AA corresponds to the standard stacking pattern. Furthermore, we analyze a two-layer stacking configuration, in which every two layers perform lateral displacement, denoted by AA_AB_AA and AA_AC_AA. The schematic diagram is shown in figure 1.

Results and discussion
The electrical structure of the T-phase MX 2 is similar with a single band across the Fermi energy (E F ) that is predominantly contributed by the M-d orbital as shown in figure 2(b) [30,[53][54][55]. Because of the quasi-twodimensional nature of the layered structure, the band around the E F is almost flat along the Γ-A direction. The sixfold petal-like Fermi surface with an electronic pocket at M exhibits weak dispersion along the K-Z direction, reflecting the quasi-two-dimensional properties of the electronic state. The group-V TMDs all similarly have a significant imaginary frequency along Γ-M that corresponds to the appropriate superlattice vector of structural reconstruction, denoted here as q CDW, and the maximally instable acoustic softened mode around at 0.25-0.33×ΓM [27,55,56], which are successively confirmed typical ´ 13 13 structural reconstruction and transform to a stable CCDW phase with SD clusters [25,27,57,58]. The optimized lattice parameters are shown in table S1 in the SM, which are consistent with experimental [17][18][19][20] and computational [21][22][23] results. Each SD cluster consists of two shells with six M atoms (green and blue) and one central atom (red), as illustrated in figure 2(d). The distorted CCDW phase are energetically more stable than the undistorted structure in terms of lower energy gain as shown in table S1. The phonon dispersion spectra without imaginary frequencies also demonstrate that the CCDW structure is kinetically stable as shown in figure 2(f). The electronic structure of all CCDW phases exhibits a large in-plane band gap and significant out-of-plane dispersion along Γ-A, as shown in figure 2(e). A pair of quasi-1D Fermi sheets are present in the CCDW phase in place of the quasi-2D flowershaped Fermi surface of the high-symmetry phase [21,30,52,59]. The highest occupied band arises dominantly from dz r 2 2 orbital of the central atom and the M atoms at the edges of the clusters barely make a difference [21,23,29,30]. The electronic structure with in-plane flat band characteristics and out-of-plane quasi-1D metallic nature, which indicating very little in-plane overlap between orbitals centered in neighboring clusters, while charge can only hop significantly along the layer stacking direction (c axis).
We calculate the electronic structure within the spin-polarized GGA for the CCDW phase. Relative to the non-spin-polarized results, the energy is only lowered by a tiny value of 2.08 meV/SD for NbS 2 , 1.13 meV/SD for NbSe 2 , 3.12 meV/SD for TaS 2 and 2.11meV/SD for TaSe 2 , respectively within spin-polarized calculations. It also suggests that the electronic correlation effect is not the primary cause of the CDW distortion due to the slight discrepancy. Spin-Orbit Coupling (SOC) also have minor impact on electronic structure as well, as shown in Figure S1(b). When electron correlation effects and SOC are taken into account together, the CCDW phase remains non-magnetic metallic in nature, despite the fact that the in-plane band gap widens due to electron-  electron interactions, as shown in figure S1(c). All suggests that the metallic behavior of CCDW phase may solely result from out-of-plane dispersion, where interlayer coupling effects delocalize in-plane states and inhibit spin polarization, ultimately resulting in a ground state that is nonmagnetic.
Therefore, it can be assumed that electron pairing can lead to the insulating state if the structure of the CCDW phase is doubled along the c-axis. The CCDW phase of the monolayer 1T-MX 2 is a ferromagnetic (FM) Mott insulator, the bilayer structure is an interlayer antiferromagnetic (AFM) insulator [21,25,27,52]. We built a 1 × 1 × 2 superstructure and set the anti-parallel magnetic moment AFM in the adjacent layers as depicted in figure S1(d). It was discovered that a narrow band gap (about 36 meV) can only be opened in steadily increasing U values until U = 3.0 eV, and the magnetic moment about 0.20 ∼ 0.30 μB/atom originates entirely from the central atom of each SD cluster. The influence of the 2D structure is primarily responsible for the ferromagnetic Mott insulating state of T-MX 2 attained in the CCDW monolayer limit [22,25,27,51,60]. Moreover, the interlayer coupling effect and stacking order are not to be disregarded and deserves further consideration.
The computed results, similar in T-MX 2 (M = Nb, Ta and X = S, Se), suggest that the lattice information and stability are significantly influenced by the stacking order, as given in figure 3. E tot , which dictates structural stability and strongly depends on stacking order, is the total energy of the optimized ground state. E vdW is the interlayer vdW action energy, expressed as the energy difference between with and without vdW correction. ΔE vdW is positively proportional with the layer spacing Δd, that is, the vdW energy decreases with decreasing layer spacing [61]. The layer separation of every layer decreases as compared to the AA stacking order, which has a significant impact on the vdW energy. In actuality, the S-S distances vary depending on the stacking configuration, and this spatial influence of the S atoms has a major effect on the vdW energy. Due to shorter S-S distances, interlayer bonding with the p z orbital results in a weaker-bond, while bonding is not feasible for larger S-S lengths [28]. It can be found that although AC stacking is the most favorable for ΔE vdW , it does not dominate in terms of ΔE tot . On the contrary, the AA_AC_AA stacking order has the lowest E tot , but ΔE vdW is not advantageous with respect to AC stacking order. To evaluate the stability of various stacking orders, we defined the average interlayer formation energy in accordance with equation total energy of a specific stacking pattern, E i denotes the energy of the corresponding isolated monolayer, and N denotes the number of layers in the stacking system. The E f can be a metric of stability and can reflect the contribution of microscopic interlayer interactions to the stacking systems [16,62]. In all cases, the AA_AC_AA stacking order is extremely stable from the formation energy standpoint compared to other configurations.
Furthermore, it is evident that the nonmagnetic electrical structure obviously relies on the stacking configuration. Only the AA_AC_AA stacking order opens a natural band gap without the electron correlation effect, as shown in figure 4. The results of NbS 2 are consistent with previous studies [16]. Furthermore, the valence band maximum (VBM) is situated at Γ, which completely different from the Mott insulator mechanism that has traditionally been considered [20,26,63]. Combined with the previous experimental work in which TaS 2 was found to dimerize in the c-axis direction by x-ray diffraction (XRD) and angle-resolved photoelectron emission spectroscopy (ARPES), an educated guess is that the c-directional electron pairing generates a band insulator [20,26,27,37]. Since only one single electron exists for every SD cluster, stability is achieved by pairing electrons in the c-direction. It appears that two lone electrons from M 1L c and M 2L c (L c stands for the central atom in a layer) pair up and fill the isolated orbitals, accompanied by an opening band gap and a decrease in the total energy, the same occurs in M 3L c and M 4L c . Additionally, the central atom contributions of two adjacent layers are nearly identical, demonstrating their degeneracy, as shown in figures 4(e) and (f). It is distinct from the investigations of band alignment [64,65] and evident barriers [66] in heterostructure stacking systems. Furthermore, although we examine the case of homojunction stacking, novelty effects may present in heterojunction stacking systems.
The AA_AC_AA is the most stable stacking configuration, according to the prior computation. The number of electrons in each chalcogen layer and interlayer spacing is indicated in figures 5(d)-(g). As can be seen, significant interlayer dimerization, the interlayer spacing of the AC interface is reduced by 0.36-0.79 compared to AA, on account of the stronger interlayer vdW effect. The conclusion is further supported by the additional electrical accumulation at the AC interface. In MX 2 , different chemical bonds cause the interlayer and intralayer interactions to differ that is why the varying degrees of interlayer dimerization. Common to MX 2 is a valency mismatch between the central M ion (3+) and the neighboring chalcogens ion (2-). There are two primary bonding criteria used to evaluate interlayer dimerization as shown in figure 5(a) [67]: one is incomplete ionic charge transfer denoted as ΔQ I between the transition metal and its neighboring chalcogens, forming M-X ionic bonds. Stronger ionic bonds make the lattice more rigid and less prone to dimerization, theoretically, it cannot be easily deformed. Another is intralayer covalent ΔQ σ , consisting of intralayer σ-electron hopping between the p z orbitals of chalcogens and interlayer σ hopping across the vdW gap, defined as the differential charge at the chalcogen sites upon migration from a monolayer to a multilayer system. A smaller value of |ΔQ σ | indicates that the interlayer transition causes each chalcogen layer to lose less charge, facilitating the interlayer dimerization. A larger value of |ΔQ σ | means that increasing the number of layers results in a softer lattice in the c direction. In comparison to NbS 2 , as the number of layers increased, the weaker Nb-Se bonding and the stronger interlayer interaction of NbSe 2 are represented in the smaller ΔQ I and smaller ΔQ σ , respectively. As a result, interlayer coupling is strengthened, and interlayer dimerization in NbSe 2 is more apparent. Which is consistent with the larger band gap of NbSe 2 , which is around 93.40 MeV in the non-magnetic state ( figure 4(b)), compared to about 45.30 MeV in NbS 2 ( figure 4(a)). TaS 2 has the strongest Ta-S ionic bonding, and the spatial extent of the ionic charge transfer along the c-axis is slightly more expansive (as show figure S2) [67] compared to NbS 2 , making the dimerization effect slightly weaker.
The AA AC AA stacking order displays band-insulating behavior due to paired SD layers. However, the ARPES and x-ray diffraction results of 1T-TaS 2 indicate that the out-of-plane stacking order is rather random and there may be unpaired SD clusters in the c-direction [31,36]. Actually, it is difficult to avoid considering of the intricate phase diagram of the CCDW phase come from stacking disorder due to the slight energy differences between stacking orders. Recent research on 1T-TaS 2 has also demonstrated that, if the continuity of the sequence parameters closely related to the interlayer transition is taken into account, the dimerization cannot be only described by the Peierls mechanism [28,29,31,32,36,68]. In order to investigate the origin of the insulation caused by the AA_AC_AA stacking, we chose three sites, BC 1 , BC 2 and BC 3 , on average, in the BC path of the double layer translation, as shown in the illustration of figure 5. In particular, disregarding the inplane structure and layer spacing, the stacking pattern is the only variable. As shown in Figures S3(b)-(f), it is worth noting that while the other bands do not significantly change throughout this process, the single band near the Fermi level gradually moves downward, eventually drops below the Fermi level at the C site, opening a band gap. The line graph in figure 6 depicts the Γ-point energy of the band structure during the transition from AA_AB_AA to AA_AC_AA stacking configuration. On a slip of approximately 2.40Å from B site, the band gap opens. This indicates that the primary reason of the band structure evolution is lateral slippage in a specific region and the stacking orders are in some way continuous.
The band insulator benefits from the stronger interlayer hopping, but strong correlation still present at the single-layer limit and cannot be absolutely ignored [69,70]. In addition, the spin moment with S = 1/2 in the center of every SD cluster offers a chance to investigate the quantum state and Coulomb correlation in the CCDW phase of T -MX 2 . We set AFM↑↓ and FM↑↑ for AA, AB and AC stacking, AFM↑↓↑↓, AFM↑↓↓↑ and FM↑↑↑↑ for AA_AB_AA and AA_AC_AA stacking, as depicted in figure 7. Different stacking interfaces display different magnetic configurations, and all have lower energy than the non-magnetic (NM) state, which indicates varying degrees of sensitivity to the magnetic order. For NbS 2 and NbSe 2 , the AA and AB interfaces tend to form AFM↑↓ states, so the ground state of AA_AB_AA is AFM↑↓↑↓; while the AA_AC_AA order prefers to form AFM↑↓↓↑ depending on the fact that the ground state of AC interface is FM↑↑. The magnetic states of the AA AC AA stacking order, AFM↑↓↑↓, AFM↑↓↓↑ and FM↑↑↑↑, almost have the same energies, proving that the AA_AC_AA stacking is not susceptible to the magnetic coupling between layers and the characteristics of the band insulator are robust. In the Ta system, the situation is much different that most of the stacking patterns do not exhibit magnetic properties. Various stacking orders offer an effective technique to customize the magnetic interaction between layers, generating fresh suggestions for the search for metamagnetic materials in practical applications [71,72].
When the electron-electron correlation is taken into account under the condition of spin polarization, all other stacking orders transition into insulators, as shown in figures S3(a)-(c). In this case, the MIT is dominated by the Mott-Hubbard mechanism because of the strong electronic correlation between the layers, different from the insulating properties of AA_AC_AA stacking that mainly derived from dimerization and interlayer hopping.  In the AA_AB_AA stacking, it is worth noting that there is a potential band gap that can be open when only considering the electron-electron correlation without spin polarization, as shown in figure S4(d). It has been previously demonstrated that there is no clear demarcation between Mott-insulated and band-insulated phases, the strong electron correlation effect that exists in the interlayer pairing stacking configuration [69]. Notably, the AA_AB_AA stacking order is positioned at the intersection of the Mott insulating state and the band insulating state.

Conclusion and comments
In summary, by means of comprehensive DFT calculations, we investigated the effects of interlayer stacking and electron correlation dependence on the electronic characteristics of the CCDW phase of T-MX 2 . Without taking into account the electronic correlation effects, AA AC AA is the natural and robust band insulator due to its energetically advantageous stacking order. The interlayer dimerization effect in this system can be predicted or explained by the complex chemical bonding. Taking into account the electronic correlation effects, differing stacking configurations result in dramatically diverse electronic structures and magnetic ground states. The AA_AC_AA stacking order transfers from natural band insulator to AFM insulator, while others undergo Slater-Mott metal-insulator transition. Our research emphasizes the critical importance of stacking order and demonstrates the significance of electron-electron correlation on the electrical characteristics of the CCDW phase in T-MX 2 . It also provides a platform to regulate the electronic structure and interlayer magnetic coupling interaction. Additionally, we suggest stacking-controllable metal-insulator transitions and metamagnetic transitions will garner a lot of interest.