Origins of bandgap bowing character in the common-anion transition-metal-dichalcogenide ternary alloyed monolayer: ab initio investigation

Density functional theory is employed to investigate the origins of bandgap bowing character in transition-metal-dichalcogenide ternary alloyed monolayers (TMD-MLs). The available experimental photoluminescence (PL) data in literature have confirmed the existence of bowing character in the common-anion ternary alloys (e.g. Mo1−x W x S2) and its complete absence in the common-cation ternary alloys (e.g. MoS2(1−x)Se2x ). Our theoretical modeling of bandgap energy versus alloy composition, Egx , in these respective alloys have yielded trends and bowing parameters in excellent agreement with the available PL data (i.e. B = 0.26 eV and zero, respectively). Calculated band structures showed that the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) states in TMD-ML to be fully attributed to the metal atoms and to follow the symmetry of the irreducible representations A 1′ (singlet dz2 state) and E′ (doublet of dx2−y2 and d xy states) of the point group D 3h , respectively. Consequently, in case of common-cation TMD-ML alloys, Egx is linear and the bowing is absent. Whereas, in case of common-anion TMD-ML alloys, Egx is quadratic and the bowing is present because of the existence of competition between the cations (i.e. metal atoms) in contributing to HOMO/LUMO states. Our theoretical findings are corroborated with the available experimental data and have direct impact in TMD-based photonic nano-device applications. Video Abstract: Origins of bandgap bowing character in the common-anion transition-metal-dichalcogenide ternary alloyed monolayer: ab initio investigation


