Atomic and electronic structure of two-dimensional Mo_(1−x)W_xS₂ alloys

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⁻⁶ mbar and sealed. The mixture was heated to 1000 °C for 3 d to form Mo_1−x W_x S₂ 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₂ (10 mg cm⁻³ of the ampoule volume). To keep the I₂ stable, the ampoule was evacuated to 10⁻⁶ 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₂ 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.

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.

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₂ 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].

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.

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 $\boldsymbol{\overline \Gamma}$ to $\boldsymbol{\overline K}$. From this, 2D energy-momentum slices were extracted showing the dispersion from $\boldsymbol{\overline \Gamma}$ to $\boldsymbol{\overline 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 $\boldsymbol{\overline K}$.

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.

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_{{\text{binding}}}}\left( m \right) =$ ${E_{{\text{cluster}}}}\left( m \right) - m{E_{1{\text{W}}}} + \left( {m - 1} \right){E_{0{\text{W}}}}$. 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.

For calculations of the optimised lattice constants of the pure materials, MoS₂ and WS₂, 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₂ and WS₂ 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.

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₀. 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₂ 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 Å⁻¹. 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].

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–54].

To determine whether the alloys are atomically ordered, atomic resolution ADF-STEM images of monolayer Mo_1−x W_x S₂ 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).

Zoom In Zoom Out Reset image size

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₂ 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.

Zoom In Zoom Out Reset image size

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: ${\alpha _i} = 1 - {n_i}/\left( {x{\text{ }}{c_i}} \right)$, 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 $\left( {1 - x} \right)$ Mo atoms) in the alloy.

Zoom In Zoom Out Reset image size

Figure 2(e) shows the order parameters for Mo_0.78W_0.22S₂ 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 ${\alpha _1} < 0$, showing a lower probability of finding a W atom neighbouring another W atom. The magnitude of ${\alpha _i}$ decreases rapidly with increasing i, indicating a finite range to the order. If J < 0 then ${\alpha _1} > 0$ (a tendency to cluster), and ${\alpha _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.

To correlate the electronic structure to the atomic structure, the valence band dispersions of the Mo_1−x W_x S₂ 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 $\boldsymbol{\overline \Gamma}$ to $\boldsymbol{\overline K}$ for three alloy compositions: Mo_0.63W_0.37S₂, Mo_0.5W_0.5S₂ and Mo_0.29W_0.71S₂. 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_{x 2−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).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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.5W_0.5S₂, 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₂ to WS₂ [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₂ and 1.62 eV for WS₂) than observed in the change in PL peak emission energy (1.85 eV for MoS₂ increasing to 1.98 eV for WS₂). 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₂), this corresponds to the VBM and CBM states. But at x = 1 (WS₂), 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₂).

Prior direct measurements of E_g from ARPES of electrostatically doped monolayers found E_g = 2.07 ± 0.05 eV for MoS₂ and E_g = 2.03 ± 0.05 eV for WS₂ [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₂ 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₂ (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₂ (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₂, 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].

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.5W_0.5S₂ 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₂ and WS₂ 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–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, $\eta = \frac{{{\text{PL}}\left( {{\sigma ^ + }} \right) - {\text{PL}}\left( {{\sigma ^ - }} \right)}}{{{\text{PL}}\left( {{\sigma ^ + }} \right) - {\text{PL}}\left( {{\sigma ^ - }} \right)}}$ 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–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.

Zoom In Zoom Out Reset image size