One-dimensional electron gas in strained lateral heterostructures of single layer materials.

Confinement of the electron gas along one of the spatial directions opens an avenue for studying fundamentals of quantum transport along the side of numerous practical electronic applications, with high-electron-mobility transistors being a prominent example. A heterojunction of two materials with dissimilar electronic polarisation can be used for engineering of the conducting channel. Extension of this concept to single-layer materials leads to one-dimensional electron gas (1DEG). MoS2/WS2 lateral heterostructure is used as a prototype for the realisation of 1DEG. The electronic polarisation discontinuity is achieved by straining the heterojunction taking advantage of dissimilarities in the piezoelectric coupling between MoS2 and WS2. A complete theory that describes an induced electric field profile in lateral heterojunctions of two-dimensional materials is proposed and verified by first principle calculations.

Introduction

Confinement of electrons along one of the spatial directions results in a two-dimensional electron gas (2DEG) that exhibits interesting physical phenomena along the side of useful technological applications. Particular examples include the field of quantum transport and mesoscopic physics[1] as well as high-electron-mobility transistors that are used in integrated circuits as digital on-off switches[2]. The advantage of 2DEG conducting channel is the high mobility of charge carriers due to the absence of deleterious effects inherent to ionised impurity scattering that allows for ballistic transport[3]. Engineering of 2DEG conventionally requires the use of a modulation doping technique[4] as in the case of (AlGa)As/GaAs heterostructures. Alternatively, the 2DEG can be achieved in undoped structures with an extreme band bending induced by the strong electric field at a heterojunction between two dielectric materials with dissimilar electronic polarisation such as (AlGa)N/GaN interface[5, 6]. It is interesting to see whether polarisation effects in two-dimensional (2D) materials can be used to achieve confinements of electrons along one spatial direction?.

2D materials become a perspective avenue for keeping up with latest trends in miniaturisation of electronics, culminating in a demonstration of the single layer MoS2 transistor[7-9]. Unlike group III-nitrides, free-standing transition-metal dichalcogenides do not possess spontaneous polarisation due to symmetry arguments. However, single-atomic-layer h-BN and monolayer transition-metal dichalcogenides have been theoretically predicted[10] and experimental confirmed[11, 12] to show piezoelectricity as a result of strain-induced lattice distortions. Two types of heterostructures that involve 2D materials are discussed in the literature: (i) multilayer heterostructures produced by stacking of different 2D materials, so-called van der Waals heterostructures[13], and (ii) lateral heterostructures, which are formed when two materials are covalently bonded within the 2D plane[14].

It will be shown that a lateral heterojunction of 2D materials with dissimilar piezoelectric properties can be used to achieve additional confinement of charge carriers along the interface, which creates conditions for realisation of a one-dimensional electron gas (1DEG). A complete theory that describes an induced electric field profile in lateral heterojunctions of 2D materials is presented and verified by first principle calculations.

First-principle model

First, we will use an ab initio model to explore the feasibility of achieving conditions for 1D confinement of charge carriers in a lateral heterojunction of two single-layer materials. For this purpose, an 80-atoms MoS2/WS2 supercell is constructed as illustrated in Fig. 1(a). MoS2 and WS2 are chosen due to an almost identical lattice parameter of two materials (less that 0.1% mismatch), which reduces the misfit strain at the interface. One would expect the heterostructure to possess no built-in electric field since transition metal dichalcogenides manifest no net polarisation unlike group-III nitride bulk semiconductors. This hypothesis can be verified by plotting the potential energy across the heterojunction (Fig. 1b). The potential energy profile shows periodic oscillations with minima in the vicinity of nuclei and maxima corresponding to interstitial regions. It is evident that maxima of the potential energy remain constant within MoS2 and WS2 domains with an abrupt step-like transition at the interface. The confinement of charge carriers resembles that in a quantum well (Fig. 1d).

Figure 1

Electron confinement in lateral MoS2/WS2 heterojunction. (a) 80-atoms model of the heterojunction. (b,c) Electrostatic potential energy profile across the heterojunction without strain and with the strain of ϵ 1 = 0.1, respectively. The scan is taken between points with the fractional coordinates (0, 1/2, 0) and (1, 1/2, 0). The built-in electric field corresponds to a macroscopic slope of the potential energy. (d,e) The electron wavefunction amplitude |ψ e()|2 represents the lowest unoccupied state in unstrained and strained heterostructures, respectively. The strain-induced electric field confines electrons forming a one-dimensional conducting channel along the MoS2/WS2 interface. The band diagrams show the spatial evolution of the conduction band edge (CBE) schematically to assist with interpretation of the wavefunction plot.

