High-throughput characterization of transition metal dichalcogenide alloys: Thermodynamic stability and electronic band alignment

Alloying offers a way to tune many of the properties of the transition metal dichalcogenide (TMD) monolayers. While these systems in many cases have been thoroughly investigated previously, the fundamental understanding of critical temperatures, phase diagrams and band edge alignment is still incomplete. Based on first principles calculations and alloy cluster expansions we compute the phase diagrams 72 TMD monolayer alloys and classify the mixing behavior. We show that ordered phases in general are absent at room temperature but that there exists some alloys, which have a stable Janus phase at room temperature. Furthermore, for a subset of these alloys, we quantify the band edge bowing and show that the band edge positions for the mixing alloys can be continuously tuned in the range set by the boundary phases.


I. INTRODUCTION
Monolayer TMDs constitute a class of two-dimensional (2D) materials, commonly of MX 2 stoichiometry, where M is a transition metal (e.g., Mo, W, Zr, Hf, Ti, Pd, Pt) and X is a chalcogen (S, Se, Te). These monolayer compounds exhibit a structure where the transition metal is sandwiched between chalcogen atoms and exhibit hexagonal or trigonal symmetry (Fig. 1a). TMDs have received a considerable amount of attention over the last decade, e.g., due to the excellent optical properties of the semiconducting group 6 TMD MoS 2 [1], the prospects of vertical integration into heterostructures [2], the emergent properties of moiré structures [3,4], and the outlook of using TMDs as catalysts for the hydrogen evolution reaction [5]. Many TMD monolayers exhibit moderate band gaps of around 1 to 3 eV, which is suitable for photochemical applications or for use as channel materials in nanoelectronics [6][7][8]. Due to the emergence of high-κ dielectrics as gate material in transistors such as ZrO 2 and HfO 2 , 2D TMDs based on the same transition metals may also play a role in future high-κ transistor designs [9].
The electrical transport properties as well as the optical response is ultimately dependent on band edge positions. The position of the band edges are particularly important in devices based on heterojunctions. Specifically, TMD heterostructures with type-II alignment can be used for photoseparation of charge carriers [10]. There are various ways to tailor the electronic properties, including functionalization [11], Coulomb engineering [12][13][14], strain engineering [15,16], and alloying [14,17,18]. While alloying may allow for a continuous tuning of the band edge position, not all alloys can be manufactured due to their inherent thermodynamical properties. The main objectives of this work are therefore to map out the phase diagrams of a large number of binary TMD alloys and the band edge variations of promising alloy combinations. * erhart@chalmers.se All of the binary alloys based on (Mo, W)(S, Se) 2 have been synthesized [19], usually in the hexagonal H structure (spacegroup P6m2, ITCA number 187; Fig. 1b), and have been thoroughly studied in the context of band gap engineering by composition control [18,[20][21][22][23][24][25]. It has been shown that these alloys are random at ambient conditions due to (very) small mixing energies [18]. Although ordered phases have been found in computational studies [19], these are not observed in experiments [26,27], most likely due to the very small energy gains on the order of 1 meV compared to the respective disordered solution. In addition, the MoS 2x Te 2(1x ) alloy has been fabricated [28]. For the optical properties, in particular, it is relevant whether alloying preserves the valley degree of freedom found in the boundary phases, i.e., the nonalloyed systems. In this regard, it is noteworthy that, e.g., Mo x W 1x S 2 has been recently shown to be robust in this regard [27].
While the group 6 TMD alloys have been extensively studied, the prospect of alloying these TMDs with transition metals from other groups of the periodic table has not been systematically addressed, although doping of e.g., MoS 2 with Ti has been reported [29]. In addition, there are various families of TMDs, including different structural prototypes, for which the alloying behaviour is largely unexplored. The TMD monolayers based on Zr, Hf, Pd, and Pt preferentially crystallize in the trigonal T phase (space group P3m1, ITCA number 164; Fig. 1a) [30][31][32][33] and exhibit indirect band gaps. The basic electronic structure and excitation spectra of Zr and Hf based boundary phase TMDs have been studied previously [34]. Some theoretical [35] and experimental [36,37] studies on the electronic properties of the layered bulk alloys exist but to the best of our knowledge the thermodynamic and electronic properties of T-TMDs monolayer alloys have not been systematically addressed.
A special subclass of monolayer alloys are Janus monolayers, i.e. ordered MXX' compounds where the chalcogens X and X' occupy the top and bottom layer sides of the TMD sheet. Janus monolayers have been suggested to be useful for band alignment engineering [14] and as catalysts for the hydrogen evolution reaction [38]. While . The alloys consisted of either (c-d) two species on the metal (M) site and one (no alloying) on the chalcogen (X) site, or (e-f) one species (no alloying) on the M site and two species on the X site. For space group 187, we considered Ti, Zr, Hf, Pd, and Pt on the M site, but we never mixed Ti/Zr/Hf with Pd/Pt. For space group 164, we considered Ti, Zr, Hf, Mo, and W. The X site was always occupied by S, Se, or Te, either just one of them (M-site mixing), or a mixture of two of them (X-site mixing).
some Janus structures have been manufactured [39], little is known about the general thermodynamical stability of these phases.
In this work, we aim to provide a more comprehensive perspective of TMD monolayer alloys. We consider 72 binary alloys of the T and H structure types (Fig. 1) and focus on their thermodynamic properties, specifically the underlying phase diagrams, and electronic properties, specifically band edge positions (for the semiconducting systems) and work functions (for the metallic systems). In the next section we outline the methodology used in this work with details provided in the Supplementary Information (SI). This is followed by a description and discussion of the phase diagrams and mixing characteristics in Sect. III A. Models for the critical temperature and the structural categorization are presented in Sect. III B while the computed band edge positions and bowing parameters are summarized in Sect. III C.

II. METHODOLOGY
Here, we study the thermodynamic and electronic properties of TMD alloys with chemical formula MX 2 ( Fig. 1) in the T structure (space group P3m1, ITCA number 164) and the H structure (space group P6m2; ITCA number 187). For the T structure, we consider mixing on the M-site involving Ti, Zr, and Hf (9 systems) and Pd and Pt (3 systems) as well as X-site mixing (15 systems). For the H structure, we consider mixing on the M-site involving, Ti, Zr, Hf, Mo, and W (30 systems) as well as X-site mixing (15 systems). These alloys were selected based on the known stability or metastability of the boundary phases as documented in the Computational 2D Materials Database [32].
We then carried out a further characterization of the electronic properties of these systems, specifically focusing on the band gap and the position of the conduction and valence band edges using DFT calculations and special quasi-random structures (SQSs) to mimic complete (fully random) mixing [40,[42][43][44][45][46][50][51][52][53][54] (Supplementary Note S4). From the latter analysis, we excluded alloy systems with very high critical temperatures, leaving us with 21 and 27 systems in structure types T and H, respectively.

A. Thermodynamic properties
The mixing behavior at zero temperature can be grouped into three different categories: (A) systems that exhibit two boundary phases separated by a miscibility gap and no long-range ordered mixed phase (Fig. 2a) as well as systems that exhibit a long-range ordered phase at 50% composition with (B) in-plane order (Fig. 2b) or (C) out-of-plane order (Fig. 2c), a situation that is com- monly referred as a Janus system [39]. Out of the 72 alloys considered, we find that 49 alloys fall into category A (Sect. III A 1), 12 into category B (Sect. III A 2), and 11 into category C (Sect. III A 3). These categories are discussed in more detail in the following sections.

Category A: Phase separation into boundary phases
M-site mixing with elements from different groups. Among the 38 M-site alloys that are non-mixing at zero temperature, 18 have critical temperatures above approximately 1200 K (Fig. 3a,b), which suggests that these combinations should be difficult to stabilize if at all. This group comprises all alloys that mix species from different groups of the periodic table on the M-site such as Hf and Mo or W and Zr. According to the Hume-Rothery rules, well established for metallic alloys [55], this behavior can be rationalized by the electronegativity differences. The group 4 elements exhibit electronegativities of 1.3 to 1.5 on the Pauling scale [56] whereas the group 6 transition metals have values from 2.2 to 2.4 (Table S4). This indicate that the bonding in TMDs involving group 4 species is much more ionic than in their counterparts from group 6. For some of these alloying systems, e.g., Zr x Mo 1x 2 and Ti x Mo 1x Te 2 , there is even a sign difference in the Born effective charges (BECs) (Table S3), which is a further indication of the incompatibility of the mixing species. The critical temperature of alloys involving combinations of Zr or Hf with Mo or W is only weakly dependent on the specific transition metal combination. Rather the critical temperature is primarily determined by the chalcogen species with S and Te-containing systems exhibiting the highest and lowest critical temperatures, respectively.
M-site mixing with elements from the same group. The remaining 20 non-mixing M-site alloys have critical temperatures below approximately 900 K (Fig. 3a,b). The 9 systems that do not contain Ti even have critical temperatures below 200 K and should thus be readily miscible. The smaller miscibility in Ti-containing systems can be explained by the much larger difference in lattice parameter between Ti and Zr/Hf-based TMDs (Table S3), which is mirrored by the electronegativity differences (Table S4). In fact, as will be discussed in Sect. III B, the lattice parameters of the boundary phases are in general strong predictors for miscibility. As in the case of M-site mixing with elements from different groups, the critical temperatures decrease monotonically from S via Se to Te. X-site mixing. The critical temperatures for nonmixing X-site alloys are generally very low, with the sole exception of the MS 2x Te 2(1x ) alloys with Pd and Pt that still exhibit values of 550 K and 870 K, respectively (Fig. 3c,d). Generally, the highest critical temperatures are obtained for S-Te alloys, which can be attributed to the large lattice mismatch between the MS 2 and MTe 2 boundary phases. We note that according to the Hume-Rothery rules a 15% difference or more in the covalent radii should indicate immiscibility. Yet, for the present systems the covalent radius is not a good predictor for how well Te mixes with other chalcogen species. This is confirmed by the analysis in Sect. III B.

Category B: Systems with in-plane ordering
The systems in this category include both M-site and X-site mixing that all exhibit an intermediate phase at 50% with in-plane ordering (Fig. 2b).
There are 4 M-site alloys, namely H-Mo x W 1x X 2 with X=S,Se, and Te as well as H-Hf x Zr 1x Te 2 (Fig. 3b), and 8 X-site alloys, including H-MS 2x Se 2(1x ) with M=Mo, W, and Zr as well as five Hf-based systems (Fig. 3c,d).
The only X-site Hf-based alloy that does not exhibit this type of long-range order is T-HfSe 2x Te 2(1x ) .
Generally, the ordering tendency is weak and the critical temperatures are all below 200 K. For almost all practical purposes these systems can be therefore considered as fully miscible.
3. Category C: Systems with out-of-plane (Janus-type) ordering The systems in this category have the general composition MXX' where one chalcogen species (X) occupies the lattice sites above the transition metal and the other chalcogen species (X') occupies the sites below.
All Ti and Zr-based X-site alloys in both T and Hstructure types with the single exception of H-ZrSSe have a Janus ground state at 50% composition. The critical (order-disorder) temperature was obtained from MC simulations by following the difference in composition in each individual layer. The thus obtained critical temperatures are above room temperature for T-TiSTe, T-TiSeTe, H-TiSTe, and H-TiSeTe with value as high as 640 K for TiSTe, indicating remarkable stability. To the best of our knowledge, these Janus monolayers have not been synthesized yet. On the other hand, we find MoSSe, which has in fact been synthesized [39], to be thermodynamically unstable in the present study (Fig. S5). While the Janus monolayer is in fact the isomolar alloy with the highest energy, its energy is only about 13 meV/site equivalent to about 150 K and thus the driving force for decomposing such a structure is very small. This supports the notion that the growth conditions play a crucial role in enabling the formation of these structures.
It is instructive to consider the factors that lead to thermodynamically stable Janus monolayers by comparing the prototypical cases of H-MoSSe and T-ZrSSe. The basic parameters of the boundary systems in these two cases are very similar (Table S3). The most striking difference is observed the BECs, which differ qualitatively. While in H-MoSSe the elements of the BEC tensor are negative as has been shown for the boundary phases also in other studies [57], in T-ZrSSe these elements are positive. This finding is supported in Sect. III B, where we construct a classifier for phase behavior, which reveals that the BECs are important indicators for the appearance of thermodynamically favorable Janus-structures.

Comparison with experiment
Experimentally very little data is available with regard to the thermodynamic properties of TMD alloys, which can be attributed to the considerable difficulties associated with such measurements. It has been observed that phase separation in MoS 2x Te 2(1x ) occurs mainly on the Te-rich side [58], which was suggested to be due to formation of the distorted 1T' phase. In the present study, we have only considered isostructural phase diagrams but nonetheless we observe an asymmetric binodal (Fig. S5). At 80% S the alloy should phase-separate at around 275 K, whereas for 20% that temperature rises to around 325 K. This suggests that the experimentally observed asymmetry could be already be explained by the thermodynamics of the isostructural system.
For the Mo 1−x W x S 2 alloy, random mixing at room temperature for crystals grown with chemical vapour transport was recently shown [27], which is what is expected from the present determination of the critical temperature for the order-disorder transition of 30 K.
The present study also has potential implications for the understanding of the formation and stability of twodimensional lateral heterostructures. It is apparent that the existence of a miscibility gap is not a prerequisite for forming lateral heterostructures with sharp interfaces. Such a configuration has, for example, been realized between MoS 2 and WS 2 [59,60], for which both the present and earlier calculations predict a clear preference for mixing (Fig. S4). The small negative mixing energy of down to −6.5 meV/site in combination with the mixing entropy should imply that such interfaces roughen and vanish at sufficiently temperatures when sufficient kinetic activation enables the system to achieve a thermodynamically more favorable configuration.
On the other hand, for phase separating systems, which can be synthesized in mixed form under nonequilibrium growth conditions, high-temperature annealing, may be used to construct temperature stable lateral heterostructures [61]. In this context, Fig. 3 and Fig. 6 provide guidance as to which heterostructures can and cannot be manufactured using annealing. For example, the proposed lateral heterostructure of T-ZrS 2 and T-HfS 2 [62] is likely to require synthesis techniques beyond thermal annealing due to the low critical temperature of 30 K for the T-Zr x Hf 1−x S 2 alloy.

Comparison with previous modeling studies
The critical temperatures of the MS 2x Te 2(1x ) and MSe 2x Te 2(1x ) alloys with M=Mo or W has been computed before in Ref. 20, which used a mean field approximation of the configurational entropy and the PBE exchange-correlation functional. Here, we go beyond a mean field approximation by using MC simulations that explicitly sample the compositional configuration space and use a more advanced exchange-correlation functional that can be expected to provide a better description of these systems [47]. Overall we find qualitative agreement, and in the case of the Se-Te alloys there is a satisfying quantitative agreement as well. We predict that the critical temperature is 270 K (200 K) to be compared with 360 K (280 K) for the W (Mo) based alloy. For the S-Te alloys, however, we find considerable lower critical temperatures of 370 K (240 K) to be compared with 690 K (490 K) for the W (Mo) based alloy. It is also apparent that the present calculations yield a more asymmetric phase diagram, as already noted above. This suggests that the explicit treatment of the configurational entropy and the exchange-correlation are likely to play a role here.
It has been previously estimated based on a mean field approximation of the configurational entropy that the critical temperature for the order-disorder transition in Mo 1−x W x S 2 alloy for x = 1/3 and x = 2/3 is around 40 K [23], which is very close to the value of 30 K obtained here.

B. Models for predicting alloying behavior
In the previous sections we already indicated how certain properties of the boundary phases can serve as qualitative predictors for the properties of the mixed system. Now we aim to analyze such relations more quantitatively. In the case of metallic (binary) alloy the Hume-Rothery rules are well established as a set of basic conditions under which one can expect mixing. According to these (qualitative) rules mixing ought to be possible if (i ) the relative difference of the covalent radii of the components is less than 15% and the components have (ii ) similar electronegativity as well as (iii ) the same electronic valency. Even though these rules have been formulated for metallic systems, they can still partially explain the mixing behavior in several of the TMD alloys considered here. The Hume-Rothery conditions are, however, based on the properties of the free atoms and as such do not account for any many-body interactions in the boundary phases. For example, in the case of MSTe alloys, the covalent radii of the chalcogens exhibit a 30% difference yet the critical temperature ranges from around 200 K to around 1000 K indicating that the contributions from interactions in the boundary phases are crucial for the alloying behavior.
It can be a daunting task to identify the features of the boundary phases that are the most suitable for predicting alloying ability. Thankfully there are machine learning methods that are very well suited for this task. For this purpose, here, we employ the sure independence screening and sparsifying operator (SISSO) approach that provides an automatized procedure for feature selection to construct predictors for the critical temperature as well as the alloy category (Fig. 3) [63]. To this end, we considered the following properties of the boundary phases: lattice parameters (a 0 ), covalent radii, BECs, electronegativities, unit cell areas (A), sheet modulus (B), ionization potential (IP), electron affinity (χ), and work function. Original features were constructed as the magnitude of the difference between the boundary phases. Since we found the the sign of the BECs to be relevant in the case of the Janus monolayers (Sect. III A 3) for the categorization model we constructed a signed BEC features as where Tr(Z * 1 ) > Tr(Z * 0 ) for the two boundary phases indicated with subscripts 0 and 1, respectively. This quantity provides the information regarding the BECs, removes the ambiguity due to large anisotropy of the monolayers, and at the same time contains information regarding the sign of the BECs.
We note that some of the considered features are highly correlated, e.g., the lattice parameter depends on the character of the bonding in the solid, which is also connected to the BECs, electronegativity, and bulk modulus.
With the exception of only three systems the critical temperatures of the X-site alloys considered here are be- low 370 K. We therefore only consider M-site alloys here, for which we seek a linear relationship between between the descriptor to be determined and the critical temperature. The descriptor that we find via SISSO is which yields a root mean square error (RMSE) of 292 K and a coefficient of determination R 2 of 0.956 (Fig. 4). The descriptor is not unique and depends on the features considered in the model as well as the set of alloys considered as discussed in the original reference [63]. More crucially it, however, provides clear indications as to which boundary phase features to consider when rationalizing alloying behavior. Specifically, the descriptor here emphasizes that ionization potential (or the valence band edge position) as well sheet modulus and unit cell area are important parameters. In addition we constructed descriptors to enable a categorization of X-site alloys according to the classification introduced in Sect. III A. The resulting model has a 93% success rate in partitioning the different categories of mixing behaviour (Fig. 5) and uses two descriptors. The first descriptor (abscissa in Fig. 5) obtained via SISSO is given by It is dimensionless and contains features directly related to bonding (γ BEC and |δχ|). The second descriptor (ordinate axis in Fig. 5) is and has units of pressure. Unfortunately, neither of these two descriptors are amenable to a clear physical interpretation.  Two systems are incorrectly categorized. HfS 2x Se 2(1x ) is categorised as both a Janus and an ordered system, while in reality it is an ordered system yet the Janus structure is only ≤ 5 meV f.u. −1 than other configurations. PdS 2x Se 2(1x ) (categorised as both non-mixing and Janus) on the other hand is a non-mixing system with a critical temperature of 80 K, but the Janus structure is 30 meV f.u. −1 higher in energy than other configurations.

C. Electronic structure
Following the analysis of the thermodynamic properties, we carried out a systematic evaluation of the valence and conduction positions as a function of composition. From this analysis we excluded systems with very high critical temperature, leaving us with 48 systems (Fig. 6).
To first order the band edge positions change linearly with composition between the boundary phases, a behavior that is often referred to as Vegard's law. The deviation from this dependence is commonly described by the bowing parameter b, which is defined via where ε represents the concentration dependent quantity in question. Note that according to this (standard) definition −b/4 corresponds to the deviation from a linear interpolation at x = 0.5.
Band edge bowing parameters for MX 2x X' 2(1−x) for M=Mo,W, and X,X'=S,Se,Te were already reported in Ref. 20. The current results are overall in good agreement with these data, albeit slightly smaller in magnitude (Table S5).
The general trend for X-site alloys is that the valence band bowing is negative and increasing in magnitude in the order S 2x Se 2(1−x) , Se 2x Te 2(1−x) , S 2x Te 2(1−x) (Fig. 6c,d). Especially for the S 2x Te 2(1−x) alloys, the magnitude of the bowing parameter can exceed values of 1 eV. For the conduction band the variations are generally smaller and there is no apparent trend that extends through all TMDs. For example, for T-type Zr and Hf-based TMDs the bowing parameters for S 2x Se 2(1−x) , Se 2x Te 2(1−x) , and S 2x Te 2(1−x) are all in range between 0.07 and 0.09 eV. Furthermore, we find that for X-site alloys, the bowing parameter of the valence band is always negative, and that for H-type X-site alloys, the conduction band bowing parameter is always negative.
For T-type X-site alloys, the conduction band bowing parameter can have either sign. For example, in T-HfS 2x Se 2(1−x) , the valence band exhibits a bowing parameter of −0.34 eV while the conduction band bowing parameter is 0.09 eV.
For T-type M-site alloys, the magnitude of the bowing of the valence band edge is generally small, with the maximal (absolute) value of −0.05 eV occurring in T-Pd x Pt 1−x S 2 (Fig. 6a,b). The conduction band edge exhibits significantly larger bowing in this class of alloys with values up to 1.10 eV (in T-Pd x Pt 1−x S 2 ). Finally, for H-type M-site alloys, the valence band bowing is larger in magnitude than for T-type M-site alloys with a maximal magnitude of −0.41 eV found in H-Ti x Zr 1−x S 2 . The largest bowing of the conduction band for the H-type Msite alloys is found in Ti x Hf 1−x S 2 with a value of 0.29 eV.
The valence band and conduction band variations that are possible in the considered alloys are illustrated in Fig. 6. Specifically, we find that the variation of the valence band can be very large in MS 2x Te 2(1x ) , e.g., in HfS 2x Te 2(1−x) the difference between the boundary phases is 1.7 eV. In general the valence band position increases in the series S-Se-Te. The typical behavior is illustrated for T-ZrS 2x Se 2(1x ) (Fig. 6e) and H-TiSe 2x Te 2(1x ) (Fig. 6f).
It has been noted before that the magnitude of the valence band bowing parameter increases in the order S 2x Se 2(1−x) , Se 2x Te 2(1−x) , and S 2x Te 2(1−x) for H-type Mo and W-based X-site alloys [20]. Our results show that this trend also holds true for Ti, Hf, and Zr-based alloys that exhibit H symmetry and Zr, Hf, Ti, Pd, and Pt alloys that exhibit T symmetry. The band edges exhibit in most cases moderate bowing with the exception of the X-site alloys MS 2x Te 2(1x ) . These are the X-site mixed systems that exhibit the largest lattice constant mismatch and largest differences in boundary phase features. Therefore, these alloys exhibit the largest variability of the considered alloys. The large bowing parameter in these compounds may however complicate band edge tuning since a small variation in composition is associated with a relatively large variation in band edge position.

IV. CONCLUSIONS
In this study, we have provided a comprehensive study of the phase diagrams and electronic structure of mono-layer TMD alloys. It has been shown that mixing systems with in-plane order are absent at room temperature but that specifically the Ti-based X-site alloys may exhibit a Janus ground state that remains the most stable structure far above room temperature. Furthermore, we have shown that M-site alloys with transition metals from different groups in the periodic table are associated with very high critical temperatures (> 1000 K) whence it is unlikely or at least very difficult for them to be manufactured. The band edges exhibit in most cases little bowing with the exception of X-site alloys that include both S and Te, implying that in many systems the band edge positions are relatively well approximated by a linear interpolation of the values of the boundary phases.