Pettifor Maps of Complex Ternary Two-dimensional Transition Metal Sulphides

Alloying is an established strategy to tune the properties of bulk compounds for desired applications. With the advent of nanotechnology, the same strategy can be applied to 2D materials for technological applications, like single-layer transistors and solid lubricants. Here we present a systematic analysis of the phase behaviour of substitutional 2D alloys in the Transition Metal Disulphides (TMD) family. The phase behaviour is quantified in terms of a metastability metric and benchmarked against many-body expansion of the energy landscape. We show how the metastability metric can be directly used as starting point for setting up rational search strategies in phase space, thus allowing for targeted further computational prediction and analysis of properties. The results presented here also constitute a useful guideline for synthesis of TMDs binary alloys via a range of synthesis techniques.

Since the discovery of graphene, 2D materials have been at the forefront of Materials Science and Discovery. In addition to fundamental research interest [1], recently their unique properties and reduced dimensionality have sparked an interest for nanoscale engineering applications. Ideas for 2D-materials-based devices can be found in tribology [2], electronics [3] and catalysis [4]. In this relatively new field, there have been so far only limited attempts to exploit the vast chemical space spanned by alloys to optimise properties. Up to now, most research efforts have focused on identifying 2D unaries and binaries both theoretically [5,6] and experimentally [7,8]. However, little is known about their thermodynamic phase behaviour. The structures and ordering of possible alloys are largely unexplored territory [9]. Only few 2D ternaries have been reported by experiments [10,11] and, while a handful of binary alloys has been studied [12][13][14], no systematical analysis has been carried out. But knowledge of thermodynamic behaviour is fundamental for advancing the engineering applications of 2D materials. When properties such as bandgap and electronic transport need to be tuned to desired values by chemical doping, the presence of miscibility gaps and competing ternaries has to be taken into account [15].
The vast crystallographic and chemical spaces need not be explored by experiments alone. Computational tools can provide guidelines to experimental synthesis, reducing the number of possible candidates by orders of magnitude. As an example, Mounet et al. [5] reduced a dataset of 1 × 10 5 bulk geometries from experimental databases to 258 easy-exfoliable monolayer (ML) candidates. As a comparison, large-scale experimental studies usually deal with dozens of candidates [7,8].
In the last century, the discovery of new metallic alloys was guided by empirical methods like the Hume- * a.silva@soton.ac.uk † d.kramer@hsu-hh.de Rothery rules [16] and Pettiford maps [17]. These rules are based on atomic proprieties like relative ionic size and electronegativity, combined through chemical intuition and experience. The somewhat surprisingly wide validity of these simple rules in metallic alloys has been proven by experiments in the 1940s. With the advent of Density Function Theory (DFT) and Cluster Expansion (CE) methods in the 1980s, the physics underpinning the phase diagram of metallic alloys was explored systematically, with a symbiotic relationship between experiments and simulations [18]. Nowadays, we are able to create large databases of materials and rationalise complex trends coupling the predictive power of DFT, the massive improvement in computation power and the availability of software tools. These capabilities, along with experimental validation, should allow us to build on the Hume-Rothery and Pettiford rules and extend their concepts to novel classes of materials. Indeed there are examples of such efforts in recent literature: the known empirical rules have been cast in terms of well defined probabilistic models trained on large computational datasets [19] or extended to include the physics of oxides [20].
Here, a framework is presented and a dataset compiled to explore alloy possibilities for the TMD family, the most widely studied 2D material family for engineering applications. The article is structured as follows. The first section defines the chemical and coordination spaces considered. Then, a metric to quantify metastability and solubility tendency in different hosts is developed. The metric is applied to the chemical and coordination space defined in the first section, yielding the host most receptive for alloying for each transition metal (TM) pair. In the third section, the CE formalism [18] is used to benchmark the predictions of our metric and to identify stable orderings. For illustration, an analysis of the phase behaviour is presented here for four representative alloys. The Supporting Information (SI) contains further examples. Finally, the article concludes with a discussion of how the framework could guide synthesis efforts. arXiv:2207.02741v1 [cond-mat.mtrl-sci] 6 Jul 2022

