Modeling the Superlattice Phase Diagram of Transition Metal Intercalation in Bilayer <i>2H</i>-TaS₂

Van-der-Waals hosts intercalated with transition-metal (TM) ions exhibit a range of magnetic properties strongly influenced by the structural order of the intercalants. However, predictive computational models for the intercalant ordering phase diagram are lacking, complicating experimental pursuits to target key structural phases. Here we use density functional theory (DFT) to construct a pairwise lattice model and Monte Carlo to determine its associated thermodynamic phase diagram. To circumvent the complexities of modeling magnetic effects, we use the diamagnetic ions Zn²⁺ and Sc³⁺ as computationally accessible proxies for divalent and trivalent species of interest (Fe²⁺, Co²⁺, V³⁺ and Cr³⁺), which provide insights into the thermodynamic phase diagram well above the paramagnetic transition temperature. We find that electrostatic coupling between intercalants is almost entirely screened, so the pairwise lattice model represents a coarse-grained charge density reorganization. The resulting phase diagram reveals that entropically favored √3 × √3 ordering and coexisting locally ordered √3 × √3 and 2 × 2 domains persist across a range of temperatures and stoichiometries. This occurs even at quarter-filling of interstitial sites (corresponding to bulk stoichiometries of <i>M</i>_0.25TaS₂; M = intercalant) where there are experimental reports of 2 × 2 and √3 × √3 coexistence but a thermodynamic preference for long-range 2 × 2 order is often assumed. Accompanying Raman spectra revealing quenching-induced √3 × √3 order in <i>2H</i>-V_0.25NbS₂ supports that deficient √3 × √3 domains may dominate 2 × 2 at high temperatures in a range of TM-intercalated transition metal dichalcogenides.