Introduction
The breakthrough discovery of graphene not only rewarded the inventors (Andre Geim and Konstantin Novoselov) Nobel Prize in Physics in 2010 but rather paved the way for the appearance of a new horizon of the physics of low dimensional 2D systems [1][2][3]. Triumphs in the syntheses of monolayers (MLs) with various compositions and applications occurred thereafter, such as silicene [4,5], germanene [6,7], BN [8,9], ZnO [10,11], and transition metal di-chalcogenides (TMD) [12,13]. The utilized growth methods customarily were at the level of chemical vapor deposition (CVD) and ranging up to the state-of-the-art molecular beam epitaxy [4][5][6][7][8][9][10][11][12][13]. Worthwhile, the achievement of reduced dimensionality cemented the way for a new revolution in materials science taking place in 2D systems in both fundamental science and technological applications. Namely, in their bulk 3D structures, TMDs possess indirect bandgap and by achieving the fabrication of MLs having direct bandgap [13], the interest in them completely changed and diversity of applications beyond photonics appeared thereafter.
In their bulk form, TMDs (denoted as MX 2 , M = Mo, W stands for metal atom and X = S, Se, Te stands for chalcogenide atom) belong to van der Waals (vdW) layered structures; where multilayers are stacked and maintained cohesively via vdW interactions. Each ML is formed by M-X covalent bonding with partial ionic character making like honeycomb triangular lattice (if seen from 'c-axis' direction perspective). Seen from side view, the structure might be described as a metal 'M' layer sandwiched between two chalcogenide 'X' layers; and the thickness of the ML is about 3 Å whereas the distance between adjacent neighboring chalcogenide layers is also about 3 Å. Both bulk and ML of TMDs are lacking the inversion symmetry but both possess horizontal mirror symmetry and belong to the crystalline point group of D 3h [14]. TMD-MLs have properties distinct from graphene; namely, (1) TMD-MLs have direct bandgap transitions of energies laying within the visible spectrum making them suitable for semiconducting applications such as transistors [15], optoelectronic devices such as light-emitting diodes (LEDs) and light detectors [16,17]; (2) lacking the inversion symmetry with charge accumulation at K-corners has opened a new field of physics called valleytronics [18]; (3) TMD-MLs possess very strong spin-orbit coupling (SOC) and splitting, which can be explored in spintronics' applications [19]; and (4) TMD-MLs can be used with other 2D materials (e.g. graphene and BN) to fabricate so called vdW hetero-structures to be used in solar cells, LEDs, photodetectors, fuel cells, photocatalytic and sensing devices [20][21][22].
In the experimental side, in 2010, two experimental groups led by Splendiani et al [16] and Mak et al [23] were among the first to achieve the fabrication of TMD-MLs based on MoS 2 , and to detect their embryonic photoluminescence (PL), which was taken as an indication of the crossover from indirect (E g = 1.26 eV) to direct bandgap (E g = 1.90 eV) transitions along the attempts of growing thinner and thinner MLs [23]. The synthesis method was based on the micro-exfoliation technique, similar to the one used by the graphene inventors, because the interactions between layers is relatively weaker than the cohesive interaction of ML. Thereafter, further development in growth techniques suggested various methods for bandgap engineering via either allowing the ML to persist in maintaining the PL or using the ML in hetero-structures [24,25]. For instance, Wang and co-workers [26] reported their success in fabricating high-quality Mo 1−x W x S 2 homogeneously-alloyed MLs with tunable Mo/W ratio using low-pressure chemical vapor-phase deposition (LP-CVD) method. The homogeneity of the Mo/W mixture in the synthesized Mo 1−x W x S 2 MLs was confirmed using both Raman spectroscopy and the high-resolution transmission electron microscopy. Photoluminescence measurements showed the optical bandgap to be strongly dependent on composition 'x'. The PL data demonstrated a pronounced bowing character with a bowing parameter of about B ∼ = 0.26 eV [26]. Similar bowing parameter was reported by other experimental groups [27,28]; for instance, Tan and co-workers reported B = 0.25 ± 0.04 eV for similar alloyed Mo 1−x W x S 2 MLs [27]. These latter researchers have further focused on studying the effects of disorder. Their modeling results demonstrated that the valence-band maximum (VBM) and the conduction-band minimum (CBM) to be controlled by the composition (i.e. E g = E g (x)) and more importantly that CBM to be sensitive to the disorder (i.e. Stoichiometry) in the alloy. The ordered phase was found to correspond to the minimum energy gap with CBM dominated by Mo contributions [27].
The bowing character persists to exist in the common-anion TMD-ML alloys. For instance, Tongay and co-workers reported their successful growth and PL data of Mo 1−x W x Se 2 ML alloys. Their reported bowing parameter B = 0.14 eV was further taken as an indicator of the sensitivity of optical gap on composition and thus revealing somehow the high quality of the ML alloys. Recently, Kopaczek and co-workers [30] reported their results of combined experimental and theoretical studies on Mo 1−x W x Se 2 ML alloys and corroborated the findings of a bowing parameter of about 0.14 ± 0.03 eV. As the band edges (i.e. VBM and CBM) being controlled by the metal atoms, the bowing character was observed to become vanishingly small in case of the common-cation TMD-ML ternary alloys, such as MoS 2(1−x) Se 2x [31,32] and WS 2(1−x) Se 2x ML alloys [33,34]. In short-term, the experimental data showed the evidence of existence of the bowing character in the common-anion ternary TMD-ML alloys [25][26][27][28][29][30] and the absence of this character in the case of common-cation ternary TMD-ML alloys [31][32][33][34]. Most of the previously mentioned experiments have presented their data while been corroborated with the ab initio simulations. Yet, the origins of the existence/absence of the bowing character are still not fully understood and deserve further investigation which will be the scope of the present work.
In contrast to the bowing behaviors in TMD-ML alloys, in III-V and II-VI semiconducting alloys, the bond is more covalent in character and the valence band (VB) of the crystal is predominantly contributed from the anion atoms. So, in the case of common-cation ternary alloys (e.g. CdSe x Te 1−x and ZnSe 1−x Te x alloys), the VBM is mainly attributed to the anion contributions, and based on electronegativity character perspective, the anion atoms compete in trapping the electronic charge and such competition would trigger the bandgap bowing effect [35]. Whereas in the case of common-anion ternary compound-semiconductor alloys (e.g. Cd 1−x Zn x Te and Cd 1−x Zn x Se alloys) the lack of electronegativity competition would result in the vanishingly small bowing effect or perhaps its complete disappearance [36]. So, the scope of the present investigation is to utilize first-principles calculations to search for the origins for existence/absence of bowing effects in TMD-ML alloys and assess the main factors causing and controlling these effects.
The present paper is organized as follows. Next section gives detailed description of the computational method used in our work. Section 3 elaborates a discussion of the results. It includes also a group theoretical analysis of the splitting of metal d-states under the ligand field of the crystalline point group (i.e. which is trigonal prismatic D 3h ). The last section summarizes our main findings.