I. CHEMICAL AND COORDINATION SPACES
The starting point to build the space of possible compounds is the 2D-materials database compiled by Mounet and coworkers [5]. The database comprises 258 mechanically stable ML structures identified from experimental bulk compounds. Thus, the following phase stability study is conducted on ML geometries only.
To reduce the computational effort, the selection of the possible prototypes and elements to mix is guided by knowledge in the literature [5,8,21,22] and the original database is filtered according to the class of materials of interest. Here, the database is scanned for compounds of the form M n A 2 , where M is a TM cation (highlighted in Figure 1a) and A is the anion, oxidising the TM (see Section I in the SI for the list of anions considered). In selecting the prototypes, the possible cations are restricted to the transition metals considered but the anions are not limited to sulphur, as layered prototypes that could host TMD alloys may not be expressed in terms of sulphides in the database (see Section I and table SII in the SI for details). This search yields the eight prototypes shown in Figure 1b-i, whose space group is reported in Table SIII of the SI. While here the symmetry of each prototype is frozen, focusing on the substitutional degree of freedom, it is in principle possible to identify pathways between these crystal structures allowing for phase transitions between the prototypes [23].
Intermediate TMs (Cr, Mn, Fe, Ru, Os) are considered here although they do not form layered sulphides on their own but might form ML alloys in combination with other TMs, e.g. Fe-doped MoS 2 ML [21]. Late transition metals from group XI onward are excluded, as they do not bind with chalcogenides to form layered materials [8]. This yields the N = 21 TMs highlighted in Figure 1a as a possible cations M in the M S 2 stochiometry.
While the methodology described here is valid for any stochiometry and cation-anion selection, our analysis will focus on M S 2 compounds, as these are the most frequently synthesised and studied compounds of the family. This selection yields TM × prototypes × chalcogenides = 168 binaries as a starting point for TM 1 × TM 2 × prototypes = 3528 substitutional alloys on the TM site. The total number of candidates, although large from an experimental point of view, allows for an exhaustive theoretical analysis rather than approximate methods based on a statistical sampling of configurational space [24].

