A first principles phase diagram calculation, that included van der Waals interactions, was performed for the 3D bulk system $(1\u2212X)\xb7MoS2\u2212(X)\xb7MoTe2$. Surprisingly, the predicted phase diagram has at least two ordered phases, at $X\u22480.46$, even though all calculated formation energies are positive; in a ground-state analysis that examined all configurations with 16 or fewer anion sites. The lower-temperature *I*-phase is predicted to transform to a higher-temperature $I\u2032$-phase at $T\u2248500\u2009K$, and $I\u2032$ disorders at $T\u2248730\u2009K$. Both these transitions are predicted to be first-order, and there are broad two-phase fields on both sides of the ordered regions. Both the *I*- and $I\u2032$-phases are predicted to be *incommensurate*, i.e., *aperiodic*: *I*-phase in three dimensions; and $I\u2032$-phase in two dimensions.

## I. INTRODUCTION

Recently there has been great interest in two-dimensional (2D) transition metal dichalcogenide (TMD) materials such as *MoS*_{2}, *MoSe*_{2}, and *MoTe*_{2}, their solid solutions, and related 2D materials.^{1,2} Traditionally, *MoS*_{2} has been used as a dry lubricant^{3} that is stable up to 623 K. Currently, interest is focused on applications as: band-gap engineering materials;^{4,5} nano-electronic devices;^{2,6–8} photovoltaic devices;^{9,10} valleytronics applications;^{11,12} 2D building blocks for electronic heterostructures;^{13} and as sensor materials.^{14}

The individual, three-atom-thick, 2D-layers of the bulk system are bonded by van der Waals forces, hence these forces influence the bulk and multilayer synthesis and therefore anion order-disorder and/or phase separation in solid solutions. The results presented below, for 3D bulk $MoS2\u2212MoTe2$, imply that van der Waals interactions may strongly affect phase stabilities, either between adjacent layers in bulk or few-layer samples, or between monolayers and heterogeneous substrates.

Of the three quasibinary solid solutions ($MoS2\u2212MoSe2$, $MoSe2\u2212MoTe2,\u2009MoS2\u2212MoTe2$) $(1\u2212X)\xb7MoS2\u2212(X)\xb7MoTe2$ has the greatest difference in anionic radii (*R _{S}* = 1.84 Å;

*R*= 2.21 Å),

_{Te}^{15}which suggests that it is the most likely to exhibit interesting solution behavior. One expects a simple miscibility gap as reported by Kang

*et al*.

^{4}for monolayer $MoS2\u2212MoTe2$, hence the prediction of two configurational entropy (

*S*) stabilized

_{con}*incommensurate*, i.e.,

*aperiodic*, phases is extraordinary (stable phases that have positive formation energies must be entropy stabilized).

## II. METHODOLOGY

### A. Total energy calculations

Total structure energies, $\Delta EStr$ were calculated for fully relaxed *MoS*_{2}, *MoTe*_{2} (2H-structure, space group $P63/mmc$, AB-stacking of three-atom-thick layers), and for 233 $Mom+n$(S_{m}Te_{n})_{2} supercells. The Vienna $ab\u2009initio$ simulation program (VASP, version 5.3.3^{16,17}) was used for all density-functional theory (DFT) calculations, with projector augmented waves (PAW) and a generalized gradient approximation (GGA) for exchange energies. Electronic degrees of freedom were optimized with a conjugate gradient algorithm. Valence electron configurations were: Mo_pv $\u20094p55s4d$; S_$s2p4$; Te_$s2p4$. Van der Waals interactions that bond the three-atom thick 2D X-Mo-X layers (X = S, Te) together were modeled with the non-local correlation functional of Klimes *et al*.^{18} Total energies were also calculated *without* van der Waals interactions, but, up to a basis of 140 structures, ground-state analyses always predicted false ground-states (supplementary material). Convergence with respect to k-point meshes was achieved by increasing the number of k-points until the total energy converged. A 500 eV cutoff-energy was used in the “high precision” option, which converges *absolute* energies to within a few meV/mol (a few tenths of a kJ/mol of exchangeable S- and Te-anions). Precision is at least an order of magnitude better. Residual forces of order 0.02 eV or less were typical. Often, convergence with respect to hexagonal c-axis length was not automatic, and it was necessary to choose an initial c-axis value that is close to the converged value. Calculated interlayer spacings in *MoS*_{2} and *MoTe*_{2} are 2.992 Å and 3.513 Å, respectively, corresponding experimental values are: 2.977 Å (Ref. 19) and 3.382 Å.^{20}

Formation energies ($\Delta Ef$) for 233 Mo_{m+n}(S_{m}Te_{n})_{2} supercells are plotted in Fig. 1, in which values for $\Delta Ef$ are normalized per mol of exchangeable anions, S and Te

Here *E _{Str}* is the total energy of the Mo

_{m+n}(S

_{m}Te

_{n})

_{2}supercell; $EMoS2$ is the energy/mol of MoS