Computational method
We employed the Spanish initiative for electronic simulations with thousand atoms (SIESTA), whose efficiency stems from the use of strictly-localized atomic orbitals as basis set in order to allow the method to deal with large systems [37,38]. Actually, we used SIESTA which is incorporated within (atomistic tool-kit (ATK) package. The accuracy and cost can be adjusted in the code to range from quick exploratory calculations to highly accurate simulations matching the quality of other approaches, such as plane-wave methods. The basis set utilizes the atomic orbitals composed of the product of a radial part and spherical harmonics. So, the basis set becomes smaller than those of linearized augmented plane wave method and a better description of finite systems such as molecules, clusters and 2D MLs is achieved in avoiding the cost of describing vacuum. The Hamiltonian has a linear scaling with the system size so it possesses the ability to deal with very large systems of a size of thousands of atoms. Norm-conserving pseudopotentials are used to represent the electron-ion interactions and to focus on valence electrons filling the atomic outer shells. Periodic boundary conditions (PBCs) are implemented within the framework of SIESTA in the solution of Poisson equation and treatment of Hartree's term of the single-particle Hamiltonian.
We choose a supercell consisting of 12 atoms in case of TMD-ML MX 2 (i.e. 2 × 2 primitive cells containing: four metal and eight chalcogenide atoms) and of 24 atoms in case of TMD bulk structure to use for comparative study (i.e. 2 × 2 × 1 primitive cells containing eight metal and 16 chalcogenide atoms) with PBCs applied in all directions. The size of this supercell would allow one to study the case of the common-anion ternary TMD-ML alloys (i.e. Mo 1−x W x S 2 , with discrete x variable 0 x 1 and step Δx = 1/8 = 0.125, while the anion atom 'S' is kept unchanged), with seven different x compositions of tungsten content. It would also allow the study of the common-cation ternary TMD-ML alloys (i.e. MoS 2(1−x) Se 2x , with discrete x variable 0 x 1 and step Δx = 1/16 = 0.0625, while the cation atom 'Mo' is kept unchanged), with 15 different x compositions of selenium content. We have also considered various possibilities of alloying configurations/realizations. So, the effect of disorder (stoichiometry) will be taken into account in our computation. Raising the size of the supercell would not be necessary as the configurational space of disorder would raise enormously [26,[39][40][41] and as stated by Tan and co-workers [26], the CBM is the only band edge sensitive to alloy heterogeneity/disorder in case of TMD-ML alloys (e.g. specifically the Mo 1−x W x S 2 alloys). The same authors showed that the minimum bandgap energy should correspond to the homogenous/ordered alloys. So, our current supercell size should be sufficient to study both the effects of composition and disorder and pursue our study of variation of bandgap energy versus composition in case of ordered alloys in modeling the PL data due to high-quality grown (homogeneous) alloyed MLs.
Furthermore, it is worth to mention about the energy cut-off used within the framework of SIESTA. The energy cut-off E cut = 150 Ry (i.e. 2040.75 eV), which is the default value, does correspond to the Fourier and real-space grids on which the Poisson equation would be solved (including the real-space electron density and the effective potential descriptions). Of course, the higher E cut would give finer real-space grid and consequently better accuracy. For the exchange-correlation functional, we used the general gradient approximation, parameterized by Perdew-Burke-Ernzerhof [42]. All atoms were allowed to relax until the Hellmann-Feynman forces on atoms and the total energy of the system become less than 0.05 eV Å −1 and 0.1 meV, respectively. The pressure tolerance is set to 0.1 GPa. For the Brillouin-zone sampling, we used the Monkhorst-Pack technique [43]. Grids of K-mesh = 3 × 3 × 1 (for ML) and 3 × 3 × 2 (for bulk) are used for density functional theory (DFT) iterations in the tasks of both the atomic relaxations and the calculations of the partial and total density of states (PDOS, TDOS), respectively. In SIESTA, the statistics of DOS is carried out with a smearing or a thermal broadening of k B T = 0.025 eV at room temperature to smoothen the DOS.
It should be emphasized that our systems are confirmed to be diamagnetic so that we excluded the spin-polarization in our calculations. Our simulations should include the atomic relaxations and the calculations of total energy, PDOS/TDOS, bands, the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) states. In addition, we have also shown group theoretical study of the splitting of metal-atom d levels under the ligand field of the trigonal prismatic point group (D 3h ), where its corresponding reducible representation gets decomposed into the irreducible representation of the point group (i.e. five-fold d level = A 1 + E + E ). The results are going to be discussed in the next section.  Figure 1 shows both the experimental PL data and the theoretical bandgap energy results versus composition of TMD-ML ternary alloys of the common-cation and common-anion cases (i.e. MoS 2(1−x) Se 2x and Mo 1−x W x S 2 , respectively). It is important to elaborate about the experimental syntheses of such alloys and their structural and optical characterizations.

Modeling of PL data
In the synthesis of common-cation TMD-ternary alloys, Bay and co-workers used CVD growth technique to synthesis MoS 2(1−x) Se 2x alloyed ML on silica, while having full control on growth parameters such as gas-flow rates, substrate temperature and precursor concentration (i.e. S/Se ratio) [32]. Besides, they had a full control on x values within the range x = 0-1 to yield bandgap energies tailing between 1.82 eV (MoS 2 ) and 1.55 eV (MoSe 2 ), where the size of the grown ML flakes has been measured to be as large as 150 μm on SiO 2 (buffer of 300 nm)/Si substrates. Actually, they reported two growth regimes: (i) S-rich alloys were grown by controlling S to Se ratio at 750 • C; (ii) Se-rich alloys were synthesized at higher growth temperatures (900 • C) by using Se rich vapor conditions, which was attributed to the chalcogen exchange mechanism. Concerning the structural properties, the full width at half maximum of PL spectra of the flakes were reported to range between 19 and 37 nm along with uniform spatial PL emission intensities were used as indicating the high crystal quality [32]. Again, CVD technique was used in the synthesis of common-anion TMD-ternary alloys. Wang and co-workers developed a low-pressure CVD (LP-CVD) method to grow MLs of Mo 1−x W x S 2 (x = 0-1) 2D crystals with a wide range of Mo/W ratios [26]. PL measurements showed the strong correlation of bandgap energy on the alloy composition (x) and clearly demonstrated the bowing effect. This in turn was considered as signature or indication of the achievement of the successful growth of uniform and crystalline alloyed samples [26].
It is very important to start the fitting of the experimental data by an assessment of the pristine MLs' photonic bandgap theoretical values. In this context, there is a strong correlation between the band and atomic structures so that the bandgap is sensitive to the lattice parameters. Furthermore, from the perspective of DFT calculation, the achievement of reliable relaxed geometrical structure has always been taken as a measure to the quality of the used pseudopotentials in the calculation. For instance, in our present work, for relaxed 3D bulk structures MoS 2 and WS 2 , we obtained the lattice parameters (a = 3.160 Å, c = 12.295 Å) and (a = 3.153 Å, c = 12.323 Å); which are in excellent agreements with the experimental lattice parameters (a = 3.15 Å and c = 12.30 Å) [44] and (a = 3.153 Å, c = 12.323 Å) [44], respectively. Furthermore, we have interest in MLs as they stand relevant for photonics applications. The MLs possess the symmetry of a triangular lattice with a single lattice parameter a = b (i.e. honeycomb structure as seen from top view perspective). Table 1 summarizes the lattice parameters and bandgap energies for the three pristine MLs of our interest and compare them to the available experimental data [45][46][47][48]. Our theoretical results show that E g = 1.79 eV, 1.46 eV, and 1.90 eV corresponding to pristine MLs MoS 2 , MoSe 2 , and WS 2 , respectively. Taking into account the usual underestimation of the DFT to the bandgap energies in semiconductors, then the respective error bars of 1.6%, 6.4% and 6.4% with respect to the experimental values to be modeled [26,32], would be considerably small to be deliberated as good achievements. These achievements are attributed to the improvements in our appropriate use of hybrid functionals in the approximation of exchange-correlation and the pseudopotentials.
The experimental PL data of Bay an co-workers [32] and Wang and co-workers [26], corresponding to the single MLs of MoS 2(1−x) Se 2x and Mo 1−x W x S 2 alloys, are presented in figure 1 (i.e. panel figures 1(a) and (b), respectively). These two respective experimental data clearly show the absence of bandgap bowing effect in the common-cation TMD-ML ternary alloys and its pronounced existence in the case of common-anion TMD-ML ternary alloys. The least-square fitting of PL data of Mo(SSe) 2 alloys yielded linear function of E g versus alloy composition (i.e. E g = −0.33x + 1.83); whereas the fitting of PL data of MoWS 2 alloys yielded a quadratic function of E g versus alloy composition (i.e. E g = 0.26x 2 − 0.12x + 1.89). This is just the opposite to what usually was observed in compound-semiconductor ternary alloys made of either III-V or II-VI ternary alloys. So, these experimental PL data induced in us a great curiosity and motivation to theoretically investigate the origins and the factors controlling the bowing effect in TMD-ML alloys.
Our fitting results for the same two alloys Mo(SSe) 2 figure 1(c)). This, in turn, is found to be consistent with the attempt of quadratic fitting of PL data of panel figure 1(a) (which showed: E g = 0.05x 2 − 0.36x + 1.8), yielding in a similar way a vanishingly small bowing parameter of B = 0.05 eV. So, our sampling of configurational space yielded theoretical results corroborating the experimental PL data of Mo(SSe) 2 alloys. As a matter of fact, the alloys which yield the optimum bandgap energy do correspond to the homogenous alloys in consistency with what was reported by Tan and co-workers [27]; whereas those which have tendency to Se clustering would yield E g of higher estimations. To further investigate the absence of bowing effect is Mo(SSe) 2 alloys, we decided to take a case of homogenous alloy with an alloy content x = 0.5 to inspect the electronic structure and optical properties, here below.  Figure 2 shows the plots of amplitudes of both the highest orbital molecular orbital (HOMO) and the LUMO Eigen-functions for both MoS 2(1−x) Se 2x (with x = 0, 0.5, and 1) and Mo 1−x W x S 2 (with x = 0, 0.25, and 1) ML alloys, in panel figures 2(a) and (b), respectively.

