Structural phase transitions in two-dimensional Mo- and W-dichalcogenide monolayers (2024)

Monolayer crystal structures

Under ambient conditions, all group VI TMDs (except WTe₂) are reported to exist in a layered bulk crystal structure composed of monolayers wherein the X atoms are in trigonal prismatic coordination around the M atoms¹⁷. The atomic stacking sequence within a single XMX monolayer is βAβ. In keeping with prior literature, we will refer to this as the 2H phase, even though the prefix ‘2’ is irrelevant in monolayers because it refers to a bulk stacking mode. Group VI MX₂ monolayers in the 2H structure are semiconducting with band gaps between 1 and 2 eV (refs 2, 33, 34). 2H TMDs are promising semiconductors for flexible electronics applications^{2,3,4,5,6,7,8}. The 2H structure gives rise to metallic edge states that are associated with electrocatalytic activity¹¹. The primitive unit cell of the 2H phase is hexagonal. For reasons of consistency between different phases, our calculations on 2H use a non-primitive rectangular unit cell whose axes align with zigzag and armchair directions of the structure. These special axes can be experimentally identified using second harmonic generation^35,36, and possibly also using the intrinsic piezoelectricity predicted to exist in these materials³⁷. Figure 1 shows the 2H structure within a rectangular unit cell having lattice constants a and b.

When one of the 2H structure’s X layers is shifted (for instance, βAβ→βAγ), the X atoms are in octahedral coordination around the M atoms, and the crystal becomes metallic. This phase is referred to as 1T and is observed in group IV and group V TMD compounds (for example, TiS₂ and TaSe₂ (ref. 17)). Its atomic structure is also shown in Fig. 1. We have calculated the atomic vibrational normal modes (Г-phonons) within the relaxed rectangular 1T unit cell of all six group VI TMDs. In all cases, one of the optical phonon modes has an imaginary vibrational frequency. This result asserts that the high-symmetry 1T structure is unstable (saddle point in 18D atomic potential energy surface), at least in the absence of external stabilizing influences.

The group VI TMDs do have a stable metallic structure with octahedral-like M–X coordination. This lower-symmetry phase, which we will refer to as 1T′, is a distorted version of the 1T structure^12,17,38. This phase has recently been shown to enhance electrocatalytic activity in WS₂ (ref. 12). A rectangular (as well as primitive) 1T′ unit cell is displayed in Fig. 1. The 1T′ phase is observed in WTe₂ under ambient conditions¹⁷, MoTe₂ at high temperature³⁹ and as a metastable phase in instances of chemically exfoliated¹⁸ and restacked³⁸ MX₂ monolayers.

Phase energetics in the absence of mechanical stress

Figure 2 shows the calculated equilibrium (that is, stress-free) relative energies of the phases (2H, 1T and 1T′) of the group VI MX₂ monolayer materials. These values are calculated using DFT under the generalized gradient approximation (GGA) for electronic exchange and correlation effects, using the Perdew–Burke–Ernzerhof⁴⁰ (PBE) functional. The results are consistent with the experimental observation for bulk crystals that only WTe₂ is 1T′ under ambient conditions¹⁷. It is also evident from these data that the energy associated with the 1T to 1T′ relaxation is considerable: several tenths of an eV per MX₂ formula unit. The three phases’ equilibrium lattice constants and energies are tabulated in Supplementary Table 1.

Phase energetics under biaxial strain

Thermodynamics asserts that a system will seek to minimize whichever thermodynamic potential is appropriate for the prevailing mechanical and thermal boundary conditions⁴¹. The simplest example of such a thermodynamic potential is the internal energy U. In the low-temperature limit, any system will seek to minimize U when it is constrained to a given shape. Under these conditions, an MX₂ monolayer constrained to be described by a rectangular unit cell with dimensions a × b is expected to be in the lower-U phase for those values of a and b.

Experimentally, relevant phase diagrams of monolayers differ from those of bulk materials at high pressure in at least one important respect: the monolayer can be mechanically coupled to a substrate with friction, enabling the independent control of a and b lattice parameters (Fig. 3). Another key distinction is that large elastic deformations in monolayers can be reached through tensile strain^42,43, whereas large elastic deformations in most bulk materials are accessible only under compression. Large compressive stresses are problematic in TMD monolayer materials due to the spontaneous ripple formation that has been studied in MoS₂ (refs 44, 45).