Next, the same heterostructure is uniformly strained in the direction perpendicular to the heterojunction, i.e., along x-axis (Fig. 1a). The magnitude of strain is deliberately chosen high (10%) in order to magnify observed effects. The Poisson’s contraction is simulated by relaxing the second lateral dimension of the supercell to eliminate the macroscopic stress σ 22, accompanied by a full relaxation of internal degrees of freedom. It is found that, after relaxation, the macroscopic strain of 10% is non-uniformly distributed among both material domains. The effective strain in MoS2 is 10.5%, while WS2 accommodates only 9.5%. This result can be attributed to differences in stiffness between two materials.

It is also noticed that the external strain breaks 3-fold rotational symmetry, which is responsible for the absence of spontaneous polarisation in MoS2 and WS2 due to the cancellation of polarisation dipoles (Fig. 2). The symmetry breaking is evident from the disparity in Mo-S bond lengths: 2.52 Å vs 2.41 Å for the bonds oriented along or tilted with respect to the strain direction. The electrostatic potential profile plotted in Fig. 1(c) reveals the presence of an electric field in MoS2 and WS2 domains of approximately equal magnitude, but the opposite direction. The magnitude of electric field varies (±1%) depending on the coordinates of the line scan (see Supplementary information for more details); the average field is approximately 8.2 ± 0.5 mV/Å. The created saw-like potential confines charge carriers in the vicinity of the MoS2/WS2 interface (Fig. 1e) producing a narrow 1D conduction channel along x-axis of the width a few interatomic spacings.

Figure 2

Strain-induced change in electronic polarisation of hexagonal MoS2. (a) Rectangular unit cell. (b) Cancellation of local dipoles p induced by of the charge transfer ±Z * due to the C3 rotational symmetry, which results in the vanishing of a spontaneous polarisation P 0. (c,d) Symmetry breaking due to uniaxial strain induces a macroscopic dipole moment giving rise to the strain-induced polarisation. (e) Under the equal plain strain condition (ϵ 1 = ϵ 2), the symmetry is preserved. Thus no change in polarisation should be observed.

Qualitatively, an origin of the electric field can be attributed to heterogeneity in polarisation induced by the strain in MoS2 and WS2 domains (see Fig. 3). To gain a quantitative understanding of the observed effects in 2D materials, a model that couples continuum mechanics and Poisson equation is developed below.

Figure 3

Continuum model of lateral MoS2/WS2 heterointerface. (a) The strain along x-axis induces a heterogeneity of electronic polarisation P due to differences in the piezoelectric response between two materials. (b) Inhomogeneities in polarisation create regions with an opposite electric field E that results in subsequent spatial bending of the conduction band minima (CBM) and valence maxima (VBM).

Continuum model

The purpose of this model is to describe the electric field profile induced due to piezoelectric effects in 2D strained heterostructures. The problem is similar to that solved by Ambacher et al.[6] for AlGaN/GaN heterostructures, however, there are peculiarities related to 2D character of the materials in question, which warrant repeating some basic steps.

The free electro-elastic energy density stored in a linear medium can be expressed as[15] where ϵ = (ϵ 1, ϵ 2, ϵ 6) are components of the strain tensor written in the Voigt’s matrix notations, E is the electric field projection along i axis, C are components of the stiffness matrix, ε are components of the electrical permittivity tensor of the material, and the range of indices i, j = 1, 2, 6, l, m = 1, 2 is adapted to 2D. Oftentimes, the macroscopic strain is found by minimising the elastic energy only[6] (first term in Eq. (1)). However, it should be emphasised that the electric field and strain are coupled through the electric displacement, which takes the formHere 0 is the permanent (spontaneous) polarisation and e are components of piezoelectric strain tensor. In the absence of free charges, the Gauss’s law requires

This implies continuity of the electric displacement at the interface of two domains (see Fig. 3a)which includes contributions from permanent, strain-induced, and field-induced electric dipoles in the material. The strain (r) and electric field (r) distributions can be found by minimising the total electro-elastic energysubject to boundary conditions, e.g., an applied macroscopic strain.