_{2}; $EMoTe2$ is the energy/mol of MoTe

_{2}.

All supercell energies are positive which suggests a miscibility gap system, unless one or more entropy stabilized phases are stable.

### B. The cluster expansion Hamiltonian

A cluster expansion Hamiltonian (CEH),^{21} for the (1-X)$\xb7$ MoS_{2}–(X)$\xb7$ MoTe_{2} quasibinary system was fit to the set of 235 formation energies, $\Delta EVASP$, solid dots (green online) in Fig. 1 with a cross validation score of (CV)^{2} = 0.00723896. Fitting of the CEH was performed with the Alloy Theoretic Automated Toolkit (ATAT)^{17,22–25} which automates most of the tasks associated with CEH construction. A complete description of the algorithms underlying the code can be found in Ref. 23. Large open squares in Fig. 1 (red online) indicate values of the 235 $\Delta EFit$ that were calculated with the CEH. Smaller open squares ($\Delta EGS$, blue online) indicate the results of a ground-state analysis in which the CE was used to calculate the formation energies for all ordered configurations with 16 or fewer anion sites, 151 023 structures.

## III. RESULTS AND DISCUSSION

A first principles phase diagram (FPPD) calculation was performed with grand-canonical, and canonical, Monte Carlo (MC) simulations using the emc2 and phb codes which are part of the ATAT package.^{22–24} Most phase boundaries were calculated with the phb program which uses equilibration tests to set the numbers of equilibration- and MC-passes.^{24} To draw high-T extensions of the two-phase fields, and to locate the $I\u21ccI\u2032$ transition, a 48 × 48 × 12 unit cell simulation box was used, with 2000 equilibration passes and 2000 MC-passes (see supplementary material for comparisons of various equilibration- and MC-pass settings in calculations of the $I\u21ccI\u2032$ phase transition). The predicted phase diagram is shown in Fig. 2; where $(MoS2)$ denotes the *MoS*_{2}-rich solution phase, and similarly for $(MoTe2)$.

Kang *et al*.^{4} performed first principles phase diagram calculations for four dichalcogenide monolayer systems: $MoSe2(1\u2212x)Te2x,\u2009WSe2(1\u2212x)Te2x$, $MoS2(1\u2212x)Te2x$, and $WS2(1\u2212x)Te2x$. Their CEs were fit to about 40 structures per system, and van der Waals interactions (with substrate) omitted. All systems were predicted to have miscibility gaps, and surprisingly, all consolute points are on the Te-rich sides. One expects the consolute point to be on the S-rich side, because it typically requires less energy to substitute a smaller S-ion into a larger Te-ion site, than vice versa. Figure 2 also has reduced solubility on the Te-rich side, but this is related to immiscibility between the **I**- ($I\u2032$)-phase, and the Te-rich phase.

It is not clear that the Kang *et al*. $MoS2\u2212MoTe2$ phase diagram would still be a simple miscibility gap had they included hundreds of structures in their CE-fit rather than about 40. Hence the monolayer vs. bulk comparison is uncertain.

Surprisingly, the phase diagram predicted here has multiple two-phase fields, separated by two ordered incommensurate phases, neither of which is a ground state. To seven digits, the calculated bulk composition of the **I**-phase, just above its 82.5 K minimum temperature of stability, is X = 0.4642857 = 13/28; i.e., $Mo14S15Te13$. Stability of the **I**-phases is a robust result, in Monte-Carlo simulations: (1) CEH fits to 128, 153, 162, 182, 225, and 235 formation energies all predict **I**-type ordering; (2) **I**-phase forms spontaneously on cooling of the $I\u2032$-phase at $T\u2272500\u2009K$, and on heating of a low-T equilibrium $(MoS2)+(MoTe2)$ assemblage to $T\u227382\u2009K$; (3) I$\u2009\u2032$-phase forms spontaneously on heating of the **I**-phase, or cooling of a disordered solid solution with $X\u22480.46$. This calculation considers only *S _{con}*, and ignores excess vibrational entropy,

*S*, which could conceivably destabilize the

_{vib}**I**-phases. In light of (1) above, however, this seems highly unlikely. Also, there is no fully satisfactory way to model

*S*for an aperiodic phase, and a reasonable approximate structure (with

_{vib}**I**-phase like ordering) would require at least a low symmetry 84-atom cell; which is beyond the scope of this study.