For the six group VI MX₂ monolayers, we use GGA–DFT to calculate the energies U(a,b) of the three monolayer crystal structures on a 7 × 7 grid in (a,b) space, giving a total of 49 points of (a,b)-values around the minimum-energy equilibrium lattice constants a₀ and b₀ (listed in Supplementary Table 1). U is obtained after allowing the ions to relax their positions within each unit cell. Intermediate values for each phases’ U(a,b) are subsequently approximated using the Lagrange⁴⁶ interpolation method:

Lagrange interpolation is chosen because it contains no physical assumptions about the shape of the U(a,b) energy surface over a large range of tensile and compressive strains. It also greatly facilitates the approximation of local derivatives (which we need later) without suffering from conditioning problems found in other high-order polynomial methods.

Using this approach, we discover that the 2H and 1T′ U(a,b) energy surfaces intersect for sufficiently large strains. Figure 4 shows the contours that follow the intersection of the 2H and 1T′ energy in (a,b) space. The changes in a and b required to change the relative energies U of the 2H and 1T′ phase range from 13% (MoS₂) to 3% (MoTe₂). Because many bulk materials begin to dissipate strain energy through fracture or dislocations at strains on the order of 0.1% (roughly 10⁻⁵ eV per atom), these threshold strains may at first appear to be prohibitively large. However, it has been experimentally demonstrated that monolayer TMDs are exceedingly strong: Bertolazzi et al.⁴² have performed atomic force microscopy (AFM) experiments where 2H-MoS₂ monolayers are shown to reversibly withstand in-plane tensile stresses up to 15 Nm⁻¹, corresponding to ~10% of the material’s in-plane Young’s modulus. Such deformations correspond to an elastic energy of order 0.1 eV per atom. From the local derivatives of the 2H phase’s U(a,b) at the equibiaxial transition strain, we find that the 2H stresses (in Nm⁻¹) at this point are 12.8 (MoS₂), 10.8 (MoSe₂), 6.9 (MoTe₂), 13.6 (WS₂) and 10.5 (WSe₂). This suggests the possibility that a transition between 2H and 1T′ might be observable below but near the breaking threshold.

From Fig. 4, we can also see that WTe₂, which is usually in the 1T′ phase, can be pushed into a 2H regime under compression, which is complementary to all the other cases where one would go from 2H to 1T′ through tension. While this is interesting, we will not focus on it since in-plane compression in monolayer WTe₂ may be experimentally challenging to achieve without incurring any buckling response.

From the contours in Fig. 4, it is also apparent that a 2H-1T′ transition might be most easily accessible in MoTe₂ under tension along the b axis. The remainder of the Results section will find that the predicted phase boundary for MoTe₂ moves even closer to ambient conditions when including thermal effects, more advanced treatments of electronic exchange and correlation, and by generalizing the governing thermodynamic constraints beyond the case of fixed lattice constants.

Thermal corrections in MoTe₂

The data shown in Figs 2 and 4 derive from the DFT-calculated potential energy U=U_crystal, and therefore omit the vibrational component of free energy. This vibrational component can be important when the energy difference between phases is <0.1 eV, on the order of k_BT. Figure 2 shows that the 1T′-2H energy offset thus calculated for MoTe₂ is sufficiently small (43 meV) such that vibrational effects could have a role in the energetic ordering of phases. For a temperature-controlled experiment, the Helmholtz free energy A(a,b,T) replaces U(a,b) as the relevant thermodynamic potential. We perform this calculation for MoTe₂ by treating the DFT-calculated vibrational normal modes as quantum mechanical harmonic oscillators⁴⁷. We shall see that these effects in MoTe₂ have a considerable impact, even in the idealized T→0 K case because of contributions from the vibrational zero-point energy, which can be regarded as a manifestation of the partially wavelike nature of the atomic nuclei.

In the case of 2H- and 1T′-MoTe₂, we calculate a temperature-dependent vibrational free-energy correction based on the frequency spectrum of the 15 nonzero Г-point phonons belonging to the rectangular cell of each phase. These phonon frequencies are obtained by applying linear-response theory to the forces in perturbed unit cell configurations for all 7 × 7 (a,b) grid points. We then add the quasi-harmonic vibrational free-energy correction A_vib=U_vib–TS_vib to the potential energy U_crystal:

Since this approach samples 15 discrete values of ω_i at Г, it does not take into account the full dispersion of phonon frequencies for arbitrary wave vectors in the Brillouin zone. Large supercell calculations that sample the vibrational spectrum more finely (up to 159 values of ω_i in a 3 × 3 supercell) reveal that the use of 15 phonons in equation (2) leads to errors under 3 meV per MoTe₂ formula unit when comparing 2H and 1T′ free energies (Supplementary Fig. 1). This error is significantly smaller than the 1T′-2H’ U_crystal energy difference of 43 meV.

