First-Principles Calculation of Optoelectronic Properties in 2D Materials: The Polytypic WS2 Case.

The phenomenon of polytypism, namely unconventional crystal phases displaying a mixture of stacking sequences, represents a powerful handle to design and engineer novel physical properties in two-dimensional (2D) materials. In this work, we characterize from first-principles the optoelectronic properties associated with the 2H/3R polytypism occurring in WS2 nanomaterials by means of density functional theory (DFT) calculations. We evaluate the band gap, optical response, and energy-loss function associated with 2H/3R WS2 nanomaterials and compare our predictions with experimental measurements of electron energy-loss spectroscopy (EELS) carried out in nanostructures exhibiting the same polytypism. Our results provide further input to the ongoing efforts toward the integration of polytypic 2D materials into functional devices.

Introduction

Two-dimensional (2D) materials of the transition metal dichalcogenide (TMD) family have attracted ample interest due to the wide range of tunability of their electronic and optical properties.[1−6] This flexibility in tailoring the physical properties of TMD materials can be traced back to their marked sensitivity with respect to the dimensionality,[7,8] specific edge configurations,[9−11] and stacking sequences.[12−15] The most common stacking sequences (polytypes) present in TMD materials are those based on the octahedral coordination (1T) and on the trigonal prismatic coordination (2H and 3R).[16−19] Furthermore, it has been reported that mixed 2H and 3R polytypes can also occur within the same TMD-based nanomaterial. This phenomenon, known as polytypism, leads to unconventional crystal phases displaying a mixture of stacking sequences and has been reported among several other TMD materials, including MoSe2, MoS2, and WS2.[20−22]

The different stacking sequences or polytypes exhibited by TMD nanostructures have associated specific variations in the resulting electronic and optical properties. For instance, 1T-MoS2 nanosheets are found to be metallic,[18] while their 2H and 3R counterparts display instead a semiconductor behavior. While extensive investigations of the optoeletronic properties of both the 2H and 3R polytypes in TMD materials have been carried out, much less is known about the implications of the 2H/3R mixed crystal phase. In this respect, achieving a deeper understanding of the optoelectronic properties associated with this 2H/3R polytypism is a key component of the ongoing efforts toward integrating polytypism-based TMD materials into functional devices.

Motivated by this, in this work we carry out first-principle calculations of the optoelectronic properties associated with the 2H/3R polytypism occurring in WS2 nanomaterials by means of density functional theory (DFT)[23,24] as implemented in the WIEN2k framework. We ascertain the semiconducting nature of 2H/3R polytypism and evaluate the corresponding band gap, the value of which is found to be in agreement with a recent experimental determination.[25] To obtain accurate estimates of the band gap and band structure associated with the different polytypes of WS2, we use the GW approximation as implemented in the GAP2 code by Jiang et al.[26,27] and that takes as input the output of the DFT calculations from WIEN2k. By combining DFT with GW calculations, we are able to evaluate the density of states (DOS), the band structure, the joint density of states (JDOS), and the energy-loss function (ELF) corresponding to 2H/3R WS2 nanomaterials. Furthermore, the calculations of the ELF can also be compared with the experimental measurements from electron energy-loss spectroscopy (EELS) carried out on WS2 nanostructures displaying the same 2H/3R polytypism.

Our analysis hence makes possible disentangling the dependence of the electronic and optical behavior of the 2H/3R mixed WS2 nanostructures on the underlying crystal structure. In turn, our results can be used to motivate further investigations of the growth mechanism of polytypes in 2D materials, and ultimately explore the possibilities in terms of achieving heterostructures based in a single TMD material.

Computational Techniques

The electronic properties of the 2H, 3R, and mixed 2H/3R polytypes of WS2 are investigated using both the linearized augmented plane wave (LAPW) and local orbitals (lo) methods as implemented in the WIEN2k Density Functional Theory package.[28] For the density of states (DOS) calculation, the van der Waals (vdW) interaction characteristic of WS2 polytypes is taken into account by implementing the nonlocal vdW functional model.[29,30] The nonlocal vdW interactions use optB88[31] for the exchange term, the local density approximation[32] (LDA) for the correlation term, and the DRSLL kernel for the nonlocal term.[33] For the nonlocal vdW integration, the cutoff density rc is set to 0.3 bohr–3, while the plane wave expansion cutoff Gmax is set to 20 bohr–1. No spin polarization is considered.