A. Lattice stability
The total energy of each compound M S 2 in all prototypes p, i.e. pairs (M, p), is obtained from Equation of State (EoS) calculations. The volume range considered in the EoS is determined using the notion of covalent radius r c of the element i. The protocol is described in Section II of the SI.
The energy above the ground state of each compound M S 2 in a given prototype p, also known as lattice stability [25], is given by the total energy per site with respect to the ground state (GS), i.e.
where E(M, p) is the minimum energy of M S 2 compounds in prototype p obtained from EoS calculations and n is the number of sites in the metal sub-lattice, i.e. the number of TM in the unit cell. The offset energy for each TMD E GS (M ) is the minimum energy across the prototype space E GS (M ) = 1 n min p E(M, p) for layered TMD and the total energy of the 3D bulk structure E GS (M ) = 1 n E 3D (M ) for non-layered TMDs. The non-layered TMDs are identified by comparing the minimum-energy 2D prototype across the considered ML geometries to the GS reported in the Materials Project (MP) database [26,27] for the given M S 2 compound. The 3D geometry has lower energy than the relative 2D GS for six metal disulphide, namely FeS 2 , CoS 2 , RuS 2 , RhS 2 , OsS 2 , IrS 2 (TM in gray boxes in Figure 1). An analysis equivalent to the one presented here but restricted to 2D geometries is reported in Section IX of the SI as it might be relevant for experimental techniques able to bias the synthesis towards atomically thin films [31].
For the layered TMDs, the binding energy between the layers (typically around 10 meV/atom for TMDs [28,29]) is neglected here, since this offset does not affect the ML phase behaviour [30]. Figure 2 reports the energy above the ground state per lattice site defined in Eq. (1) for the selection of TMs and prototypes shown in Figure 1. Each column shows the energy above the ground state of the given TM in the eight prototypes with respect to the identified 2D GS. Green squares mark the GS of layered TMDs and orange squares mark the lowest-energy 2D prototype of non-layered TMDs. As a guide to the eye, each entry is colored according to its energy above the ground state, as reported by the colorbar on the right, and periodic table rows are separated by vertical dashed lines.
The ground states of known layered compounds are identified correctly according to the MP database: d 2metal TMDs (TiS 2 , ZrS 2 and HfS 2 ) display octahedral CdI 2 coordination, Figure 1c. The MoS 2 prismatic prototype, Figure 1d, is the GS of d 4 TMDs, while the d 10 metals Ni and Pd are found to favour the square planar PdS 2 prototype, Figure 1b [26,27]. A systematic comparison of the predicted ground state for the 21 pristine compounds in Figure 2 with experimental and computational data available in the literature [26,27,[51][52][53][54][55][56][57][58][59][60] indicates that our protocol correctly describes the energetics of the considered chemical space. The comparison is reported in Section III and Table SIV of the SI. Moreover, the larger steric hindrance of heavier TMs in the same group raises the energy above the ground state of unstable prototypes. This can be observed by following the row relative to prototype PdS 2 in    Figure 1a. The colorbar on the right reports the energy above the ground state in eV over lattice sites. Green squares mark GS prototypes, defined by EF = 0. Orange squares mark the lowest-energy 2D prototype of transition metals displaying a 3D GS (gray boxes in Figure 1a). Vertical dashed black lines separate rows of the periodic table, see  Finally, it is important to appreciate the scope of validity and the possible sources of errors in the dataset presented here. The DFT calculations performed are spin-polarised, thus non-magnetic and ferromagnetic groundstate are correctly described. Antiferromagnetic (AFM) orderings are not considered, as calculations are performed in cells comprising a single TM site. The only AFM orderings for the considered stoichiometry are reported for NiS 2 and MnS 2 [61]. While important for materials properties, AFM GS in layered TMDs are usually almost degenerate in energy with FM states [61]. Moreover, no Hubbard correction (GGA+U) is included here. The effect of Hubbard U on the relative total energy for the considered TMD stoichiometry is negligible [61], but a detailed benchmark must be carried out when applying our protocol to different stoichiometries, as discussed in the Methods section.

II. IDEAL SOLID SOLUTION LIMIT
Starting from the lattice stability matrix in Figure 2, a question arises naturally: is it possible to identify which metals are likely to mix in a given prototype? A straightforward approach to explore this question is the ideal solid solution limit, a non-interacting model based on the relative energy of pristine TMDs defined in Eq. (1). Given a binary alloy in a prototype p, M x Q 1−x S 2 | p , the ideal solid solution represents a model with negligible interactions between the fraction x of sites occupied by M and the remaining 1 − x sites occupied by Q. In the energy-composition space, the system behaviour is represented by the line connecting the energy above the ground state of QS 2 in prototype p at x = 0 with the energy above the ground state for M S 2 at x = 1 in the same prototype, i.e. the element (Q, p) and (M, p) of the matrix in Figure 2, respectively. Hence, in the ideal solid solution model, the energy above the ground state of a mixed configuration at concentration x is given by: By construction, this energy is exactly zero everywhere if M and Q share the same ground-state structure p, In any other case, the energy will be positive: suppose the metal M has a ground-state geometry p = p, the fraction x of material M S 2 | p would transform into p to reach equilibrium at zero temperature. The model effectively quantifies the metastability at zero temperature of alloys in a selected prototype p as a function of concentration x. By construction, this model cannot predict stable mixtures, i.e. negative formation energies, but can be used to estimate the likelihood of solubility and phase separation in a system: the lower the metastability of the solid solution model, the smaller any entropic or chemical stabilising mechanisms must be to stabilise alloys under synthesis conditions. As an example, let us consider the effect of finite temperature in the solid solution model. The equilibrium of an alloy in the prototype p at temperature T is determined by the free energy where the substitutional entropy of a binary alloy is a function of the concentration x only, independent of the elemental pairs: which counts possible configurations of the two atom types on the metal sub-lattice [32]. If there exists a concentration and temperature (x * , T * ) at which T * S(x * ) > E 0 Q,M,p (x * ), then the free energy becomes negative and the mixture is thermodynamically stable (see Section IV of the SI for an example). Note the free energy F Q,M,p (x, T ) of different hosts p intersect at the same composition x found for the energy above the ground state E 0 Q,M,p (x), the entropy S(x) being a function of concentration only. Thus, the simpler linear energy model in Eq. (2) yields the same relative energy ordering of the prototypes, as shown in SI Section IV. For an example of an electronic-driven stabilisation mechanism present also at zero temperature, see the discussion in the SI Section X.A.3 of the ternary GS of (Mo:Nb)S 2 and (Mo:Ta)S 2 shown in Figure 3a,b.

A. Metastability Metric
In order to make the relative metastability between prototypes quantitative, a metric in the compositionenergy space is needed to compare different combinations. Consider a prototype p and two metal sulphides M S 2 and QS 2 with GS prototype p M and p Q , respectively. The convex hull across all phases in the concentration-energy space is the line E = 0 connecting the energies of the end-members in their respective GS prototypes, dashed gray lines in Figure 3. A point on this line at the fractional concentration x = 0, 1 represents a phase separating system where the fraction x of M S 2 is in its GS prototype p M and the remaining 1 − x is in its own GS p Q . For a configuration to be stable, its energy must be lower than this hull. As our model by definition cannot break this hull, we characterise the metastability of a model alloy by its positive energy above the ground state, i.e. its distance from the hull [33].
We define a descriptor intended to capture the energetic "disadvantage" of a particular prototype (p, Q, M ) relative to the relevant binary ground states as follows. The metastability window of the (p, Q, M ) triplet is defined as the range of concentration x where the distance from the hull in Eq. (2) within the prototype p is lower or equal to the distance from the hull within the groundstate prototypes p M and p Q , as shown by blue regions in Figure 3.
Let us apply this construction to an example: consider the energy above the ground state in the solid solution model of the (Pd:Nb)S 2 alloy in Figure 3a. The blue line refers to the energy above the ground state of the CdI 2 prototype, while the red and green dashed lines refer to the ground state of the PdS 2 end-member x = 0 (PdS 2 prototype) and NbS 2 end-member x = 1 (MoS 2  prototype). The metastability of the prototypes varies as a function of the concentration. Near the respective end-members, the ground-state prototypes are favoured, e.g. the PdS 2 prototype has lower distance from the hull in the range x ∈ [0, 0.2]. The CdI 2 prototype lies closer to the hull in range x ∈ [0.2, 0.9], suggesting that a metastable solution in this range in this prototypes is more likely than in either of the two ground-state prototypes. When the two TMDs share the same prototype GS, the distance from the hull in that prototype is zero everywhere, like in Figure 3b. In this case the metastability window extents from 0 to 1, suggesting that solubility is likely. When the prototype p is the ground-state for one of the metals, the metastability window extends from the extremal concentration, x = 0 or x = 1, up to the intercept with the energy above the ground state in the other prototype, as shown in Figure 3c. Finally, a metastability window might not exist for a given triplet, as shown in Figure 3d: the distance from the hull in the FeO 2 prototype is higher than in either ground-state prototypes for any concentration. In this case, the formation of alloys within this prototype is unlikely.
Applying the construction depicted in Figure 3 to all TM pairs yields a N × N matrix, for each prototype p. Each entry of these metastability matrices are a 2 × 2 matrix containing the bounds of the metastability window and the energy above the ground state in Eq. (2) evaluated at the metastability limits, i.e. minimum and maximum hull-distance within the window. The matrices associated with each prototype are reported in Section V and the dataset of the SI.

B. Optimal Prototypes
Given a pair of TMs, the prototype most receptive for alloying can be identified by comparing the metastability windows in different prototypes build from the metastability metric in the previous section. A function associating a score to each metastability window needs to be defined in order to rank different prototypes. This ranking has to assign a single value to the metastability windows of TM 1 -TM 2 -prototype triplets. The following parametric function is chosen as goal function where w is the width of the metastability window and the energy penalty is the hull-distance of the centroid defined by the window in the energy-concentration space, i.e. blue points in Figure 3. Thus, the function encourages large metastability windows w and discourage large energy penalties . Details regarding the goal function and the selection of the appropriate weight ζ for the present dataset are reported in Section VI of the SI. The optimal prototypes for each pair of transition metals, selected by f ζ with ζ = 0.080 eV/site, are shown in Figure 4. The symbol assigned to each entry refers to the optimal prototype, as shown in the lower legend; symbols on the diagonal mark the 2D-GS prototype for that transition metal. The size of each marker shows the width of the metastability window associated with that metal pair in that prototype. The colour code of each Q, M entry shows the energy above the ground state of Q 1−x M x S 2 at each end of the metastability window, as indexed by the metal on the horizontal axis. For example, consider the Ti 1−x Ta x S 2 binary in the CdI 2 prototype, whose energy landscape is reported in Figure 3c. Follow the green lines in Figure 4 to the entry in the upper triangle, Ta row and Ti column. This entry shows the energy above the ground state on the Ti-side of the metastability window, left-hand-side in Figure 3c. Since the CdI 2 prototype is the GS of TiS 2 , the energy above the ground state on this side is zero, indicated by deep blue color. Conversely, the entry in the lower triangle, Ti row and Ta column, shows the energy above the ground state on the Ta-side of the metastability window, right-hand-side in Figure 3c. Since the CdI 2 prototype is not the TaS 2 native prototype, the energy on this side is positive, light-blue color.
The same procedure, following the yellow lines, applies for Mo 1−x W x S 2 binary in MoS 2 prototype, whose energy landscape is reported in Figure 3b. The end-members share the same GS, hence the plot shows two large, deep blue symbols with zero energy penalty. Figure 4 provides a visual tool to navigate the possible mixtures of transition metals within the sulphur planes. Large blue marks in Figure 4 indicate a small energy penalty in the metastable window, and, thus, that miscibility between the two metals within the S host is likely. For example, in the case of TiS 2 (GS prototype octahedral CdI 2 ) and TaS 2 (GS prototype prismatic MoS 2 ), Figure 4 indicates good miscibility in the CdI 2 prototype, that can be traced back to the relatively low energy above the ground state of TaS 2 in the TiS 2 native prototype, E F (Ta, CdI 2 ) = 0.06 eV/site, see the lattice stability in Figure 2. On the other hand, a high energy penalty and small metastable window likely results in miscibility gaps. These likely phase-separating systems constitute the missing elements in Figure 4.
The distinction between likely-mixing and likelyseparating systems can be made more quantitative by extending the Hume-Rothery rules to our case. Following the original rules, miscibility between transition metals within the sulphur host is expected if the lattice mismatch between the pristine compounds is less than 15 % [16] (see SI Section VII for definition and values of the mismatch in these compounds). Moreover, we extend the original rules using the metastability metric of the prototype. Following the work by Sun et al. [33] on metastability of inorganic crystals, we set a threshold of E = 120 meV/site as an upper limit for the energy above the ground state of the optimal prototypes, as metastable compound within this range have been observed experimentally. As a result, Figure 4 features missing elements where the optimal prototypes are unlikely to be receptive to alloying due to large lattice mismatch or high energy above the ground state. Since experimental formation energies on these compounds are scarce, the threshold proposed here are tentative values that can easily be updated with novel experimental data. The unfiltered matrix is reported in Section VIII of the SI.
As a first benchmark, the information in Figure 4 can be compared with alloys reported in the literature. We focus on alloys of MoS 2 , as many alloys for this wellknown system are reported; consider the relevant column in Figure 4, highlighted by the leftmost yellow line. Zhou and coworkers [7] recently reported synthesis of ML of (Nb:Mo)S 2 , which is shown as likely to mix in Figure 4. However, the same work reports a (Mo:Re)S 2 ML alloy, while the metastability window of this TM pair is small and high in energy (≈ 350 meV/site in Figure 4 (and Figure S13b in the SI). Another recent work [34] reports the experimental characterisation of (V:Mo)S 2 ML, which is also a TM pair likely to mix according to our analysis.
Onofrio and coworkers [22] compiled a dataset of possible substitutional alloys of 1H-MoS 2 ML throughout most of the periodic table using DFT methods. According to the authors' analysis, based on substitution in the smallest possible unit cell (roughly x = 0.5), compounds based on all early TMs between group III and group VI show negative formation energy. The authors prediction for metals of group V (V, Nb, Ta) and group VI (Cr, W) agree with our metastability metric. In contrast, group IV elements (Ti, Zr, Hf) show a low likelihood of miscibility according to Figure 4, while Ref. [22] report negative formation energies. The case of (Mo:Ti)S 2 is discussed in more detail below, showing that the prediction of our metric agrees with CE analysis and available experimental data.

C. Polymorphism
The information in Figures 2 and 4 can be coarsegrained to understand the tendency of different TMs to stabilise foreign hosts in mixtures. Given a metal M , the energy cost of forming meta-stable phases as pure M S 2 is given by the columns of Figure 2, that report energy above the ground state of each M S 2 compound in the considered hosts p. For example, consider the first column in Figure 2. TiS 2 , whose GS is the perfectly octahedral CdI 2 , exhibits a low energy penalty for the distorted octahedral coordination of WTe 2 , E F = 0.02 eV/site. For MoS 2 , whose GS is the prismatic coordination, the lowest-energy meta-stable prototype is distorted WTe 2 (E F = 0.55 eV/site) and perfect CdI 2 octahedral displays a higher energy above the ground state of E F = 0.84 eV/site. The WTe 2 polymorph has indeed been observed experimentally [4] and the CdI 2 one has been reported in simulations of MoS 2 layers at high temperature [35].
Similarly, the metastability metric helps to evaluate the tendency of a metal M to stabilise non-native hosts when alloyed with a second metal Q. Purple-shaded marks in Figure 5 report the minimum centroid energy penalty across all possible combinations TM 1 -TM 2 -p, for each TM 1 -p pair. A low centroid energy of a given prototype p (x axis) suggests that the considered metal M (y axis) could potentially stabilize this prototype when mixed with another metal in the sulphur host. Figure 5 confirms the meta-stable tendencies highlighted in the previous paragraph. The lowest-lying prototype for both Ti and Mo is WTe 2 , meaning that alloys in this prototype could be stabilised by the presence of these metals. A relatively low energy penalty for the CdI 2 prototype is observed in group V TMDs (VS 2 , NbS 2 , and TaS 2 ). This suggests that these TMDs could be receptive for alloys in these meta-stable coordinations, alongside the native MoS 2 prototype.

III. METAL SITE ORDERINGS
The phase behaviour predicted by the metastability metric reported in Figure 4 can be benchmarked by exploring the stability of possible orderings and miscibility regions using a many-body expansion based on electronic-structure calculations. The formation energy of a pseudo-binary system M x Q 1−x S 2 is modelled with the CE formalism [18]. The interaction between different species on the TM site sub-lattice, like the triangular one formed by orange and blue circles in Figure 6, is modelled via a set of many-body interactions, termed clusters, e.g. the pairs α and β and the triplet γ in Figure 6. The sulphur atoms, yellow circles in Figure 6, are spectators, i.e. they are considered in the DFT total energy calculations but not in the CE interaction figures. The GS end-members are taken as reference to compute the formation energy of the ordered configuration σ(x) at concentration x in M x Q 1−x S 2 : (5) where E(σ(x))| p is the total energy of the configuration σ(x) in the host lattice defined by the prototype p. E(M, p M ) and E(Q, p Q ) are the total energies of M S 2 and QS 2 in their GS prototypes, p M and p Q , respectively. This chemical reference assures that the formation energy in Eq. (5) at end-member concentration x = 0 and x = 1 corresponds to the energy above the ground state reported in Figure 2.
The set of geometrically distinct orderings is generated using CASM [36][37][38]. The geometries are fully relaxed, including cell shape and volume. The dataset is updated iteratively with stable orderings suggested by the CE model until predicted and computed convex hulls coincide. For details see the Methods section.
The following section reports our benchmark results, which cover the cases of highly-miscible TMs within the same GS host, a phase-separating system and a system with finite-miscibility of a TM in a non-native prototype.
Two other examples, one of perfect miscibility and one showing the limitation of the CE model, are presented in Sections X.C and X.D of the SI.
A. High miscibility: (Mo:group V)S2 Pseudo-binary Alloys The metastability metric in Figure 4 predicts high miscibility for mixtures of Mo-Group V elements. This class of alloys attracted interest as a possible realisation of MoS 2 -based devices. In particular, (Nb:Mo)S 2 alloys have been indicated as a viable p-doping solution for MoS 2 ML transistors [3,12]. Ta-doped MoS 2 composite coatings have been identified as a promising fatigueresistant material for tribological applications [39].
The computed alloy of both (Mo:Nb)S 2 and (Mo:Ta)S 2 , reported in Figure 7a,b respectively, show novel ternary GS that break the convex hull and low zero-temperature formation energy across the whole concentration range. In particular, on the Mo-rich side (lefthand-side in Figure 7a,b) substantial doping should be achievable at finite temperature, due to the absence of competing ternary ordered configurations. On the Nband Ta-rich side (right-hand-side in Figure 7a,b) the phase diagram is dominated by the ternary compounds breaking the convex-hull (solid lines). These ternaries, Mo 1/3 Nb 2/3 S 2 , Mo 1/3 Ta 2/3 S 2 , and Mo 1/9 Ta 8/9 S 2 , are reported here for the first time to the best of the authors knowledge. However, the small energy scale formally stabilising these ordering at zero temperature make it likely that long range order might be destroyed at room temperature and above (see Section X.A.4 in the SI). A good understanding of the phase behaviour of these systems is needed, especially as the doping concentration needed in p-doped devices may reach 20% [12] and the competition with ternary phases might make synthesis problematic.
One would expect similar behaviour from Nb and Ta dopants, as the two have the same covalent radii, electronic configuration [40] and same lattice parameter in TMD compounds. Indeed the qualitative behaviour is the same for both systems, as predicted by the metastability metric. Quantitative behaviour differs slightly: a single ternary Nb 2/3 Mo 1/3 S 2 breaks the hull in the Figure 7a while the Ta system displays a richer landscape with competing ternaries Ta 2/3 Mo 1/3 S 2 and Ta 8/9 Mo 1/9 S 2 . This quantitative difference arises from subtle electronic differences in the Nb and Ta ions. Modelling these alloys present a double challenge, as one needs to capture at the same time the many-body, nonlocal character of phase stability and long-range elastic interactions due to lattice mismatch between NbS 2 or TaS 2 and MoS 2 . The CE formalism is suited to handle the first task, while the description of elasticity is problematic [41]. Since the CE expansion is performed on a complete representation of the energy landscape of the lattice model, the CE can describe small elastic displacements, at the cost of increased complexity. Indeed, more than a hundred orbits, up to five-vertex clusters, must be included in the model to appropriately describe the convex hull in Figure 7a,b, far more than for the nearcommensurate (Mo:W)S 2 case, as reported in Section X an Table SV  different contributions and of the ternary ground states are reported in Section X.A.3 in the SI.
While the system shows miscibility gaps between stochiometric GS at zero temperature, the small formation energies in the computed configurations, typically E(σ(x)) < k B T room = 0.025 eV, suggests that these miscibility gaps close below usual synthesis temperature T synth ≈ 600 K (see Section X.A.4 in the SI). The metastability metric in Figure 4 can help identify metal pairs that would phase separate rather than form alloys in TMDs. As an example of this behaviour, Figure 7c reports the formation energy of the (Mo:Ti)S 2 alloys. This system has been analysed in detail in our previous computational work in Ref. [30] and characterised experimentally [42].
A high lattice stability energy of Mo in the TiS 2 ground-state prototype and vice versa results in a low score in the metastability metric; see the corresponding missing entry in Figure 4 (or the small light-blue triangle in the unfiltered matrix in Figure S13b in SI). This prediction is confirmed by the CE model in Figure 7c. No configurations in the MoS 2 prototype (blue symbols) display lower formation energy than the solid solution limit (solid blue line). Within the CdI 2 prototype, some configurations display a lower energy compared to the solid solution limit, red crosses between the solid red line and dashed gray line, respectively. This energy gain, however, is not enough to break the inter-prototype convex hull (dash-dotted gray line at E = 0), resulting in an overall phase separating system. The origin of this phase behaviour lies in the different electronic structure in the local environment of the TM, as explained in terms of crystal field levels in Ref. [30]. The CE model trained on DFT data have been used to estimate solubility limits in the phase space as a function of temperature, predicting low miscibility at high temperature, in line with experimental observation [30,42].
C. Cross-host miscibility: (Ti:Ta)S2 Pseudo-binary Alloys Finally, we report an example of cross-host miscibility, i.e. an alloy system between two TMDs that do not share the same GS prototype. This case is identified by combining all the information presented here. The starting point is the polymorphism plot in Figure 5. Group V elements (V, Nb and Ta) show low formation energy in the CdI 2 prototype, which is the ground state of many TMDs (see Figure 2), e.g. group IV elements (Ti, Zr and Hf). Consulting the metastability metric in Figure 4, possible alloying combinations of VS 2 and TaS 2 with any group VI elements stand out as promising candidates, while NbS 2 displays a slightly larger formation energy and can be set aside. Taking also the mismatch into account as stated by the adapted Hume-Rothery rules, the (Ti:Ta)S 2 system is the most promising candidate: the mismatch for (Ti:V)S 2 l VS2 /l TiS2 = 0.928 is larger than for (Ti:Ta)S 2 l TaS2 /l TiS2 = 0.990 (see Section VII of the SI). The full metastability metric construction leading to the high score of (Ti:Ta)S 2 in Figure 4 is also reported in Figure 3c for reference. From Figure 3c, it is also clear that the MoS 2 prototype is unfavorable for TiS 2 probably resulting in phase separation between the two metals in this prototype. This tendency is also visible in the MoS 2 prototype metastability matrix in Figure S3b in the SI.
We now benchmark the prediction from the metastability metric and the updated Hume-Rothery rules against actual alloy configurations from DFT. Figure 7d reports the formation energy of (Ti:Ta)S 2 alloys in the CdI 2 (red symbols) and MoS 2 prototypes (blue symbols). As predicted by the metastability metric, TiS 2 and TaS 2 segregate in the MoS 2 prototype: no configuration lies below the solid solution limit (straight blue line). In the CdI 2 prototype, native host for TiS 2 but not for TaS 2 , the alloyed configurations lie below both the solid-solution line (dotted gray line) and the cross-host solid-solution hull (dash-dotted gray horizontal line) from x ≈ 0 up to x ≈ 0.7. While at zero temperature only the GS on the convex hull (red solid line) are stable, the energy scale is small compared to room temperature, suggesting that solid-solution alloys in the CdI 2 prototype should be possible to synthesise in experiments, e.g. with CVD methods. Indeed, there are reports of (Ti:Ta)S 2 solid solution alloys in the literature [43], although no crystallography data or solubility limits are available to date. This experimental confirmation validates the exploration approach outlined in this section.

IV. CONCLUSIONS
We presented a systematic analysis of possible alloys in the TMD chemical space. The metastability metric provides a simple yet useful picture to guide in-depth computational studies and experimental synthesis. Predictions by the metastability metric are in good agreement with alloy systems reported in literature. Moreover, many-body expansion based on electronic-structure methods of selected binary alloys confirm the predictive power of the metric both in identifying phase separating and highly miscible systems. While this work focused on TMDs, the methodology developed here can be transferred to any stochiometry and composition.
The optimal prototype matrix and the other tools can help to identify viable alloy candidates minimising the trial-and-error attempts, speeding up the progress of nanotechnologies. Section III C demonstrates a possible protocol that could be followed to aid CVD synthesis of novel ML alloys to stabilise TMs in non-native local environments.
In a wider context, the framework developed here fits in the effort of making chemical intuition quantitative. The exploration of a large dataset, easily produced with modern DFT methods, allows to rationalize trends across the periodic table and refine the known empirical rules. In particular, attempting to transfer the Hume-Rothery rules for metallic binaries to the class of 2D TMDs seem attractive. Here, we propose to replace the ionic size with the lattice parameter of the M S 2 crystal. The rules on electron counts and electronegativity are implicitly embedded in the lattice stability differences, along with other more complex descriptors, like the d-band overlap and crystal field effects, as shown in the Section X.A.3 of the SI and in Ref. [30]. This last rule generalisation is based on the predictive power of DFT, that has been the cornerstone of Computational Material Science in the past decades. Here we propose that miscibility is likely if the formation energy of the metastability window defined here is lower than 120 meV/site.
To summarise, we presented a set of tools and ideas that will guide computational chemists and experimentalists in charting the under-explored chemical space of TMDs.

METHODS
a. First principles calculations The total energy calculations are carried out with the Vienna Ab Initio Software Package (VASP) [44][45][46], within the PAW framework for pseudo-potentials [47]. The generalisedgradient-approximation to DFT as parametrised by Perdew, Burke, Ernzerhof [48] was used in this work. The Kohn-Sham orbitals are expanded in a plane-wave basis with a cutoff of E cutoff = 650 eV and the BZ is sampled with a 17 × 17 × 1 mesh. The electronic density was computed self-consistently until the variation was below the threshold of 1 × 10 −6 eV. We perform spin-polarised calculation; the electronic structure can converge to nonmagnetic or ferromagnetic states, as we consider only primitive unit-cells in our calculations. The position of the ions in the unit cell were relaxed until the residual forces were below the threshold 1 × 10 −2 eV/Å. To ensure no spurious interactions between the periodic images, a vacuum of 20Å was added along the c axis.
Note that while error cancellation in the stoichiometric analysis carried out here makes the Hubbard U correction not necessary, Ref. [61] shows that this becomes fundamental in modelling thermochemical reactions involving valance changes, as the reaction enthalpy of most sulphurisation reactions is not correctly described at U=0. b. CE model training The fitting procedure is carried out within the CASM API [36][37][38]. Each configuration σ i is weighted according to its distance from the convex hull: where E(σ i ) is the formation energy of the configuration σ i , E hull (x i ) is the formation energy of the convex hull at the concentration x of the configuration σ i and k BT is a fictitious temperature set according to the energy scale of the problem. These weights bias the fitting towards reproducing more accurately low-energy configurations, which are the relevant ones to capture the phase behaviour of the system. Orbits included in the CE model are selected with a genetic algorithm based on the Distributed Evolutionary Algorithm in Python (DEAP) suite [49]. A population of 100 individuals, each starting with five random-selected orbits, evolves for 20 generations. The best 50 models are selected from five repetitions of the evolution process. The evolution is driven by the cross-validation score of each individual, computed using the ten-split K-fold algorithm as implemented in Scikit-learn [50]. In order to favour low-complexity models with fewer orbits φ, a penalty p(c) = γΣ c is added to the cross-validation score of each individual c. Σ c denotes all the cluster functions defining the model c, i.e. all the orbits φ associated with non-null effective cluster interaction J. A value γ = 1 × 10 −6 has been found to yield a good compromise between reducing the number of orbits in the selected models and retaining satisfying accuracy.

DATA AND CODE AVAILABILITY STATEMENTS
Metastability metric data and the code used to generated it are included in this published article as supplementary information files as JSON database and NumPy binary files for the former and Jupyter notebooks/Python3 scripts for the latter. CE data and code used in this study are available from the corresponding author on reasonable request.