Atomic and electronic structure of two-dimensional Mo(1− x )W x S2 alloys

Alloying enables engineering of the electronic structure of semiconductors for optoelectronic applications. Due to their similar lattice parameters, the two-dimensional semiconducting transition metal dichalcogenides of the MoWSeS group (MX2 where M = Mo or W and X = S or Se) can be grown as high-quality materials with low defect concentrations. Here we investigate the atomic and electronic structure of Mo(1−x)W x S2 alloys using a combination of high-resolution experimental techniques and simulations. Analysis of the Mo and W atomic positions in these alloys, grown by chemical vapour transport, shows that they are randomly distributed, consistent with Monte Carlo simulations that use interaction energies determined from first-principles calculations. Electronic structure parameters are directly determined from angle resolved photoemission spectroscopy measurements. These show that the spin–orbit splitting at the valence band edge increases linearly with W content from MoS2 to WS2, in agreement with linear-scaling density functional theory predictions. The spin–orbit splitting at the conduction band edge is predicted to reduce to zero at intermediate compositions. Despite this, polarisation-resolved photoluminescence spectra on monolayer Mo0.5W0.5S2 show significant circular dichroism, indicating that spin-valley locking is retained. These results demonstrate that alloying is an important tool for controlling the electronic structure of MX2 for spintronic and valleytronic applications.


