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.

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 MoS 2 and WS 2 domains with an abrupt step-like transition at the interface. The confinement of charge carriers resembles that in a quantum well (Fig. 1d).
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 MoS 2 is 10.5%, while WS 2 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 MoS 2 and WS 2 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 MoS 2 and WS 2 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 (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 (r)| 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 MoS 2 /WS 2 interface. The band diagrams show the spatial evolution of the conduction band edge (CBE) schematically to assist with interpretation of the wavefunction plot.
Qualitatively, an origin of the electric field can be attributed to heterogeneity in polarisation induced by the strain in MoS 2 and WS 2 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.
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 i is the electric field projection along i axis, C ij are components of the stiffness matrix, ε lm 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 form

∑∑ ∑∑
Here P 0 is the permanent (spontaneous) polarisation and e li are components of piezoelectric strain tensor. In the absence of free charges, the Gauss's law requires I I I which includes contributions from permanent, strain-induced, and field-induced electric dipoles in the material. The strain ϵ(r) and electric field E(r) distributions can be found by minimising the total electro-elastic energy subject 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 MoS 2 and WS 2 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 D 3h ) 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.
The static dielectric permittivity is one of the least studied properties of single-layer transition metal dichalcogenides. The present calculations yield the value of ε ε = .
/ 45 for the in-plane relative dielectric permittivity of a single-layer MoS 2 , 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 V that is used for separation between periodical images in the direction perpendicular to the planar structure. To represent a free-standing layer of MoS 2 , the value of H V = 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 relationship   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 V contributes approximately 25% to the value of ε 2D . Minimization of the total free energy W for the 2D strained lateral heterostructure of MoS 2 and WS 2 was performed using a Lagrange multiplier approach with respect to the strain tensor ϵ I,II and electric field E 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 = . , which is in excellent agreement with DFT results. The greater strain in MoS 2 (domain I) is due to its lower stiffness C as compared to WS 2 (see Table 1). The continuum model also properly captures magnitude of the electric field |E| = 8.2·10 7 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. MoS 2 /WS 2 lateral heterostructures usually have a morphology of equilateral triangular flakes of the size of a few micrometres 14,19 . MoS 2 forms an inner core surrounded by the WS 2 outer layer 17 . Gong et al. 19 reported achieving an atomically sharp MoS 2 /WS 2 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 MoS 2 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 MoS 2 /WS 2 lateral heterostructures, the PL intensity is enhanced at the MoS 2 /WS 2 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.

Conclusions
One-dimensional conductivity channel is obtained in a lateral MoS 2 /WS 2 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 MoS 2 and WS 2 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][27][28] . The structure was represented by a single layer of MoS 2 or WS 2 with a vacuum separation, which is approximately equal to a quadruple value of the equilibrium spacing between layers of the bulk 2H-MoS 2 . 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 as where 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 11 12 which is inherent to D 3h 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 phase Here a i is the lattice parameter associated with the crystallographic axis i, Ω 0 is the volume of the unperturbed unit cell, and φ i 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.
Free energy minimization. The objective is to find a set of strains and electric fields 1  I  2  I  6  I  1  I  2  I  1  II  2  II  6  II  1  II  2  II that 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 matrix 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 I and L II (Fig. 3), linear coefficients of the partial derivatives ∂w/∂x k form a matrix (15) Now the following boundary conditions need to be incorporated (16a)   11  I  1  I  12  I  2  I  11  I  1  I  11  II  1  II  12  II  2  II  11  II  1  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 − a a 0 II 0 I corresponds to a lattice mismatch between two materials. The left-hand-side of Eq. (16a-c) can be transformed into a matrix form  11  I  11  I  11  I  11  II  11  II  11  II   I  I I   0  I  0  II where 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 (18) that represents Lagrangian of the problem . Unknowns  T  1  I  2  I  6  I  1  I  2  I  1  II  2  II  6  II  1  II  2  II  1  2  3       are obtained by solving a linear equation 20) with the right hand side being a column vector L a a ( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ( ) , ) T 1 I I I 0 II 0 I  The first ten elements of  are zero due to the requirement of  ∂ ∂ = x / 0 k 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).