Proceeding with the results of equation (2), the Lagrange interpolation procedure is carried out again at intermediate values of a and b. The top three curves in Fig. 5 show that the vibrational correction A_vib moves the PBE-calculated (a,b)-space free-energy crossing in MoTe₂ to smaller strains. Indeed, a 0-K zero-point free-energy correction lowers the threshold strain by up to 0.01 in some regions, and increasing the temperature from 0 to 300 K produces another shift of up to 0.02 in the same regions.

Kinetic aspects of the MoTe₂ phase transition

Thus far, our purely thermodynamic analysis does not provide information regarding the timescales required to observe a transition between the 2H and the 1T′ phase. Experimental observation of a thermodynamically allowed transformation might not be feasible if a large kinetic barrier renders it inaccessible. A Climbing-Image Nudged Elastic Band⁴⁸ calculation within GGA–DFT reveals a 2H-1T′ energy barrier of 0.88 eV per formula unit in MoTe₂ at the equilibrium lattice constants of the 2H phase. In the Arrhenius kinetics picture with a characteristic frequency of 10 THz, the timescale associated with this barrier is 50 s at room temperature. Although other factors such as interfaces, substrate, temperature, strain and impurities are likely to alter the transformation kinetics, this timescale suggests that this transformation is likely to be observable in the room temperature laboratory.

Hybrid functional applied to MoTe₂

Like many contemporary computational materials studies, our results presented thus far rely heavily on the generalized gradient PBE⁴⁰ functional for electronic exchange and correlation effects. The development of new exchange–correlation functionals is a highly active field, and a Hybrid PBE/Hartree–Fock approach known as the Heyd–Scuseria–Ernzerhof HSE06 (ref. 49) functional has recently shown superior agreement with experimental results in structural metal–insulator transition phase boundaries in Si⁵⁰. To explore how HSE06- and PBE-based predictions differ in a MoTe₂ context, we recalculate the U_crystal component of equation (2) using the HSE06 functional on the 7 × 7 grid of PBE-relaxed geometries. The HSE06-based phonon spectrum is not recalculated for this work because of the formidable computational demands presently posed by high-quality HSE06 calculations. Instead, the PBE-calculated frequencies are reintroduced on top of the HSE06-calculated U_crystal. The bottom three curves in Fig. 5 show that use of the HSE06-calculated U_crystal brings the 2H-1T′ threshold strains even closer to the origin. At 300 K, the transition is predicted to be within 2% strain of the equilibrium 2H lattice constants.

Alternative thermodynamic constraints

Some careful consideration of the relevant thermodynamic constraints is warranted here for the case of monolayers. The thermodynamic constraint of fixed lattice constants used in Figs 4 and 5 is an unusual one (different from constant volume or area) that is perhaps most applicable when crystallographic coherence is maintained with a strongly binding substrate. At some fixed temperature T, the Helmholtz free energy A(a,b,T) is the thermodynamic potential whose minimization determines which crystal phase (or coexistence of phases) will exist in a crystal possessing a rectangular unit cell with dimensions a × b. However, this fixed-cell constraint certainly need not apply when friction with the substrate is weak, or when the monolayer is freely suspended. In analogy with the important distinction between constant-pressure and constant volume experiments in bulk materials, we now show how to generalize our MoTe₂ results to other thermodynamic constraints with stable two-phase coexistence regions.

Perhaps the simplest thermodynamic constraint occurs when both in-plane principal tensions are equal, such that σ_xx=σ_yy. This isotropic-tension condition is analogous to the isotropic constant-pressure case in three dimensions. The natural thermodynamic potential governing this isotropic-load system (that is, the potential that is minimized) is a ‘hydrostatic’ Gibbs free energy G_hydro, where the surface tensions σ_xx=σ_yy≡σ and σ_xy=0:

The previously independent variables (a,b,T) used in A=A(a,b,T) are now a function of (σ,T) through the definition of a 2D hydrostatic contour {a,b}={a(σ,T),b(σ,T)}. This contour is determined directly from the interpolated A(a,b,T) surface and its local derivatives (for example, using σ_xx=(1/b)∂A/∂a).

Another physically relevant constraint is that of constant area. This is closely related to a constant volume constraint for bulk materials. Such a constraint can lead to mixed-phase regimes. Constant area is a macroscopic constraint that might be applicable when the edges of a freely suspended monolayer are clamped in place to fix the area, independent of stress.