HOMO and LUMO states
It is important to emphasize that the ML of TMD (i.e. denoted by MX 2 , M = metal such as Mo and W; X = chalcogenide such as S, Se and Te) should exhibit one of two polymorphs: (i) trigonal prismatic phase, which belongs to the D 3h point group (i.e. it is referred as 1T or 2H); or (ii) octahedral phase, which belongs to the O h point group (i.e. it is referred as 1H). In both 1H and 1T phases, the non-bonding d bands of the TMDs are located within the gap between the bonding (σ) and antibonding (σ * ) bands of M-X bonds [49] (i.e. they might originate sometimes form non-bonding states). Under the ligand field, the d orbitals of M atom split in different forms depending on both the coordination environment and the d-electron count: (1) in case of octahedrally-coordinated transition metal centers of point group symmetry O h , the d orbitals split into two groups of ascending energy order as follows: triply-degenerate group t 2g (containing d xy , d yz , and d zx orbitals) and doubly-degenerate group e g (containing d z 2 , and d x 2 −y 2 orbitals).
(2) In case of transition metal centers with trigonal prismatic coordination and having point group symmetry D 3h (like our present case), the d orbitals split into three groups in the following ascending energy order: a singlet A 1 (composed of d z 2 ), a doublet E (composed of d x 2 −y 2 and d xy ), and a doublet E (composed of d yz and d zx ), with a sizable energy gap of about E g ∼ 1 eV between the first two groups of orbitals [49]. The electronic properties of TMDs are mainly controlled by the metal atom and very dependent on the filling of the non-bonding d bands from group 4 to group 10 (i.e. when d orbitals are partially filled TMDs exhibit metallic conductivity; whereas when d orbitals become fully filled TMDs become semiconducting). On the other hand, the effect of the chalcogenide atoms on the electronic structures is minor (e.g. in case of semiconducting TMD, scaling in an increasing size like S to Se to Te would gradually reduce the band-gap energy) [49,50].
The preferred phase adopted by TMD is very much dependent on the d-electron count of the transition metal 'M'. Most of TMDs crystallize in octahedral phase in exception of two groups. TMDs with M atom from group 5 (i.e. so named d 1 , because five electrons are used in sp 3 d 1 hybridization) can crystallize in both octahedral and trigonal prismatic phases. Whereas, TMDs with M atom from group 6 (i.e. so named d 2 , because six electrons are used in sp 3 d 2 hybridization) are found to crystallize in just trigonal prismatic phase (e.g. M = Mo and W), which is our present study case.
It is worth mentioning that figure 2 displays both the relaxed atomic structures and the HOMO and LUMO eigen-states. Two views are presented (i.e. a top view and a side view in which the xy-plane looks vertical). All the relaxed structures confirmed their stability in the trigonal prismatic phase 2H (i.e. of point group symmetry 'D 3h '). This crystalline point group is lacking the inversion symmetry operation and has lower symmetry than the octahedral group 'O h '. Under the point group D 3h , the five-fold degeneracy of d-states of the metal atom should be lifted and should exhibit a strong splitting into three groups of states in the following ascending energy order: a singlet A 1 (d z 2 orbital), a doublet E (d x 2 −y 2 and d xy orbitals), and a doublet E (d yz and d zx orbitals). In TMDs, there exists a gap between A 1 and E groups. By inspecting the HOMO/LUMO states presented in figure 2, one notice two main trends: (1) trend #1: the common feature between all LUMO states is that they are all confined to the metal-plane as they are likely attributed to the E -group (i.e. d x 2 −y 2 and d xy orbitals). It is evident that the LUMO state has the characteristics to be a non-bonding state as it is not distributed/shared along the bonds but rather being confined and curdled on the metal atoms. This aspect will be further corroborated by the results of PDOS in next sub-section. Furthermore, one exception, that should be emphasized, is that in case of Mo 1−x W x S 2 with x = 0.25 (shown in figure 2(b)), the LUMO state does not have contribution from W atom at all. This means that the conduction band (CB) edge of MoS 2 is lower in energy than that of WS 2 (i.e. E MoS 2 C < E WS 2 C ) and the LUMO is accommodated solely by Mo atoms in case of MoWS 2 alloys. Such competition should be considered as one of the main reasons triggering or causing the occurrence of bandgap bowing behavior in common-anion TMD-ML ternary alloys such as in the case of MoWS 2 alloys. (2) Trend #2: the HOMO state in all panels of figure 2 has the characteristics of non-bonding state, as being confined to the metal plane but to follow the symmetry property of the A 1 irreducible representation (composed of a single component of d z 2 orbital). These results are corroborating the group theoretical predictions that d states split under the crystal field of point group D 3h with a singlet A 1 laying the lowest, a doublet E laying in the middle and E laying being the highest. So, our present figure 2 is in excellent agreement with the predictions of the group theory (see subsection 3.5 for more details). Figure 3 displays the results of PDOS and TDOS for the two alloys of concern: (a) MoS 2(1−x) Se 2x ; and (b) Mo 1−x W x S 2 with x = 0, 0.5, 1 and x = 0, 0.25, 1, respectively. The color of PDOS curves are as follows: Mo (in red), S (in blue), Se (in green), W (in orange); and TDOS curves (in black). In all panels, although the ratio of metal to chalcogenide atoms should be 1:2, the states near band-edges within 1 eV (i.e. in the ranges E V − 1.0 E E V and E C E E C + 1.0) are more dominated by contributions from metal atoms. Besides, both VB and CB edges are fully attributed to metal atoms, in consistency with the previously shown results in figure 2 of HOMO and LUMO states. This makes the metal atoms have a control on the band-gap energy. Consequently, one can deduce the following behaviors: (a) in case of common-cation TMD-ML ternary alloys (e.g. MoS 2(1−x) Se 2x alloys), there is no competition in the contribution to the HOMO/LUMO states as these latter are made solely from the same metal (i.e. Mo atoms). So, the band-gap energy E g (x) should vary linearly versus the anion composition x. (b) On the other hand, in case of common-anion TDM-ML ternary alloys (e.g. Mo 1−x W x S 2 alloys), the competition in the construction of HOMO/LUMO states exists. Such competition between the two metal atoms would be the reason triggering the bowing effect in the bandgap energy characteristics (i.e. E g (x) = (1 − x)E g (MoS 2 ) + xE g (WS 2 ) − Bx(1 − x), which is a quadratic variation, where B is the 'bowing parameter'). Figure 4 shows the band structures of both (a) the common-cation TMD-ML ternary MoS 2(1−x) Se 2x alloys, and (b) the common-anion TMD-ML ternary Mo 1−x W x S 2 alloys. All the structures yield direct energy bandgaps at the special K-point in the Brillouin zone. Similar to figure 3, we focus on the energy interval [−3.0, +3.0] eV (i.e. around Fermi level, where E F = 0). By inspecting the band structures, one can make   the following remarks: (i) without the inclusion of spin-polarization, the bands possess degeneracies not more than two at the high-symmetry points in the Brillouin zone, which follows the sizes of the irreducible representations of the D 3h point group, shown in table 3. (ii) The highest VB is clearly flat especially near the BZ center (i.e. at Γ-point). The top of VB is located at K-point of the BZ and the difference between top of VB at K-point and Γ-point is very small (i.e. of the order of 'SOC' energy) (see table 2). The bands shown are due to the ab initio calculations without spin. Yet, the inclusion of spin-polarization results have shown that the order of the bands remain the same and furthermore these materials are confirmed to be diamagnetic. (iii) The lowest CB lays at K-point for all these three compounds to yield direct-band gap transitions. Actually, such direct bandgap transition was discovered for 1 ML systems in TMD in just 2010 [16,23]. So that these materials, thereafter, attracted enormous interest as becoming relevant for photonic applications. (iv) It is well-known trend that E g decreases with increasing size of either cation or anion atoms in all semiconductor compounds (i.e. see the cases of bulk MoS 2 , MoSe 2 , WS 2 having E g = 1.79 eV, 1.46 eV, 1.90 eV, respectively). Yet the trend here is not quite alike III-V and II-VI compound semiconductors because of the strength of the spin-orbit interaction is fairly large.