2D materials pose a challenge related to defining the integration volume required to evaluate the total free energy in Eq. (5). There are attempts in the literature[16] to assign an effective thickness to atomically thin monolayers to compare their properties (strength, elastic modulus, or dielectric constant) to bulk materials. However, such analysis always bares the element of ambiguity. Alternatively, it seems more logical for 2D materials to use area rather than volume for normalising their specific properties. As a result, the stiffness coefficients C acquire units of N/m, whereas the piezoelectric coefficients e are expressed in units of C/m in 2D[10]. To remain consistent, an effective 2D dielectric permittivity ε needs to be defined. Then Eqs (1)–(5) can be readily extended to 2D materials, provided the free energy in Eq. (5) is integrated over the surface area, which eliminates ambiguities associated with the layer thickness.

Structural, elastic, piezoelectric, and dielectric properties of monolayer MoS2 and WS2 are gathered in Table 1. The structural unit and orientation of coordinate axes are illustrated in Fig. 2(a). The calculated lattice parameters are in agreement with experiment and other calculations reported in the literature. The hexagonal symmetry of a single layer (point group D3h) reduces the number of independent coefficients in the stiffness matrix down to two: C 11 and C 12 [17]. Our values of C 11 and C 12 listed in Table 1 agree with those obtained in previous DFT calculations. The piezoelectric tensor is characterised by a single independent element e 11, due to symmetry arguments. The calculated values agree well with prior theoretical studies. However, approximately 20% deviation from existing experimental data is observed. This deviation is acceptable giving the large uncertainty of experimental measurements.

Table 1

Structural parameters and effective 2D elastic, piezoelectric and static dielectric properties of single-layer hexagonal MoS2 and WS2 from self-consistent DFT calculations (relaxed-ion approximation).

Parameter	Units	MoS2		WS2
		Calculated	Other sources	Calculated	Other sources
a 0	Å	3.185	3.16a, 3.19b	3.188	3.15a, 3.19b
C 11	N/m	133	130b	146	144b
C 12	N/m	33	32b	32	31b
e 11	pC/m	359	290 ± 50c, 364b	249	247b
x 11	F	7.4·10−20	7.5·10−20,d	7.0 × 10−20	7.0 × 10−20,e

(1) Experimental[36].

(2) Calculated with DFT/GGA[10].

(3) Experimental[12].

(4) Obtained using Eq. (6) based on MoS2 bulk in-plane relative dielectric permittivity of 15 and the interlayer separation of 6.02 Å[37].

(5) Obtained using Eq. (6) based on WS2 bulk in-plane relative dielectric permittivity of 14 and the interlayer separation of 6.06 Å[37].

The static dielectric permittivity is one of the least studied properties of single-layer transition metal dichalcogenides. The present calculations yield the value of for the in-plane relative dielectric permittivity of a single-layer MoS2, with ε 0 being the permittivity of free space. It should be emphasised that ε 3D is an extensive property, which is determined by the thickness of the vacuum layer H that is used for separation between periodical images in the direction perpendicular to the planar structure. To represent a free-standing layer of MoS2, the value of H  = 24.6 Å was chosen, which is approximately by a factor of four greater than the spacing between layers in bulk. Berkelbach et al.[18] proposed evaluation of the effective 2D polarizability χ 2D of planar materials using the following relationshipwhich yields the effective in-plane polarizability of F, as compared to the value of F obtained for bulk MoS2 (see Table 1).

Potential energy profile scans similar to those shown Fig. 1 reveal the presence of a zig-zag electric field even in the middle of the vacuum region due to periodic boundary conditions along z-axis (see Supplementary information). To capture the energy stored in the vacuum due to the finite electric field, the effective 2D dielectric permittivity used in calculation of the free energy density in Eq. (1) is expressed as

The additional term ε 0 H contributes approximately 25% to the value of ε 2D.

Minimization of the total free energy W for the 2D strained lateral heterostructure of MoS2 and WS2 was performed using a Lagrange multiplier approach with respect to the strain tensor I,II and electric field I,II in both domains (see Methods for details). The quasi-2D continuum model with material parameters listed in Table 1 yields the strain distribution of and , which is in excellent agreement with DFT results. The greater strain in MoS2 (domain I) is due to its lower stiffness C as compared to WS2 (see Table 1). The continuum model also properly captures magnitude of the electric field |E| = 8.2·107 V/m, which coincides with the average slope of the electrostatic potential profile obtained from first-principle calculations.