The equilibrium lattice parameters for each of the three polytypes considered are found by volume and force optimization of the different unit cells, such that the force on each atom is less than 1.0 mRy/bohr. The total energy convergence criteria is set to be 0.1 mRy between self-consistent field (SCF) cycles, while the charge convergence criteria is set to 10–3e, with e being the elementary unit charge. The core and valence electron states were separated by an energy gap of −6.0 Ry. Furthermore, the calculations used a Rkmax value of 7.0, where R is the radius of the smallest Muffin Tin sphere and kmax is the largest k-vector.

Geometry Optimization

A geometrical optimization is carried out in order to find the equilibrium lattice parameter corresponding to the three 2H, 3R, and mixed 2H/3R WS2 polytypes. For each polytype, a specific amount of k-point sampling is employed until convergence is achieved with respect to the k-points, as indicated in Table . The k-points sampling of the first Brillouin zone for the lattice parameter calculations and all subsequent calculations using the DFT formalism are carried out using the tetrahedron method of Blöchl et al.[34]

Table 1

Values of the k-Mesh Used for the DFT Calculations and Crystal Structuresa

polytype	k-mesh	calculated lattice parameters
2H	16 × 16 × 3	a = b = 3.194 Å, c = 12.458 Å
3R	14 × 14 × 14	a = b = 3.199 Å, c = 18.733 Å
2H/3R	24 × 24 × 3	a = b = 3.205 Å, c = 19.057 Å

We also indicate, for each WS2 polytype, the calculated lattice parameters a, b, and c. For the 2H and 3R polytypes these values are consistent with the experimental values reported in the literature.[35].

The DFT calculations presented in this work are based on the structural atomic models displayed in Figure . Figure schematics a and b display the atomic model of the 2H and 3R polytypes, respectively. The 2H polytype exhibits the characteristic hexagonal stacking order (AA′), where the adjacent layers are rotated 180° and stacked directly upon one another. The 3R polytype is characterized by the rhombohedral stacking order (BA), in which the adjacent layers are slightly displaced from each other without any rotation. Figure c shows the atomic model of the mixed 2H/3R WS2 polytypes, which is characterized by a layer stacking order of the type BAA′. This stacking order arises from the mixture of the 2H (AA′) and 3R (BA) polytypes.[36]

Figure 1

Schematic for the different polytypes of WS2. In the top panels, it is shown how all three polytypes exhibit a hexagonal structure when viewed from the [0001] direction; (a) the 2H crystal phase exhibits a honeycomb lattice, while the (b) 3R and (c) 2H/3R polytypes both have an atom in the middle of their honeycomb structures. The stacking sequences of the 2H, 3R, and their mixed 2H/3R polytypes can be better assessed when viewed from a lateral viewpoint with respect to the layers, as illustrated in the bottom panels.

The equilibrium lattice parameters obtained (see Figure S1 in the Supporting Information) for the three polytypes are shown in Table . The lattice parameters found for 2H and 3R are consistent with previously reported values in the literature.[35] For the mixed 2H/3R polytype, the lattice parameters are found to be a = b = 3.205 Å, and c = 19.057 Å. Furthermore, using the symmetry package available in WIEN2k, one can identify that the mixed 2H/3R phases belongs to the P3m1 space group.

GW Approximation

After geometry optimization, the equilibrium lattice parameters obtained are used for the calculations using the GW approximation as implemented in the GAP2 code. There, the Green’s function G and the screened Coulomb interaction W are calculated in the RPA framework. For the 2H and mixed 2H/3R crystal phases, we use a 6 × 6 × 1 k-mesh sampling of the Brillouin zone, while a 6 × 6 × 6 k-mesh sampling is used instead for the 3R crystal phase case. Fourier interpolation is used to interpolate a fine k-mesh from the sparse-mesh originally used in the GW calculations to obtain quasiparticle energies.[37] Spin–orbit coupling effects are taken into account in a perturbative one-shot manner. Table displays the values of the k-mesh used for the GW calculation and for the calculation of the DOS, band structure, and ELF. These k-meshes are kept the same for the calculations with and without spin-orbit coupling.