Band structures
In table 2, the SOC is reported to rather increase with the size of either anion or cation atoms (i.e. SOC = 150-430 meV). The interplay between crystal-field splitting (see next sub-section) and the SOC would produce the energy gap. One should further emphasize that Rashba and Sheka discovered the SOC effects in bulk crystals and low-dimensional systems in 1959 [57]. Thereafter, in 1984, Rashba and Bychkov discovered stronger SOC in 2D systems as enhanced by the crystalline anisotropy (also called Bychkov-Rashba effect) [57]. Such effect to be added to make SOC in TMD ML more valuable but yet to remain less than crystal splitting effect. Table 3. Character table of D 3h , where the symmetry operations and the irreducible representations are denoted in Schoenflies' notation.

Crystal-field splitting of d levels
TMD (MX 2 ) based on metal atom (M = Mo, W) and two chalcogenide atoms (X = S, Se, Te) do crystalize into 'trigonal prismatic phase' (i.e. from top-view, the bulk structure looks like honeycomb/graphitic with AB stacking where M-and X-bases are alternating from layer to layer along c-axis). The bulk MX 2 of 'trigonal-prismatic phase' (2H) has lower symmetry than octahedral as it is missing the inversion symmetry so that the structure belongs to non-symmorphic crystallographic point group. The bulk MX 2 structure belongs to the space group D 4 6h (P6 3 /mmc, #194) [58], whereas, the ML structure had a furtherly reduced symmetry to belong to the space group D 1 3h (P6m2, #187) [58]. The trigonal prismatic phase has a lower symmetry than the octahedral group O h and, consequently, the degeneracy of the triplet group (T 2g ) would be furtherly lifted by its splitting into two sub-groups (i.e. a singlet state of irreducible representation A 1 and doublet states of irreducible representation E ). Furthermore, as the ML of MX 2 becomes quasi-two dimensional, more possible occurrence of more fine splitting would be plausible as it exhibits stronger SOC at reduced dimensionality [59].
The aim of this sub-section is to demonstrate the effects of both the reduced symmetry and the reduced dimensionality in TMD MX 2 ML. Table 3 shows the character table of the irreducible representations of the crystallographic point group D 3h [60][61][62] corresponding to the trigonal prismatic structure. The character table contains six classes of symmetry operations and consequently six irreducible representations, which are denoted in Schoenflies' notation. The order of the group is h = 12 (i.e. total number of symmetry operations). Table 4 shows the character representations of the d-states (i.e. atomic angular momentum L = 2) under the symmetry operations of the point group D 3h . Recall that the electron would not move inside the atomic orbital when the atom alters the hybridization to accommodate the crystalline structure, but rather it moves in the solid-state eigen-state. The spherical harmonics alter the transformation caused by the crystal field. We further emphasize that the splitting caused by the crystal field is much larger than that caused by the 'SOC' (which was elaborated in the previous sub-section).
Among the well-known properties of the character table is the normalization condition: where summation is carried out over all the symmetry classes (symbolized 'R'); g i is the number of symmetry operations in the class R, χ i (R) is the character value of the irreducible representation at the symmetry operation R; h is the order of the point group (i.e. h = 12 in our present case); and δ i,j is the Kronecker symbol. The reducible representation Γ(L = 2), corresponding to five-fold degenerate d-states of the metal (e.g. Mo or W atom), should be expanded under the action of the crystal ligand-field in terms of the irreducible representations: where Γ is the reducible representation and Γ j stands for the irreducible representations shown in table 3 (e.g. A 1 , A 2 , E , etc), and n j is an integer usually 1 or 0. There exists a formula to find n j as follows: Where χ j (R) and χ (R) stand for the irreducible and reducible representations shown in tables 3 and 4, respectively. By applying this formula (3), one can write the decomposition of d-states as follows: So the results yield the following decomposition: In our present work, it is confirmed that A 1 has the lowest energy among these three split groups. The ascending order in energy is determined to be: A 1 (singlet m = 0 state: d z 2 ), E (doublet m = ±2 states: d x 2 −y 2 and d xy ), and E (doublet m = ±1 states: d yz and d zx ). Our investigation of HOMO/LUMO states presented in subsection 3.2 have shown that A 1 and E are represented by the HOMO and LUMO states, respectively. This is consistent with literature [56,63]. If we consider the VB as an energy reference, then: E HOMO = E(A 1 ) = 0, so E LUMO = E(E ) = E g = 1.79 eV. If we assume that the splitting occurred in the proportion of ΔE 1 = E(E ) -E(A 1 ) and ΔE 2 = E(E ) -E(E ), so within the framework of perturbation theory one would expect: ΔE 2 = ΔE 1 /2. In this case, one can predict that the original five-fold degenerate atomic d-states would have an energy E d = E g = 1.79 eV and that E(E ) = 1.5 × E g = 2.685 eV. The energy diagram presented in figure 5 illustrate the energy splitting of d orbitals and make a comparison between In DOS plots, one should be able to identify VBM and CBM to originate from sulfur atoms in case of bulk and from Mo atoms in case of ML. four cases of: (i) isolated Mo atom as a start; (ii) under octahedral crystal field; (iii) under trigonal prismatic field but in case of 3D bulk; and (iv) under trigonal prismatic field but in case of 2D ML.

Nature of transitions in ML versus bulk
We have performed the calculations of the HOMO/LUMO states in case of pristine bulk MoS 2 (see figures 6(a) and (b)) to compare it with our current case of MoS 2 ML. The supercell of bulk MoS 2 contains two layers with a total number of 24 atoms (i.e. it is of size double the ML case). Both top view and side view are displayed for better clarity. If looked from top, the system is bi-layer with the hexagon of lower layer being rotated by π/3 about a z-axis passing through the center of the hexagon (i.e. the Mo atom looks on the top of S atom). So, for a better judgment, it is better to inspect the side views. Both HOMO/LUMO of bulk MoS 2 are found to be mainly constructed by anion states with mimic or none contribution from metal M. This is completely different from the ML case discussed in figure 2. The band structures and DOS are shown in figures 6(c) and (d) for both bulk and ML MoS 2 , respectively. The bandgap of bulk MoS 2 is smaller (E g = 0.88 eV) and indirect. Both our direct bandgap energy of ML E ML g = 1.79 eV and indirect bandgap energy of bulk E BULK g = 0.88 eV are consistent with ab initio results which was reported by Naumis [64]. The VBM becomes located at the BZ center (i.e. Γ-point) but the CBM remains located at the K-corner of the BZ. In case of bulk 3D, it seems that the anion atoms are the ones contributing to both VBM and CBM and the transition is indirect as likely to occur from S → S but the two anions might be located in two different but adjacent layers as the indirect transition usually requires the assistance of a phonon.
The experimentally observed two transitions (so-called A & B) [65,66] in TMD-ML can be understood to occur between localized states of metal atoms: (i) transition 'A' to occur from VB maximum of symmetry A 1 and CB minimum of SOC-split-down from E group; (ii) transition 'B' to occur from VB maximum of symmetry A 1 and CB state of the second SOC-split-up state of E group [60]. Our theoretical predictions are in favor of the occurrence of these transitions between such localized states with high transition probability and high recombination rate.

Conclusions
In contrast to the III-V and II-VI semiconductor alloys [35,36], the bowing character in TMD-ML alloys exists only in the of common-anion ternary alloys as the HOMO and LUMO states are fully controlled by the metal atoms. From theoretical point of view, spanning the configurational space of ordered/disordered alloys, we found that the optimum (minimum) bandgap energy should correspond to the ordered alloys. So, we focused on studying the variation of E g versus composition in the two cases of common-anion and common-cation TMD-ML ternary ordered alloys. The main findings of this work can be summarized as follows: (a) In the common-cation TMD-ML ternary alloys (MoS 2(1−x) Se 2x , with 0 x 1), the HOMO and LUMO states are found to be originated from non-bonding states residing on Mo layer of respective A 1 and E symmetry groups. Being controlled solely by the metal Mo atom, the bandgap energy E g is found to vary linearly with anion composition x; (b) In the common-anion TMD-ML ternary alloys (Mo 1−x W x S 2 , with 0 x 1), the two metal atoms compete in contribution to the electronic structures of HOMO and LUMO states. Such competition triggered the bowing effect and the bandgap energy E g is found to vary quadratically with cation composition x. Our theoretical estimation of the bowing parameter (B = 0.26 eV) is in good agreement with the experimentally reported values.
While the crystal-field splitting is the most important in our systems of TMD-ML alloys, the band structures can be posteriorly corrected by other parameters such as the many-body, the SOC and the excitonic effects to further understand the fine A and B excitonic transitions. As it stands, our work is expected to contribute to the understanding of the optical properties of TMD-ML alloys which should be of great interest in the optoelectronic devices' applications.