Finally, we would like to comment on a practical realisation of the strained heterostructures discussed in this paper. MoS2/WS2 lateral heterostructures usually have a morphology of equilateral triangular flakes of the size of a few micrometres[14, 19]. MoS2 forms an inner core surrounded by the WS2 outer layer[17]. Gong et al.[19] reported achieving an atomically sharp MoS2/WS2 in-plane interface. The interface is preferentially formed along “zigzag” direction (the y-axis in Fig. 3a), which is consistent with the structural model studied here. The strain can be applied employing a setup shown in Fig. 4 previously used by Conley et al.[20] to measure the band gap shift of MoS2 with strain. The method involves clamping of a specimen at the surface of a mechanically bent substrate, which allows applying of a uniform strain up to 2% in a highly controlled manner. The strain magnitude much less than 10% can be sufficient giving a much larger length of real heterostructures in comparison to that modelled here. The presence of a strain-induced electric field can be verified by measuring a photoluminescence (PL). In unstrained MoS2/WS2 lateral heterostructures, the PL intensity is enhanced at the MoS2/WS2 interface[14, 19] due to the type-II band alignment[21]. The PL intensity at the interface that develops 1DEG is expected to diminish when the strain is applied due to the induced electric field that separates charge carriers.

Figure 4

(a) Schematic illustration of a four-point bending setup for straining a triangular MoS2/WS2 lateral heterosctructure. (b,c) 1DEG is formed at the zigzag interface oriented perpendicular to the applied strain.

Conclusions

One-dimensional conductivity channel is obtained in a lateral MoS2/WS2 heterojunction. Conducting electronic states are confined along the interface by an inhomogeneous electric field that is induced by differences in the piezoelectric and elastic response of two materials thereby creating a one-dimensional electron gas. An effective model that captures interactions between electric and elastic degrees of freedom in low-dimensional heterostructures is developed. The model accurately predicts the magnitude of macroscopic electric field induced in the strained heterostructure as verified by ab initio calculations. This realisation of 1D electron gas creates an alternative to a quasi 1D conducting channel formed in the 2D electron gas of GaAs/(AlGa)As heterostructures by electrostatic gating[22, 23] that can be potentially used for low-power switching applications.

Methods

Calculation of structural, elastic, and dielectric properties

Electronic structure calculations of single-layer hexagonal MoS2 and WS2 have been performed in the framework of the density functional theory (DFT)[24] using Perdew-Burke-Ernzerhof generalized gradient approximation (GGA-PBE) for the exchange-correlation functional[25]. Structural, elastic, and dielectric properties were modelled using the Vienna ab initio simulation program (VASP) and projector augmented-wave (PAW) potentials[26-28]. The structure was represented by a single layer of MoS2 or WS2 with a vacuum separation, which is approximately equal to a quadruple value of the equilibrium spacing between layers of the bulk 2H-MoS2. The structural optimisation was carried out in conjunction with relaxation of the in-plane lattice parameter α. The structure was considered optimised when the magnitude of Hellmann-Feynman forces acting on atoms dropped below 2 meV/Å. The Brillouin zone of the primitive unit cell was sampled using 16 × 16 × 1 Monkhorst-Pack grid[29]. The mesh was appropriately scaled when supercells are considered.

A hard PAW potential was used to represent sulphur (S_h). Semi-core electrons were included as valence electrons in molybdenum (Mo_sv) and tungsten (W_pv). The cutoff energy for a plane wave expansion was set at 500 eV, which is 25% higher than the value recommended in the pseudopotential file. The higher cutoff energy was essential for obtaining accurate, converged lattice parameters.

The elastic tensor was determined using a finite differences technique from the strain-stress relationship calculated in response to finite distortions of the lattice taking into account relaxation of the ions. The total of eight strained structures that represent various permutations of the strain ϵ 1,2 = {−0.02, 0, +0.02} were considered.

The relaxed-ion dielectric tensor was calculated using the finite external electric field of the magnitude 1 meV/Å. The tight energy convergence of 10−9 eV was required to achieve the accuracy of 0.1 for the relative dielectric permittivity.

Calculation of piezoelectric coefficients