Table 2

Values of the k-Mesh Used for the GW Calculation and for the Calculation of the DOS, Band Structure, and ELFa

crystal phase	GW calculation	DOS, band structure, ELF
2H/3R	6 × 6 × 1	23 × 23 × 3
2H	6 × 6 × 1	20 × 20 × 4
3R	4 × 4 × 4	16 × 16 × 16

These k-meshes are kept the same for the calculations with and without spin-orbit coupling.

Results and Discussion

Density of States and Band Structure

Figure panels a and b display the total density of states (TDOS), the projected density of states (PDOS), and the electronic band structure corresponding to the mixed 2H/3R polytype, together with their counterparts for the 2H and 3R polytypes, Figure panels c and d and panels e and f, respectively. The TDOS corresponding to the three polytypes are similar with and without SOC effects. By comparing the TDOS with the PDOS of the mixed 2H/3R, 2H, and 3R, one can observe how the W d-orbitals are the primarily allowed occupied and unoccupied states near the Fermi energy. This finding is consistent with related DFT studies based on the MoS2 bulk crystal structure.[2,38]

Figure 2

Left panels: the calculated density of states associated with the (a) 2H/3R, (c) 2H, and (e) 3R polytypes with and without spin–orbit coupling (SOC) taken into account. Right panels: the resulting band structures of the (b) 2H/3R, (d) 2H, and (f) 3R polytypes evaluated with and without SOC.

For the band structure calculations, the k-paths in the primitive Brillouin zones are chosen in a way such that the complete irreducible set of symmetry lines is obtained. In the case of 2H and mixed 2H/3R polytypes, identical k-paths (Γ → M → K → Γ → A → L →H → A) in the Brillouin zone are chosen since their primitive Brillouin zone have the same symmetry points, see Supporting Information Figure S2. Notice that for the 3R polytype we follow the prescriptions of refs (39−44) for the selection of the k-path in the primitive Brillouin zone, see also Supporting Information Figure S3.

The DOS and band structure of the 2H and 2H/3R polytypes (with and without SOC effects) reveal that the conduction band (CB) minima is found to be at the same point in the Brillouin zone, located between the K and Γ points. The valence band (VB) maximum for the 2H polytype happens at the Γ point, while for the mixed 2H/3R polytype it lies at the A point. In the band structure of the mixed 2H/3R polytypes, near the Fermi energy an extra band appears (Figure b) as compared to the two bands in the 2H case (Figure d). This extra band near the Fermi energy in the mixed 2H/3R polytype is the origin of the shift of the VB maximum from the Γ point to the A point. As a consequence of the extra band near the Fermi energy, an additional band can be found near the K point, leading to additional bands appearing when taking spin–orbit coupling into account (see Supporting Information Figure S4). Therefore, the mixed 2H/3R polytype in bulk form exhibits an indirect band gap like in the 2H bulk form. Concerning the band gap energy of the mixed 2H/3R polytype, our calculations indicate a value of 1.40 eV (1.48 eV) with (without) SOC effects. These values are about 5% higher than for the 2H polytype, where the calculated band gap energy is instead 1.34 eV with SOC and 1.39 eV without SOC.

In the case of the 3R crystal structure, the VB maximum lies at the Z point, while the CB minimum lies between the Γ and X points in the Brillouin zone. Furthermore, the band near the Fermi energy is split into two bands at the B, X, and Q points in the Brillouin zone when including spin–orbit coupling. The CB minimum lies between the Γ and X points and is pushed down in energy when including spin–orbit coupling, while the maximum of the VB is not modified significantly by the inclusion of the spin–orbit coupling. Therefore, the resulting band gap for the 3R polytype turns out to be 1.54 eV without SOC and 1.39 eV with SOC. Hence one observes that the SOC effects have a more significant effect for the 3R polytype as compared to the 2H and 2H/3R cases, and indeed in Figure e one can see how the band gap value markedly decreases once SOC is taken into account. This effect might originate from the one-shot postinclusion of spin–orbit coupling in the GW calculations.