A very useful constraint applies when one of the principal tensions is zero, for example, when σ_xx=σ_xy=0 and σ_yy is nonzero. This condition might hold for a ribbon suspended over a trench. Considering the 2H-1T′ energy landscape in Fig. 4, it would appear promising to apply a uniaxial stress along the b axis of MoTe₂ in order to observe a phase transition. We therefore study a system that is subjected to a specified uniaxial-load F_y=∂A/∂b along the crystal’s special y axis, whereas the x-face is treated as a free surface (that is, ∂A/∂a=0). Applying the appropriate Legendre transform to A for this case yields another Gibbs-like free energy G_y that acts as the governing thermodynamic potential.

In our data, the independent variables (F_y,T) can be mapped to a uniaxial-load contour {a,b}={a(F_y,T),b(F_y,T)}, derived from the interpolated A(a,b,T) surface and its local derivatives. This approach carries the microscopic assumption that both the 2H and the 1T′ phases’ b axes point in the y direction.

A closely related case is that of fixed macroscopic uniaxial strain. In the macroscopic uniaxial-strain case, one assumes that the lattice constant a remains fixed at the same value for both phases, such that strain occurs along both phases’ b axis. Macroscopic uniaxial strain allows for two-phase coexistence regimes with crystallographic coherence between the 2H and 1T′ phase. Macroscopic uniaxial strain is also relevant if the substrate–TMD interaction is strongly anisotropic. The strain ε_yy=‹b›/b₀−1 is termed macroscopic because b is the weighted average over two different lattice constants ‹b› when strained 2H and 1T′ phases coexist. In an experiment, changes in b are proportional to the macroscopic extension of the sample.

Additional thermodynamic ensembles may be applicable for monolayers. For example, when a monolayer is weakly bound to a substrate with friction, the atoms are allowed to move to some limited extent and restricted by contact with the surface. The thermodynamic potential of the TMD monolayer will be some intermediate case between the constant-force case (applicable to a frictionless substrate) and the fixed-(a,b) case where friction is large enough to inhibit significant movement of the metal atoms.

Alternative thermodynamic constraints applied to MoTe₂

Both the hydrostatic and uniaxial-load contours are displayed in Fig. 6a (PBE-calculated U_crystal) and Fig. 6b (HSE06-calculated U_crystal). When increasing F_y or σ, at some point the Gibbs free energies of the 2H and 1T′ phases cross. When that happens, the resulting phase transition can be seen as a discontinuous jump in (a,b)-space from the 2H to the 1T′ contour. It is interesting to note that at precisely this transition load, any mixed-phase coexistence of both phases ranging from 100% 2H to 100% 1T′ is thermodynamically stable. Figure 7 illustrates the general concept of mixed-phase regimes when a monolayer is progressively tensioned, for example, by increasing F_y or ‹b› (for the uniaxial stress/strain case), or by increasing hydrostatic stress σ or the unit cell area ab.

In the uniaxial-load case for MoTe₂, the calculated strains ε_yy=b/b₀–1 required to enter the 2H end of this transition regime are shown in Fig. 8a. It is striking that the choice of exchange–correlation functional and the treatment of atomic vibrations have a large impact on this number: ε_yy ranging from 0.5% (HSE06, 300 K) to 2.4% (PBE, no vibrational free energy). Figure 8b shows that the strain on the 2H end of the hydrostatic transition is similarly sensitive to temperature and electronic exchange–correlation functional. The discrepancy between PBE and HSE06 shows a clear path towards experimentally addressing the question of which functional is better suited for 2D materials studies.

The accuracy of these predictions of small transition strains based on 2D interpolation is validated using additional calculations in the case of MoTe₂. We run a 10 × 10 grid of PBE, phonon and HSE06 calculations in a region in (a,b)-space with axial strains ranging from −5 to 10%. This new data set is three times as fine as the 7 × 7 grid and also includes PBE-calculated phonons. Coexistence regimes are extracted from a series of 1D fits on this refined grid, as illustrated in Supplementary Fig. 2. The typical difference between uniaxial-load transition strains from this approach and the equivalent strains from the 2D Lagrange-interpolating method used throughout the main text is 0.002 for HSE06 and 0.005 for PBE (Supplementary Fig. 3).

In the case of macroscopic uniaxial strain at 300 K (not shown in Fig. 8), the HSE06 functional predicts an onset of the coexistence regime at 1.6% strain, whereas PBE predicts 3.0%. These strains are obtained from 1D cubic fits (illustrated in Supplementary Fig. 4) on the refined 10 × 10 grid of DFT-calculated energies, at fixed lattice constant a=3.550 Å and varying b. Supplementary Figure 5 summarizes all uniaxial-load and uniaxial-strain results based on this 10 × 10 grid.