Calculations of piezoelectric coefficients were performed using a full potential linear augmented plane wave method implemented in Wien2k package[30] in conjunction with BerryPI extension[31] that utilises a Berry phase approach[32] for computing macroscopic polarization. Piezoelectric strain coefficients are conventionally defined aswhere dP i is the change in macroscopic polarisation along i-axes observed in response to the increment in j’s strain component dϵ j. It seems straight forward to evaluate the coefficients using a finite difference, which involves computing the polarisation of strained and unstrained structures. However, this approach introduces complications related to the choice of a reference structure that must remain commensurate with the strained cell to serve as a reference. A similar approach was introduced by Posternak et al.[33] to assess the spontaneous bulk polarisation of wurtzite BeO, where the zinc blende structure served as a reference due to symmetry arguments.

In the case of hexagonal transition metal dichalcogenides, the polarisation of an unperturbed layer can be taken as a reference zero due to the cancellation of local dipoles resulted from the 3-fold rotational symmetry as illustrated in Fig. 1. Any strain tensor that preserves this symmetry (e.g., ϵ 1 = ϵ 2) produces no change in polarisation. This result translates into a symmetry of the piezoelectric coefficients[34] which is inherent to D3 point group. It turns out that no change in the Berry phase results from the strain ϵ 1 = ϵ 2. However, there is a sizeable change in polarisation originated from the increment in the cell volume that is incompatible with symmetry-imposed constraints in Eq. (9). To resolve this contradiction, the piezoelectric coefficients were redefined in terms of the Berry phaseHere a is the lattice parameter associated with the crystallographic axis i, Ω0 is the volume of the unperturbed unit cell, and ϕ is the Berry phase along direction i that includes both ionic and electronic components. A least square fit technique was used to calculate piezoelectric coefficients for the total of eight strained structures that represent various permutations of the strain (the same as for elastic properties). Additional relaxation of atomic positions was performed for each stained structure.

Visualization of atomic structures was performed using VESTA 3 package[35].

Free energy minimization

The objective is to find a set of strains and electric fieldsthat minimise the internal energy of the system W defined by Eq. (5) for a specific case of the strained lateral heterostructure shown in Fig. 3. The optimization is subject to constrains, such as an applied macroscopic strain ϵ 1 = 0.1, continuity of both the electric displacement (Eq. 4) and matter. From the mathematical standpoint, it is a constrained optimisation of an objective function represented by a quadratic form (Eq. 1). The problem can be solved using a method of Lagrange multipliers.

First, a matrix is constructed to represent linear coefficients of the partial derivatives ∂w/∂x k, where x k is any variable from the list (11). When strain variables in the first domain are concerned, the linear coefficients are simply components of the elastic stiffness matrixwhich is written taking into account symmetry imposed by the lattice. Similarly, the dielectric permittivity tensorrepresents the linear coefficients of the partial derivatives ∂w/∂x k for variables that correspond to the electric field components. Generalising to all optimisation variable related to the domain I, the matrix of linear coefficients takes the form

Our objective function is not the energy density w, but rather the total internal energy of the system W, which takes into account the individual area occupied by each domain. For the lateral junction of two domains that share the same width but may have different length L and L (Fig. 3), linear coefficients of the partial derivatives ∂w/∂x k form a matrix

Now the following boundary conditions need to be incorporated The first condition stems from Eqs (2) and (4) that capture the essence of piezoelectric coupling between the strain and electric field. It is implied that the spontaneous polarisation is zero in both materials (P 0 = 0) when unstrained. The second and third requirements account for the continuity of the heterostructure along the direction of the applied strain and perpendicular to that. The difference corresponds to a lattice mismatch between two materials. The left-hand-side of Eq. (16a–c) can be transformed into a matrix formwhere columns correspond to the optimisation variables in Eq. (11). The symmetry of piezoelectric strain coefficients (e 11 − e 12) is taken into account during this transformation.

Finally, the energy terms and constraints are combined in a matrix

that represents Lagrangian of the problem . Unknownsare obtained by solving a linear equationwith the right hand side being a column vector

The first ten elements of are zero due to the requirement of at the optimum for each variable listed in Eq. (11). The remaining elements represent the right hand side of Eq. (16). Here λ’s are Lagrange multipliers.

Data availability

Crystallographic information files (CIF) with atomic structures used in calculations can be accessed through the Cambridge crystallographic data centre (CCDC deposition numbers 1520213–1520216).

Supplementary information