Figure 3 shows how the CE-calculated *I*-phase formation energy $\Delta Ef(I\u2212phase)\u2261\Delta EI$, varies as a function of MC-supercell size and shape. The flat minimum at $\Delta EI\u22720.03$ eV/anion is the calculated *I*-phase formation energy which is plotted as the solid black diamond in Fig. 1. Supercell dimensions were chosen to accommodate $Mo14S15Te13$ stoichiometry and **I**-phase ordering. Note that many of the $\Delta Ef$ plotted in Fig. 1, are lower in energy than $\Delta EI$, but that they are for *periodic* structures with 16 or fewer anion sites in which $Scon\u21920$ as $T\u21920\u2009K$, and clearly (Fig. 4) the *I*-phase is incommensurate, i.e., *aperiodic* with $Scon>$ 0. Figure 4 exhibits the S:Te (yellow:brown online, respectively) ordering at: (a) 200 K; and (b) 575 K, i.e., below and above the $I\u21ccI\u2032$ phase transition. Mo-atoms are omitted for clarity, and labels (001)_{D}, (100)_{D}, and (010)_{D} refer to corresponding crystallographic planes in the high-T $P63/mmc$ *disordered*-phase.

In the (001)_{D}- and (100)_{D}-planes,. $SmTen$. chains in the ⟨010$\u27e9D$ direction, most often have $m=5\u2009or\u20096$ and $n=4\u2009or\u20095$. Also, in (100)_{D}, the $\u2026SmTen\u2026$ chains exhibit irregular alignments relative to one another. Note however, that $SmTen$-chains in (001)_{D}-planes *order* along the ⟨001$\u27e9D$ direction, such that *S _{m}*-units alternate with

*Te*-units in adjacent 3-atom thick 2D-layers; i.e., …$SmTenSm\u2032Ten\u2032$… chains are stacked on top of …$TenSmTen\u2032Sm\u2032$… chains with inescapable misfits, owing to the different and variable values of

_{n}*m*and

*n*. Thus

*I*-phase ordering is inevitably imperfect, aperiodic, and incommensurate, which suggests that the $I\u21ccI\u2032$ phase transition is first-order.

The difference between *I*- and $I\u2032$-phases appears to be a distinction between 3D-ordering in the low-T *I*-phase and 2D-ordering in the high-T $I\u2032$-phase. Clearly, ordering is stronger in the *basal* (001)_{D}-plane than in the (100)_{D}- or (010)_{D}-planes, and striped order within (001)_{D} persists above the $I\u21ccI\u2032$ transition.

Figure 5(a) is a Monte-Carlo T-scan (heating) of the total energy $ETOT(T)$ which confirms first-order character for the $I\u21ccI\u2032$ transition: a critical (continuous) transition^{26} would not exhibit the sharp change at $T\u2248503\u2009K$; also a transition from the lower-T, higher-*E _{TOT}* phase to the higher-T, lower-

*E*phase requires that the lower-T phase be superheated, i.e.,

_{TOT}*metastable*before it transforms. Figure 5(b) is an idealized schematic that compares equilibrium- and metastable-transition paths (solid black line, and blue line with red dots, online respectively). This transition is more subtle in cooling simulations, but still evident in snapshots. In the Fig. 5 inset, no change of slope at the transition is evident in the Helmholtz energy, F(T), which suggests that the transition is

*weakly*first-order.

^{27}Also, when the MC temperature-increment was decreased from 1.0 K/MC-step to 0.1 K/MC-step (not shown), the predicted transition temperature decreased by about 10 K, which indicates more superheating at 1.0 K/MC-step than at 0.1 K/MC-step, hence the first-order character. This transition is shown as a dotted line in Fig. 2 because the two-phase fields that a first-order transition implies are too narrow to resolve in the MC-simulations.

## IV. CONCLUSIONS

To summarize, a first principles phase diagram calculation for the 3D bulk system $(1\u2212X)\xb7MoS2\u2212(X)\xb7MoTe2$, that includes van der Waals interactions, predicts the formation of two entropy stabilized *incommensurate*, i.e., *aperiodic* phases: the *I*- and $I\u2032$-phases, at X $\u22480.46$. Above the minimum temperature for stability of the *I*-phase, $T\u224882\u2009K$, the calculation predicts broad two-phase fields between the *I*- or $I\u2032$-phase and disordered S- or Te-rich solution phases, (*MoS*_{2}) and (*MoTe*_{2}), respectively. Both the $I\u21ccI\u2032$ and $I\u2032\u21ccdisordered$ transitions are predicted to be first-order. Dramatic changes in phase relations can be induced by arbitrarily small differences in energy, hence van der Waals interactions should not be ignored in layered 2D-systems such as $MoS2\u2212MoTe2$.

## SUPPLEMENTARY MATERIAL

See supplementary material for the ground state analysis of cluster expansions based on energies estimated without including van der Waals interactions, calculated transition points obtained from MC simulations for the high-T extensions of the two-phase fields and the comparisons of various equilibration- and MC-pass settings in calculations of the $I\u21ccI\u2032$ phase transition.

## ACKNOWLEDGMENTS

This work was supported by NIST-MGI, and B. Burton wishes to thank R. Selinger, J. Douglas, K. Migler, A. Davydov, B. Campbell, J. Lau, and L. Bendersky for useful discussions.

## References

The phrase “weakly first-order” is used to describe a transition with a sufficiently small discontinuity in order parameter that some first-order characteristics are undetectable: e.g., a kink in F(T), or **I** + **I**′ two-phase fields.