Introduction
Semiconducting transition metal dichalcogenide monolayers, such as the MoWSeS group (MX 2 where M = Mo or W and X = S or Se), are of great interest both for their practical applications and for the fundamental science that can be studied in them [1]. Much of the interest has focused on the optical properties of MX 2 : they exhibit a transition from indirect gap in bulk to direct gap in monolayers, and strong Coulomb interactions lead to high exciton binding energies [2]. In the monolayers, spin-orbit coupling (SOC) splits both the valence band maximum (VBM) and conduction band minimum (CBM) to give spin-polarised bands at the Brillouin zone corners [3,4]. This leads to spin-valley locking [5][6][7][8], with the potential for optically generated spin-polarised currents, and the resultant interest in valleytronics [9] and spintronics [10] in these materials.
Band engineering through alloying has been essential to the development of III-V semiconductors, giving control over band parameters such as band gaps and band alignments, leading to their application in commercial optoelectronics. Following this, theoretical and experimental works have shown that for MX 2 alloying allows not only continuous tuning of the band gap over a large range, from 1.55 to 1.99 eV [11][12][13][14][15][16][17][18], but also control over the magnitude of SOC in both the conduction and valence band [13,14,19].
Combined, these band parameters have dramatic influences on the optical and electrical properties of MX 2 . However, due to the difficulty in accurately determining the single-particle electronic structure, there is still uncertainty over how the electronic structure changes with alloying. The high exciton binding energies complicate the determination of electronic structure parameters from conventional optical spectroscopy measurements, such as photoluminescence (PL). This leaves some outstanding questions: for example the SOC is predicted to change linearly with composition, but inferring the SOC from the energy difference between A and B excitons suggested a non-linear bowing in Mo (1−x) W x Se 2 monolayers [20].
With the comparatively small difference in lattice parameters, high-quality MoWSeS alloys are relatively easy to synthesise, either through direct growth of monolayers by chemical or physical vapour deposition [11,17,19,[21][22][23], or by mechanical exfoliation of bulk alloy crystals grown by chemical vapour transport [14,[24][25][26]. In such alloys, the local atomic arrangement can play an important role in determining the resultant properties of the material. Both atomically random [14,19,25] and ordered [27][28][29] alloying have previously been observed in transition metal dichalcogenide monolayers and theoretical studies have suggested the possibility of order-disorder phase transitions [30,31]. Study of the ordered alloy phases revealed sensitivity of the bandstructure at both the VBM and the CBM to composition, but specifically also sensitivity at the CBM to ordering [31].
In this report we study exfoliated flakes of chemical vapour transport (CVT) grown Mo 1−x W x S 2 alloys. By comparing quantitative analysis of the atomic structure determined by scanning transmission electron microscopy (STEM) to Monte Carlo simulations, we reveal that the atomic arrangements are consistent with those expected from thermodynamic considerations. Starting from these atomic structure models, and using large unit cells to avoid simulation artefacts, we apply linear scaling density functional theory (LS-DFT) to predict electronic structure changes with alloying. Comparison with angle-resolved photoemission spectroscopy (ARPES) valence band measurements shows that the SOC at the VBM does indeed scale linearly with stoichiometry (x). The predictions indicate that at intermediate compositions the SOC at the CBM should be less than the disorder potential in the alloys. Despite this, polarisation-resolved PL measurements of a Mo 0.5 W 0.5 S 2 monolayer show that spin-valley locking is retained and hence that such alloys have a promising future in spintronic applications.

Crystal growth and exfoliation
Single crystals were synthesised by CVT in a two-step process [32]. First, Mo (purity 99.9%, Sigma-Aldrich), W (purity 99.9%, Sigma-Aldrich) and S (purity 99.9%, Sigma-Aldrich) element powders were mixed stoichiometrically into an ampoule. The ampoule was pumped down to a pressure of 10 −6 mbar and sealed. The mixture was heated to 1000 • C for 3 d to form Mo 1−x W x S 2 powder. In the second step, crystals were grown from this powder. The synthesised compounds were transferred to a new quartz ampoule with larger diameter and mixed with the transport agent, I 2 (10 mg cm −3 of the ampoule volume). To keep the I 2 stable, the ampoule was evacuated to 10 −6 mbar in ice and sealed. The ampoule was placed into a three-zone furnace as shown in supplementary material section 1 (suppmat S1) (available online at stacks.iop.org/ JPMATER/4/025004/mmedia), the charge zone was kept at 1050 • C for 20 d with the growth zone at 950 • C. After cooling to room temperature, the ampoule was opened in air and single crystals were collected at the growth end as shown in the schematic figure in suppmat S1.
For PL measurements, the as-grown crystals were mechanically exfoliated onto a Si/SiO 2 wafer using scotch tape. Prior to peeling off, the substrates were heated for a few minutes on a hot plate at 130 • C [33]. Monolayer flakes were identified by their optical contrast, by their Raman spectra, and by their PL emission. PL and Raman spectra were acquired at room temperature with a confocal microspectrometer (Renishaw inVia Reflex Raman microscopes) with 442 nm and 532 nm excitation.

Compositional analysis
Crystal composition was analysed by a combination of x-ray photoelectron spectroscopy (XPS) and energy dispersive x-ray analysis (EDX), further details are given in suppmat S2. EDX showed uniform composition across individual crystals, indicating homogeneity at the microscale. The compositions determined by XPS and EDX were consistent with each other and with those measured from analysis of the atomic resolution STEM images, indicating homogeneity also at the nanoscale.

Atomic resolution imaging
Flakes were mechanically exfoliated onto chemical vapour deposition (CVD) grown graphene on copper [34], then a layer of CVD graphene was wet transferred on top [35]. The bottom layer of copper was etched away using ammonia persulphate, the stack gently washed with deionised water, and then the graphene encapsulated Mo 1−x W x S 2 monolayer flakes were transferred to a transmission electron microscopy (TEM) support grid. STEM analysis was performed in a JEOL ARM200F TEM with CEOS probe and image aberration correctors at 80 kV. The annular dark-field (ADF) images were recorded at a probe current of ∼23 pA and a convergence semi-angle of ∼25 mrad using a JEOL annular field detector with an inner and outer collection semi-angle of 45 and 180 mrad, respectively [36,37]. The scanning rate was typically 20 µs per pixel and each image consists of 2048 × 2048 pixels [36].

Determining the atom positions
Image analysis of the STEM data was used to analyse the atomic configuration in the alloys. Thresholding, based on integrated intensity within a region centred on each atom, was applied to determine the atomic identities and this was then verified visually. The number of W-W pairs at nth-nearest-neighbour distances was counted for n = 1-18, based on the positions of the W atoms in the alloy. Next, the total number of transition metal-transition metal pairs at each distance was counted. The number of W-W pairs is divided by the product of the total number of transition metal-transition metal pairs and the W atom composition x.

Electronic structure measurements
ARPES data were acquired at the Spectromicroscopy beamline [38] of the Elettra synchrotron, using 27 eV photon illumination and a sample temperature of around 100 K. Spectra were acquired from crystals which were cleaved in ultrahigh vacuum (UHV) immediately before analysis. Linearly polarised light was focused to a submicrometre spot and uniform regions of the crystals were identified by scanning photoemission microscopy to ensure that ARPES data were collected from representative regions within single domains. Three-dimensional data sets of the photoemitted intensity as a function of kinetic energy and emission angles were acquired around the line in reciprocal space along the high symmetry direction from Γ to K. From this, 2D energy-momentum slices were extracted showing the dispersion from Γ to K. The dispersion of the upper valence band was determined by fitting the corresponding peaks in the energy distribution curves (EDC) with Lorentz functions. The effective mass was found by fitting a parabola to the resultant band dispersions around K.

Low-temperature polarisation-resolved PL spectroscopy and reflectance measurements
The PL was excited using a continuous wave (CW) dye laser with a wavelength of 625 nm via a 50 × 0.55 NA microscope objective producing a 2 µm laser spot. The power incident on the sample was kept below 100 µW. The measurements were performed in backscattering geometry with the sample in vacuum at 4 K. In order to ensure good circular polarisation purity, a quarter wave plate (Thorlabs AQWP05M-600) was placed just before the objective so that in the rest of the optical system the input and output polarisations corresponded to vertical and horizontal linear polarisations. The output polarisation was determined using a linear polariser (Thorlabs LPVIS100) and a half wave plate (AHWP05M-600) used to ensure a constant input polarisation for the spectrometer. The spectrometer used was a Princeton Instruments TriVista 555 with a liquid nitrogen cooled CCD. Measurements were repeated on multiple areas of the sample. The reflectivity measurements were taken on the same sample areas, also at 4 K with the same microscope objective. The light source used was a Fianium supercontinuum laser producing a spot size of 3 µm. A Princeton Instruments FERGIE spectrometer was used to capture the spectrum.

Monte Carlo simulations
We model the equilibrium distribution of the cations W and Mo using an Ising-like model of nearest-neighbour interactions between pairs of cations on a hexagonal lattice. The model comprises a single layer with 120 × 120 in-plane unit cells, a size corresponding to several times the largest apparent feature size in the final distributions. This was initialised to a random distribution constructed to achieve the target concentration x. The growth process was simulated by Monte Carlo sampling on proposed site swaps to achieve an equilibrium of the system after 20 000 cycles at a temperature of 800 K. To parameterise the Monte Carlo simulations, DFT calculations were performed on a selection of specifically-chosen models, as shown in suppmat figure S3(a). These models contained between 1 and 6 substitutional W atoms, arranged in small clusters in a 4 × 4 supercell. A binding energy was calculated for each cluster containing m atoms of W, relative to m well-separated substitutional W atoms, according to E binding (m) = E cluster (m) − mE 1W + (m − 1) E 0W . Then these binding energies were used to find the average binding energy per bond between adjacent W atoms, for each cluster [39]. To validate this, DFT calculations were performed for a selection of atomic configurations extracted from those observed in the experimental STEM images (see suppmat figure S3(c) and (d)). For illustration and comparison, Monte Carlo simulations were also performed with interaction energies of −50, 0 and 50 meV.

Plane-wave DFT
For calculations of the optimised lattice constants of the pure materials, MoS 2 and WS 2 , the plane-wave DFT package CASTEP was used [40]. A cut-off energy of 424.5 eV and a 10 × 10 × 1 k-point grid were used, and the generalised gradient approximation of Perdew, Burke and Ernzerhof (PBE) was employed [41]. Since the resulting lattice constants for MoS 2 and WS 2 are very similar (within 0.3%), the same value, 3.18 Å, was subsequently used for both materials, thereby neglecting the very small lattice mismatch. CASTEP calculations were also used to provide the interaction energies for pairs of W atoms for the Monte Carlo parameterization, with commensurate adjustments of the k-point grid to 3 × 3 × 1 and 1 × 1 × 1 as the supercell size increased to 4 × 4 × 1 and 12 × 12 × 1 respectively.

LS-DFT
The ONETEP LS-DFT package [42] was applied to simulate large supercells with random alloy configurations. ONETEP uses a representation of the single-particle density matrix in terms of a minimal number of support functions denoted non-orthogonal generalised Wannier functions (NGWFs). The NGWFs are themselves represented in terms of a grid of psinc functions [43] (a systematic basis equivalent to plane waves), which here use a cut-off energy of 1200 eV, and an NGWF cut-off radius of 10a 0 . The PBE functional was used, with projector-augmented wave (PAW) [44,45] potentials from the JTH library [46]. The electronic properties of different random configurations of atomic structures were calculated within an 8 × 8 × 1 supercell model. Mo 1−x W x S 2 monolayer models were created for various concentrations, x, with cation distributions determined according to the procedures described below. A separate geometry optimisation was performed for each structure, with force tolerance of 0.1 eV Å −1 . Since the band structures obtained from the disordered alloy supercell calculations have a heavily-folded band structure that is not directly comparable to the ARPES results, the effective band structure was calculated in the primitive cell using the spectral function unfolding method [14,47], the adaptation of which to the NGWF representation is described in [48].

Virtual crystal approximation (VCA) calculations
Calculations were performed within the virtual crystal approximation (VCA) [49] with the Quantum Espresso code [50,51]. PAW datasets representing the W and Mo potentials mixed in the appropriate ratio x were constructed using the 'virtual' tool within the Quantum Espresso package. SOC was included using the treatment described by dal Corso [52][53][54].

Results and discussion
Mo 1−x W x S 2 alloy crystals were grown by chemical vapour transport, as described in section 2, across the range from x = 0 to 1. Their compositions, determined by a combination of XPS and EDX analysis (see section 2 and suppmat S2 for further details), were found to be consistent within 3% between crystals grown in the same batch, and across individual crystals when averaged at the micrometre scale. When exfoliated to monolayers, they exhibit room temperature PL with peak emission frequency varying nonlinearly from 1.85 eV at x = 0 to 1.98 eV at x = 1, consistent with prior literature reports [14] (suppmat S4).

Atomic structure
To determine whether the alloys are atomically ordered, atomic resolution ADF-STEM images of monolayer Mo 1−x W x S 2 flakes were acquired. Example images are given in figure 1(a). In the magnified view, figure 1(b), the transition metal atoms can be resolved but the chalcogen atoms give only weak contrast, confirmed by the corresponding simulated STEM image of the same area in figure 1(c). Due to the large difference between their atomic numbers, there is also a clear contrast difference between Mo (Z = 42) and W (Z = 74), again verified by comparison between the line profiles of the experimental and simulated images, as shown in figure 1(d).
For comparison, we performed Monte Carlo simulations made with varying interaction energies, J, and compositions, x. Details of the simulation methodology are given in the methods and suppmat S3. To predict J for these alloys, binding energies for a range of small cluster configurations were calculated (see table S2). From these results it can be seen that the binding energy per bond is not strongly sensitive to cluster size for small clusters (m = 2-6), and indeed clusters of this size corresponds to those observed in Mo 1−x W x S 2 for low through intermediate values of x. Therefore we can use the binding energy associated with the m = 2 cluster as a representative value, hence we apply J = 7.6 meV for the generation of configurations for comparison to experiment. This was further verified by comparison of the Monte Carlo model to DFT calculations of a range of larger cells (suppmat S3). To elucidate the effect of this binding energy, we also test the impact of stronger interactions (J = ±50 meV). Figure 2(a) shows representative atomic arrangements   from Monte Carlo simulations, comparing atomic arrangements for the predicted value of J = 7.6 meV to those with stronger interactions (J = ± 50 meV) at a simulation temperature of 800 K. At negative J, clusters of W tend to form, while at positive J, isolated W atoms are energetically preferred.
For J = 0, when there is no energy difference between the W and Mo, they are randomly arranged; qualitatively, Monte Carlo simulations with J = 7.6 meV look similar to those with J = 0. Although close visual inspection of the images can suggest some ordering, quantitative statistical analysis is required to rigorously test whether the Mo/W atoms are randomly distributed within the sheets. suppmat S3 shows that in the Monte Carlo model, strong ordering disappears beyond 50 K and is completely absent by temperatures comparable to the growth conditions. For a given composition, ADF-STEM images with areas of 120 × 120 atoms were analysed, fitting the intensity at each atom coordinate and thresholding to determine the arrangement of Mo and W atoms, as described in section 2. With the atom coordinates identified, statistical analysis can be used to investigate whether the arrangements are ordered. Following Cowley [55], we define the short range order parameter: where n i is the number of W atoms out of the c i atoms in the ith coordination shell surrounding each W atom (as defined schematically in figure 3), and x is the proportion of W atoms (to (1 − x) Mo atoms) in the alloy. Figure 2(e) shows the order parameters for Mo 0.78 W 0.22 S 2 determined from the ADF-STEM images alongside those calculated from analysis of the Monte Carlo simulations. Order parameters for the other compositions followed similar trends. For the simulated data, if J > 0 then α 1 < 0, showing a lower probability of finding a W atom neighbouring another W atom. The magnitude of α i decreases rapidly with increasing i, indicating a finite range to the order. If J < 0 then α 1 > 0 (a tendency to cluster), and α i decreases monotonically to 0 with increasing i. The order parameters determined from the experimental data agree, within uncertainties, with data from both J = 0 and J = 7.6 meV simulations, consistent with a random distribution of W atoms. We conclude that under the conditions used for the growth of these crystals, the resultant atomic distributions are indistinguishable from a random distribution and are consistent with those expected at thermodynamic equilibrium at the elevated growth temperature of this system.

Electronic structure
To correlate the electronic structure to the atomic structure, the valence band dispersions of the Mo 1−x W x S 2 alloy crystals were measured by ARPES. Since the measurements are surface sensitive, the as-grown bulk crystals were mounted flat onto a sample plate and cleaved in the UHV chamber immediately prior to analysis, to generate the clean surface required. Cleaving naturally results in exposure of the (002) planes to the beam, so that ARPES probes the electronic structure in the plane of the 2D layer.
Typical spectra are given in figure 4, showing the dispersion from Γ to K for three alloy compositions: Mo 0.63 W 0.37 S 2 , Mo 0.5 W 0.5 S 2 and Mo 0.29 W 0.71 S 2 . We concentrate our analysis on the structure around the zone corners at K where, for the monolayers, the VBM is found. Around K, the spin-split upper valence band has in-plane orbital composition from the metal d x2−y2 + d xy and, to a lesser extent, the chalcogen p x+y [4]. This results in a weak k z dispersion such that the key features of the spectra in that region are quantitatively similar for bulk and monolayer materials [56,57]. From this dispersion, band parameters can be determined such as the spin-orbit splitting of the bands at K, the effective mass, and the band width (see the schematic in figure 5).
The qualitative trends are apparent by visual inspection of the ARPES spectra: as the W content increases, the spin-orbit splitting between the upper valence bands near K (SOC VB ) increases; the band width (D vb ) increases and the effective mass decreases. To quantitatively evaluate the evolution of the band structure with composition, measurements were performed on all the compositions (x = 0.22, 0.37, 0.49, 0.71, 0.87 and 1) and band parameters extracted quantitatively by fitting the ARPES spectra, as described in section 2. The data are presented in figure 6, and compared to first principles calculations as described below.
Ab initio calculations of the electronic structure across the alloy compositions were made using two approaches: the VCA [49], where fractional atom compositions are considered within the primitive cell of pure materials; and the more computationally intensive approach of calculating the electronic structure of large supercells with varying ratios of Mo:W atoms. The details of each approach are described in Methods and supplementary materials. In brief, the Quantum Espresso package [50,51] was used for the VCA calculations, and LS-DFT calculations were performed on 8 × 8 × 1 supercells (192 atoms) by the ONETEP package [58]. The LS-DFT calculations were unfolded into the primitive cell of the pure monolayers [48]; an example spectral function, for Mo 0.5 W 0.5 S 2 , is shown in figure 5, and spectral functions for the other concentrations can be seen in suppmat figure S7. Band structure parameters were determined from fitting the spectral functions. Figure 6 compares the band structure parameters measured from the experimental data to those calculated by LS-DFT and VCA simulations, as a function of the W content in the alloy. As the ARPES spectra probe the valence band, direct experimental measurements of the conduction band edge and band gap (E g ) are not possible from this data. Instead, the predicted E g are compared to the PL emission peak energy in figure 6(a). DFT calculations normally underestimate E g [6], so we concentrate on the trends with composition. As previously reported [14], the PL peak energy changes non-linearly with composition, with the bowing parameter here determined to be 0.17 ± 0.01 eV. This non-linear dependence of E g on x is expected due to the change in orbital characteristic of the states at the conduction band edge from MoS 2 to WS 2 [14]. LS-DFT predicts qualitatively similar changes in E g , with a bowing parameter of 0.12 ± 0.01 eV, but with a smaller change in E g (E g = 1.63 eV for MoS 2 and 1.62 eV for WS 2 ) than observed in the change in PL peak emission energy (1.85 eV for MoS 2 increasing to 1.98 eV for WS 2 ). The magnitude of E g determined from VCA is similar to the LS-DFT calculations, but the bowing parameter becomes negative (−0.02 ± 0.01 eV), showing that the VCA predictions deviate from the experimental observations. Caution should be taken when comparing the PL emission peak energies to the predicted band gaps (E g ) which are determined from the energy difference between the VBM and the CBM. The large exciton binding energy (several hundred meV) means that the PL peak emission energy is not a good measurement of the single-particle electronic band gap [2]. Interpretation of the PL is further complicated by the change in order of the spin-polarised conduction and valence bands, as shown schematically in figure 5(b). The PL emission is from the A exciton (bright exciton) in the K valley, corresponding to an electron transition between states of the same spin. For x = 0 (MoS 2 ), this corresponds to the VBM and CBM states. But at x = 1 (WS 2 ), this has changed and the bright A exciton emission comes from the second-lowest lying conduction band which is greater in energy than the CBM by the spin-orbit splitting at the conduction band edge (around 30 meV for WS 2 ).
Prior direct measurements of E g from ARPES of electrostatically doped monolayers found E g = 2.07 ± 0.05 eV for MoS 2 and E g = 2.03 ± 0.05 eV for WS 2 [59], only a small change in band gap as also seen in the LS-DFT results. The difference between the LS-DFT and PL data then suggests a consistent decrease in exciton binding energy as the W content increases.
ARPES provides a direct and visual measurement of the valence band dispersion of the alloys, allowing easier comparison to calculated spectra than is possible from optical spectroscopy data. Measured values of SOC VB from the Mo 1−x W x S 2 single crystals (blue data points) are compared in figure 6(b) to those calculated for the monolayer alloys by LS-DFT (black) and VCA (red). The measured value for WS 2 (SOC VB = 458 ± 10 meV) is consistent with prior reports for the bulk crystals (450 ± 30 meV [39]), and for a monolayer flake of WS 2 (450 ± 30 meV [59]). For the experimental results, SOC VB increases linearly as a function of the W content, within the error bars. The LS-DFT predictions show a similar linear trend, with a systematic under-estimate of the experimental results by around 50 meV. The VCA calculations of SOC VB are similar in magnitude to the experimental measurements but show a non-linear dependence on composition. An LS-DFT simulation was also made using the atom positions determined directly from a STEM image: the band parameters for this simulation are shown as the purple data points in the graphs in figure 6 and are consistent with the other LS-DFT simulations. The band width (D vb ), figure 6(c), follows similar trends to SOC VB : the agreement between the experimental data and the LS-DFT data is excellent, other than a systematic offset (under-estimate) between the two; and the VCA data shows a bowing that is not present in either the LS-DFT or experimental data. Figure 6(d) compares the evolution of the valence band effective mass (around K and in the Γ to K direction) calculated from ARPES spectra, the LS-DFT and VCA calculations. It is difficult to accurately determine the effective mass from the ARPES spectra, as a result the uncertainties are large and there are no clear trends. From the DFT calculations it is clear that the effective mass decreases monotonically with x, consistent with the increasing curvature of the bands and increasing band widths. This decrease in effective mass is also consistent with a decrease in exciton binding energy with increasing W content.
The ARPES spectra do not give information about the conduction band as they only probe the occupied states. However, the spin-orbit splitting at the conduction band edge (SOC CB ) can be extracted from the DFT spectral functions and this is plotted as a function of composition in figure 6(e). It is an order of magnitude lower than SOC VB . At small x, closer to pure MoS 2 , the order of the conduction bands changes, see figure 5(b). Hence SOC CB initially decreases to zero and then increases with x. Although the LS-DFT and VCA calculations do not quantitatively agree over the composition at which SOC CB = 0, it is clear that at intermediate compositions SOC CB is small. At x = 0 and x = 1, the small differences between the results are due to the different approaches to SOC taken by the two simulation packages: ONETEP uses a perturbative approach to SOC, whereas Espresso uses a full SOC treatment integrated into the self-consistency procedure.
Calculations were performed on three different random realisations of the 8 × 8 supercells for each x, allowing a rough estimate of the disorder potential to be made, on a similar lengthscale to the exciton radius, based on the variation of the CBM position between these random realisations. From these we observe that the disorder potential is larger than SOC CB at least for 0.125 ⩽ x ⩽ 0.5, becoming smaller again as SOC CB increases (see suppmat S5). This casts doubt over the usual picture of separated states split by SOC at the conduction band edge and merits further investigation of the spin-valley locking in these alloys [60].

Spin-valley locking
To test whether spin-valley locking is retained in the alloys, optical spectroscopy measurements were performed at 4 K on an exfoliated heterostructure consisting of a monolayer flake of Mo 0.5 W 0.5 S 2 encapsulated in hBN. An optical image, and schematic, of the sample is shown in figure 7(a). Polarisation-resolved PL spectra are presented in figure 7(b) and the differential reflectivity spectrum is presented in figure 7(d). Figure 7(d) also presents a fit to the differential reflectivity spectra assuming a single Lorentzian oscillator with a centre energy of 1.954 ± 0.001 eV and width of 41 ± 1 meV. This energy falls between the reported energies for free A1s excitons in encapsulated monolayers of MoS 2 and WS 2 reported in the literature which are ∼1.94 [61] and ∼2.09 eV [62] respectively. In non-alloy monolayers defect bound excitons are not observed in reflectivity spectra even when they are observed in PL spectra [63,64]. This justifies the assignment of the reflectivity feature in the alloy sample with free excitons. The width of the free exciton peak determined from the reflectivity peak is significantly broader than the width of the same peak in non-alloyed samples, ∼5 meV. We attribute this to alloy fluctuations on length scales longer than the Bohr radius of the excitons. The PL spectra show a broad feature from ∼1.85 eV to ∼1.97 eV. Equivalent PL spectra in non-alloy monolayers show a series of separated peaks including free excitons at higher energies, trions/polarons at intermediate energies and defect bound excitons at lower energies [5,[64][65][66]. It is likely that such features would be broadened in alloy samples and thus we attribute the feature observed in the alloy sample to a range of different excitonic states whose character changes from bound to free exciton with increasing energy. The PL spectra show clear circular dichroism across the full spectral range of the feature: excitation with righthand circularly polarised light (σ+) leads to emission of mainly σ+ light, whilst σ− excitation leads to σ− emission. Excitation with vertical linear polarisation (corresponding to (σ+ + σ−)/√2) results in emission with an equal amount of σ+ and σ− circular components. The PL spectra are scaled so that the total intensity for both circular polarisations at the peak of the feature is one.
The valley polarisation, η = PL(σ + )−PL(σ − ) PL(σ + )−PL(σ − ) where PL(σ±) is the intensity of the corresponding circularly polarised component of the PL emission [9], is plotted in figure 7(c). We observe a degree of circular polarisation of above 0.4 across the energy range 1.87-1.97 eV and particularly in the vicinity of the free excitons at 1.954 eV. These results are fully consistent with the expectations from the optical selection rules for the spin-polarised monolayer materials and with previous observations of spin-valley locking in pure monolayer MoWSeS flakes, [5][6][7][8] suggesting that spin-valley locking is retained for the monolayer alloys despite the atomic-scale heterogeneity in the alloys and the small SOC CB . This indicates that locally, over the length scale of the exciton, the spin texture is retained.

Conclusions
No ordering or segregation was found in CVT grown Mo 1−x W x S 2 alloys, consistent with Monte Carlo simulations parameterised by ab initio calculations. ARPES measurements showed that the SOC at the valence band edge increases linearly with increasing x (increasing W content). LS-DFT simulations using large unit cells gave predictions consistent with the ARPES measurements, but DFT simulations using fractional atom concentrations in a primitive unit cell were found to be inconsistent with the experimental results. This conforms to expectations of the known severe limitations of the predictive power of VCA, which becomes even more apparent for properties that are nonlinear with atomic number Z, such as SOC. This is especially apparent for materials which retain some localised character in the valence band, and strongly suggests instead the need for large-scale disordered models to study alloy properties. From the LS-DFT models, for intermediate alloy compositions the SOC is predicted to be smaller than the spatially localised changes in the conduction band edge due to the atomic-scale heterogeneity of the alloy. Despite this, polarisation-resolved PL measurements of a monolayer Mo 0.5 W 0.5 S 2 flake show circular dichroism consistent with spin-valley locking. The ability to rationally and continuously tune the SOC, effective mass and band gap suggests a promising future for MX 2 alloys in valleytronic and spintronic applications.