The indirect band gap value of the mixed 2H/3R polytype is consistent with the experimental one reported in refs (22 and 25). Furthermore, such a value lies in between those of the 2H and 3R polytypes in the cases with and without spin–orbit coupling.

Energy-Loss Function

An improved understanding of the optical response of the unconventional mixed 2H/3R polytype can also be achieved by means of the calculation of the energy-loss function (ELF). With this motivation, we now compare the ELF calculated at the DFT and at the GW levels, with and without taking into account the spin–orbit coupling interaction effects.

Figure panels a and b display the ELF calculated along the out-of-plane direction (c-axes), a configuration which matches that of the underlying experimental EELS measurements. In the low-loss region below 10 eV, the DFT calculations exhibit distinctive features at 7 and 10 eV, while the GW calculations exhibit a well-defined peak at around 8 eV. For the region with energy losses between 10 and 30 eV, the main feature for both the DFT with (without) SOC and the GW with (without) SOC calculations is a peak located at 20.9 eV (21.4 eV) and 21.1 eV (21.7 eV), which can be identified with the bulk plasmon of WS2. Interestingly, the bulk plasmon locations differ by no more than 0.25 eV for both the DFT and GW calculations with and without SOC effects. This result can be understood since all the band structures contribute to the bulk plasmon feature, therefore the effect of the GW calculation is less significant when compared to that of the low-loss region where the features depend more sensitively on bands near the Fermi energy.

Figure 3

(a,b) Out-of-plane energy-loss function of the mixed 2H/3R polytype calculated on the DFT and GW level with and without taking spin–orbit coupling into account. (c,d) Experimental EELS measurements acquired on a WS2 nanostructure characterized by the same 2H/3R polytypism. Note that the zero-loss peak (ZLP) has not been subtracted from the EELS data in the bottom left panel.

Figure panels c and d display then experimental EELS measurements acquired on a WS2 nanostructure also characterized by the same 2H/3R polytypism[22] in the same energy-loss range that the corresponding theoretical calculations determined. Note that the zero-loss peak (ZLP) has not been subtracted from the EELS data in the bottom left panel. In these experimental measurements, one observes in the low-loss region two peaks located at 3.5 and 8 eV. The peak at 8 eV, which can be associated to the interlayer coupling, is consistent with the features of the GW calculation both including and excluding SOC effects. On the other hand, the feature at 3.5 eV is not visible in the calculated ELF, but the calculated DOS in Figure a indicates that this peak is associated to a electronic transition from the occupied d states to the unoccupied d states of W. The absence of the feature at 3.5 eV energy loss can be explained due to limitations of the GW approach. Indeed, GW accurately describes a single-particle process; however, for optical excitations the effect of interactions between holes and electrons needs to be taken into account. Properly describing electron–hole interactions could be achieved by solving the Bethe-Salpeter Equations (BSE).[45−47] Further details about the origin of these two observed peaks in the low-loss region can be found in the JDOS discussion below.

Concerning the properties of the bulk plasmon, it is found to be located at around 23 eV in both the experimental EELS measurements and in the DFT and GW calculations, though in the latter case it appears at somewhat smaller energies. The minor difference between the calculated and experimentally measured positions of the bulk plasmon peak of the mixed 2H/3R polytype can be attributed to the exclusion of local field effects and nonzero momentum transfer effects while calculating the dielectric response of the material, as implemented in the optic package.[48−50]

Joint Density of States

We consider now the results for the joint density of states (JDOS) calculated at the GW level with and without taking spin–orbit coupling effects. We also evaluate the contribution of separate bands to the joint density of states to elucidate which ones are important in specific energy losses regions. We note that in the cases including spin–orbit coupling, we double the amount of bands needed with respect to the case without spin–orbit coupling to properly describe the energy-loss spectra. To properly describe the bulk plasmon peaks appearing in the energy-loss regions in the range between 10 to 30 eV range, a larger amount of bands needs to be incorporated into the joint density of states calculation; this is consistent with the concept of plasmons arising as collective excitations of the electrons in the material.

If one starts with the calculation of the joint density of states by including one valence band (VB), the top valence band, and one conduction band, specifically the bottom conduction band (CB), no contributions are observed, see Figure . As then one includes more valence and conduction bands in addition to the top valence band and the bottom conduction band, we can clearly assess the contributions of the different bands with respect to the features observed in the low-loss region (see Figure a–d). In particular, one can reproduce the features at 3 and 3.5 eV when adding the contributions from a total of 20 VBs and CBs in the JDOS calculation. The feature at around 4 eV does not arise in the experimental EELS data, perhaps due to limitations in the energy resolution. In addition, features at 5 and 6 eV are also observed in the JDOS but not in the experimental EELS measurements.

Figure 4

Calculated joint density of states of the 2H/3R crystal structure with the GW framework without (a,b) and with (c,d) spin–orbit coupling effects taken into account. Here VB and CB stand for valence band and conduction band, respectively.

In the energy-loss region below 10 eV, we observe that the top 20 valence bands and lower 20 conduction bands can almost completely describe the energy-loss region for the 2H/3R mixed phase. To describe the energy-loss region for energy loss near the bulk plasmon, all the conduction and valence bands available are needed, consistent with the collective excitation behavior of bulk plasmons. Furthermore, as we proceed toward smaller energy losses we can see that a decreasing amount of valence and conduction bands are needed to explain the peaks in the corresponding energy region for all the polytypes.

Discussion

Our first-principles calculations for the value of the band gap and the location of the bulk plasmon peak for the 2H crystal structure of WS2 are found to be in good agreement with previous theoretical and experimental work.[51,52] Therefore, we are confident that we can reliably apply the same theoretical infrastructure to predict the electrical and optical properties of the corresponding 2H/3R polytype.

In the case of 2H/3R polytypism, the location of the bulk plasmon peak ascertained in WS2 nanostructures by means of EELS[22,25] is rather similar to that of the bulk plasmon peak in the 2H polytype, in agreement with the theoretical predictions in this work. Furthermore, here we have also calculated that the 2H/3R polytype should exhibit an indirect band gap with a value in the range between 1.40 and 1.48 eV. This prediction is in excellent agreement with the experimental result reported in refs (22 and 25) in which a band gap value of 1.6–0.2+0.3 eV was extracted after subtracting the EELS zero-loss peak from the EELS data using machine learning techniques.

Finally, in the case of the 3R polytype we have calculated an indirect band gap with a value in between 1.39 and 1.54 eV. The large discrepancy for the band gap calculations in the cases with and without taking into account the spin–orbit coupling effects could be related to an insufficiently dense k-mesh used in the GW calculation for the 3R polytype. For this crystal structure, we also calculate a bulk plasmon peak, the position of which is close to that of the 2H and 2H/3R polytypes. To the best of our knowledge, no experimental measurements of the position of the bulk plasmon peak, the band gap type, or the band gap value have been reported for the 3R polytype, and hence it is not possible to compare our predictions to experimental data.

Conclusions

In this work we have carried out ab initio calculations using Density Functional Theory and the GW approximation of the optoelectronic properties of the 2H, 3R, and 2H/3R polytypes of WS2. We have demonstrated how the band gap value of the 2H/3R polytypism lies between the band gap values of the 2H and 3R polytypes, where it is observed to be closer to the 2H band gap value compared to the 3R band gap value. Furthermore, this first-principle calculation is in good agreement with corresponding experimental determination from EELS measurements. Comparable band structures were found of the 2H and 2H/3R polytypes, where the top of the valence band of the 2H/3R polytypism lies at the A high symmetry point. We have shown that the bulk plasmon peaks of all the polytypes occurs at similar energy-loss values. For the energy-loss function, we have determined the contribution of the different valence and conduction bands to the energy-loss intensity for different energy-loss regions. Our results provide key input toward assessing the feasibility of TMD-based heterostructures composed of a single material with mixed crystalline phases.