Electronic and transport properties of TMDC planar superlattices: effective Hamiltonian approach

Two dimensional materials with atomic thickness such as graphene, Boron Nitride, phosphorene present themselves as promising materials for next generation electronics and photonics due to unique optoelectronic properties [1–19]. One of most important members of 2D materials are transition metal dichalcogenides (TMDCs) that has been emerged as complementary to graphene based electronics. 2D-TMDCs (MX₂, with M = Mo, W, etc., and X = S, Se, etc.) are quasi layered materials with strong covalent bonding which are typically produced through liquid-phase exfoliation from their bulk counterparts. Using them with graphene in device application is not also because of similar lattice geometry of monolayer TMDCs (ML-TMDCs) but also due to comparable electronic, optical, chemical, thermal and mechanical features with graphene [20–24].

2D features of TMDCs are basically different from bulk ones, such as crossover from indirect to direct band gap as the thickness drops to one monolayer which is composed of single layer transition metal atom (M) embedded between two atomic layers of chalcogenide atom (X) in a trigonal prismatic structure [16].

Inversion symmetry in these structures is broken due to occupation of two sublattices by single M atom and two X atoms. Besides, broken central symmetry causes strong spin-orbit coupling that splits valance bands and partially conduction bands in monolayer TMDCs [25]. Properties of TMDCs can be modulated by different methods for example, strain [26, 27], external electric field [28] and doping [29].

Recently, heterostructure of TMDCs provide a new route for tuning electronics properties [5, 30, 31]. Bringing different monolayer TMDCs features together in a heterostructures paves a path way for next generation nanoelectronics and nanophotonics devices [32, 33]. However, there is a drawback in Vander wales heterostructures composed of TMDCs monolayers that creates a challenge in implementing them for applications: direct band gap characteristic disappears in these structures. Among vertical heterostructures reported so far, only MoS₂/WSe₂ has direct band gap that easily switch to indirect band gap by very small disorder [34]. In fact, all MoS₂/MX₂ superlattices illustrate indirect band gap characteristics. In order to overcome this inefficiency, planner heterostructure has been attracting interests. Importantance of this new kind of heterostructures is demonstrated as essential needs to laterally integrating 2D materials in 2D electronics.

Planner superlattices with two dimensional confinements can be realized by stitching 2D semiconductors or alternatively sandwiching them between two materials in a plane [35–45]. In particular, structural similarities in TMDC semiconductors makes the high quality planar heterostructure with atomic sharp interfaces [36–42]. Multi-heterostructures or planar superlattices with almost a few nanometers width and compassing lots of interfaces in a small length has not been reported yet and has been remained a big challenge for researchers who are exploring high efficiency solar cells or other electronic devices based on 2D TMDCs [46].

Recently, embedded sub-nanometer 2D channels have shown a great potential of proceeding future electronics and enhancement of solar cells efficiency. Two research groups independently reported the growth of planar quantum wells with less than 2nm width and high quality inside TMDCs semiconductor [47, 48]. Partial difference in atomic bonding length between different neighboring monolayers of TMDCs creates different hexagons in MoS₂ and WSe₂ in the way that by increasing the number of adjacent dislocations in WSe₂, channels with width of less than 2nm are created and planar superlattice pattern is appeared including WSe₂ and MoS₂ [42, 47]. Therefore, it can be expected that similar growth mechanism can be employed and adapted for making 2D multi quantum wells systems or planar superlattices with the ability of tuning the width of ribbons. To the best of our knowledge, there are a few works on studying planar TMDCs superlattices, theoretically. For example, in Ref. [49] authors studied tunability of direct and indirect band gap of MoS₂/WSe₂ planar superlattice through density functional theory approach (DFT).

Motivated by these studies, in this paper we study the electronic band gap structure and transport properties of planar TMDCs structures with main focus on conductivity and Fano factor. By adopting transfer matrix formalism, we will obtain electronic band dispersion and wavefuctions from governed 2D Dirac like Hamiltonian. We will evaluate the effect of incident charge carriers angle, and geometrical parameters. In order to bring our model closer to reality [47, 48], we consider disorders and defect effects. We also find wavefunctions or probability amplitude that provides a tool to predict transport and even optical response of TMDCs superlattices. Our analytical approach not only is simple to follow but also gives results with very good consistency with that obtained by atomistic calculations within DFT.

Paper is organized as follows. The model superlattice and analytical method to get band dispersion and Fano spectrum are presented in section2 then in section 3 discussions on obtained results is provided in details. Finally, we will conclude and summarize our work in section 4.

2.1. Extended transfer matrix formalism

The behavior of electrons in ML-TMDCs is governed by 2D Dirac-like equation that reads [50]:

$$\begin{eqnarray}&&\hat{{\rm{{\rm H}}}}=at\left(\tau {k}_{x}{\hat{\sigma }}_{x}+{k}_{y}{\hat{\sigma }}_{y}\right)\otimes \hat{{\mathbb{I}}}+\displaystyle \frac{{\rm{\Delta }}}{2}{\hat{\sigma }}_{z}\otimes \hat{{\mathbb{I}}}-\displaystyle \frac{{\rm{\lambda }}{\rm{\tau }}}{2}\left({\hat{\sigma }}_{z}-\hat{1}\right)\otimes {\hat{S}}_{z},\end{eqnarray} \tag{ 1 }$$

${\hat{\sigma }}_{x.y.z}$ are well known Pauli psedu spin matrix and ${\boldsymbol{P}}=\hslash \left({k}_{x}.{k}_{y}\right)=\hslash \left(-i{\partial }_{x}.-i{\partial }_{y}\right)$ is momentum operator in x−y plane. $\hat{{\mathbb{I}}}$ is unit operator, a and t are lattice constant and effective hopping integral. Δ, 2λ, s and τ correspond to the direct band gap, spin-orbit splitting at the band extrema, s = +1(−1) indicate spin index and τ = +1(−1) is valley index attributed to K (K').

Spin orbit coupling (SOC) is appeared in single layer of TMDCs due to breaking symmetry results in removing degeneracy in band extrema and creating splitting that is significant in valence band extremum (about a few hundreds of meV). Third term in equation 1 shows SOC effect.

Hamiltonian operates on a two component pseudo-spinor wavefunction as ψ = (ψ_A , ψ_B ), where ψ_A = (x, y), ψ_B = (x, y) are the envelope functions for two sublattices in TMDCs. One can write ψ_A,B (x, y) = ψ_A,B (x)e^ikyy because of translational invariance in y direction. Equation (1) can be decomposed to two coupled following equations as:

$$\begin{eqnarray}&&\left\{\displaystyle \frac{d}{dx}{\psi }_{A}\left(x\right)-\displaystyle \frac{1}{\tau }{k}_{y}{\psi }_{A}\left(x\right)=i{k}_{A}{\psi }_{B}\left(x\right),\right.\end{eqnarray} \tag{ 2a }$$

$$\begin{eqnarray}&&\left\{\displaystyle \frac{d}{dx}{\psi }_{B}\left(x\right)+\displaystyle \frac{1}{\tau }{k}_{y}{\psi }_{B}\left(x\right)=i{k}_{B}{\psi }_{A}\left(x\right),\right.\end{eqnarray} \tag{ 2b }$$

where k_A,B is defined as ${k}_{A}=\tfrac{\left(E+\tfrac{{\rm{\Delta }}}{2}\right)-\lambda \tau s}{at\tau }$ and ${k}_{B}=\tfrac{\left(E-\tfrac{{\rm{\Delta }}}{2}\right)}{at\tau },$ E is incident energy of electrons. Figure 1 depicts schematic structure of our model for a planar superlattice. Regions A, B indicate two different ML-TMDCs that are MoS₂ and WSe₂ in our work. Structural parameters and incident angle of carriers are shown clearly.

Zoom In Zoom Out Reset image size

From equation (2) it is straightforward to obtain following equations:

$$\begin{eqnarray}&&\left\{\displaystyle \frac{{d}^{2}}{d{x}^{2}}{\psi }_{A}\left(x\right)+\left({k}_{j}^{2}-{k}_{y}^{2}\right){\psi }_{A}\left(x\right)=0,\right.\end{eqnarray} \tag{ 3a }$$

$$\begin{eqnarray}&&\left\{\displaystyle \frac{{d}^{2}}{d{x}^{2}}{\psi }_{B}\left(x\right)+\left({k}_{j}^{2}-{k}_{y}^{2}\right){\psi }_{B}\left(x\right)=0,\right.\end{eqnarray} \tag{ 3b }$$

${k}_{j}^{2}={k}_{A}{k}_{B}$ in jth layer. Assuming ${q}_{j}=\left\{\begin{array}{ll}\sqrt{{k}_{j}^{2}-{k}_{y}^{2}}, & {k}_{j}^{2}\gt {k}_{y}^{2}\\ i\sqrt{{k}_{y}^{2}-{k}_{j}^{2}},\, & otherwise\end{array}\right.$ is the x component of the wavevector in the jth layer. One can suggest following answers for set of two coupled equations (3a ) and (3b ) as:

$$\begin{eqnarray}&&\left\{{\psi }_{A}\left(x\right)=a{e}^{i{q}_{j}x}+b{e}^{-i{q}_{j}x}.\right.\end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray}&&\left\{{\psi }_{B}\left(x\right)=c{e}^{i{q}_{j}x}+d{e}^{-i{q}_{j}x}.\right.\end{eqnarray} \tag{ 4b }$$

a and c (b and d) are forward (backward) spinor components. From combining equations (4) and (2) it is straight to obtain wavefunctions ψ_A,B (x) as:

$$\begin{eqnarray}&&\left\{{\psi }_{A}\left(x\right)=a{e}^{i{q}_{j}x}+b{e}^{-i{q}_{j}x}.\right.\end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray}&&\left\{{\psi }_{B}\left(x\right)=a\displaystyle \frac{i{k}_{-j}}{i{q}_{j}+\tfrac{1}{\tau }{k}_{y}}{e}^{i{q}_{j}x}-b\displaystyle \frac{i{k}_{-j}}{i{q}_{j}-\tfrac{1}{\tau }{k}_{y}}{e}^{-i{q}_{j}x}.\right.\end{eqnarray} \tag{ 5b }$$

In order to connect the wavefunctions, ψ_A,B (x) at arbitrary points of x direction in jth layer, we will consider the basis as $\varphi \left(x\right)=\left(\begin{array}{c}{\varphi }_{1}\left(x\right)\\ {\varphi }_{2}\left(x\right)\end{array}\right)$ and assume that these basis can be written as:

$$\begin{eqnarray}&&\left\{{\varphi }_{1}\left(x\right)=a{e}^{i{q}_{j}x}+b{e}^{-i{q}_{j}x}.\right.\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&\left\{{\varphi }_{2}\left(x\right)=a{e}^{i{q}_{j}x}-b{e}^{-i{q}_{j}x}.\right.\end{eqnarray} \tag{ 6b }$$

Now, wavefunctions ψ_A,B (x). can be rewritten as:

$$\begin{eqnarray}&&\left(\begin{array}{c}{\psi }_{A}\left(x\right)\\ {\psi }_{B}\left(x\right)\end{array}\right)={R}_{j}\left(E,{k}_{y}\right)\left(\begin{array}{c}{\varphi }_{1}\left(x\right)\\ {\varphi }_{2}\left(x\right)\end{array}\right),\end{eqnarray} \tag{ 7 }$$

where:

$$\begin{eqnarray}{R}_{j}\left(E,{k}_{y}\right)=\left(\begin{array}{cc}1 & 0\\ \displaystyle \frac{i}{\tau }{P}_{j}\,\sin \,{\theta }_{j} & {P}_{j}\,\cos \,{\theta }_{j}\end{array}\right),\end{eqnarray} \tag{ 8 }$$

P_j and θ_j ${\theta }_{j}$ are ${P}_{j}=\tfrac{{k}_{Bj}}{{k}_{j}}$ and $\sin \,{\theta }_{j}=\tfrac{{k}_{y}}{{k}_{j}}$ $\left(\cos \,{\theta }_{j}=\tfrac{{q}_{j}}{{k}_{j}}\right),$ respectively. The relation between basis at two points of ${x}_{j}+{\rm{\Delta }}x$ and ${x}_{j}$ can be easily obtained as:

$$\begin{eqnarray}&&\left\{{\varphi }_{1}\left({x}_{j}+{\rm{\Delta }}x\right)=\,\cos \left({q}_{j}{\rm{\Delta }}x\right){\varphi }_{1}\left({x}_{j}\right)+i\,\sin \left({q}_{j}{\rm{\Delta }}x\right){\varphi }_{2}\left({x}_{j}\right).\right.\end{eqnarray} \tag{ 9a }$$

$$\begin{eqnarray}&&\left\{{\varphi }_{2}\left({x}_{j}+{\rm{\Delta }}x\right)=i\,\sin \left({q}_{j}{\rm{\Delta }}x\right){\varphi }_{1}\left({x}_{j}\right)+\,\cos \left({q}_{j}{\rm{\Delta }}x\right){\varphi }_{2}\left({x}_{j}\right).\,\right.\end{eqnarray} \tag{ 9b }$$

From equations (9) and (8), adopted or extended transfer matrix, ${M}_{j}\left({\rm{\Delta }}x,E,{k}_{y}\right),$ can be obtained that connects wavefunctions ${\psi }_{A,B}\left({x}_{j}\right)$ and ${\psi }_{A,B}\left({x}_{j}+{\rm{\Delta }}x\right):$

$$\begin{eqnarray}{M}_{j}\left({\rm{\Delta }}x,E,{k}_{y}\right)=\left(\begin{array}{ll}\displaystyle \frac{\cos \,{\theta }_{j}\,\cos \left({q}_{j}{\rm{\Delta }}x\right)+\tfrac{1}{\tau }\,\sin \,{\theta }_{j}\,\sin \left({q}_{j}{\rm{\Delta }}x\right)}{cos{\theta }_{j}} & i\displaystyle \frac{\sin \left({q}_{j}{\rm{\Delta }}x\right)}{{P}_{j}\,\cos \,{\theta }_{j}}\\ i\displaystyle \frac{{P}_{j}\,\sin \left({q}_{j}{\rm{\Delta }}x\right)}{\cos \,{\theta }_{j}} & \displaystyle \frac{\cos \,{\theta }_{j}\,\cos \left({q}_{j}{\rm{\Delta }}x\right)-\tfrac{1}{\tau }\,\sin \,{\theta }_{j}\,\sin \left({q}_{j}{\rm{\Delta }}x\right)}{\cos \,{\theta }_{j}}\end{array}\right).\end{eqnarray} \tag{ 10 }$$

Given the components of the matrix M, it can be shown that determinant of Mj(Δx, E, k_y ) is unique, det[Mj(Δx, E, k_y )] = 1.

2.2. Wavefunction, transmission and reflection coefficient

In order to determine two-component wavefunction at each arbitrary point of x, it is necessary to have information about it at initial point of x₀, $\left(\begin{array}{c}{\psi }_{A}\left({x}_{0}\right)\\ {\psi }_{B}\left({x}_{0}\right)\end{array}\right).$ Having boundary conditions in hand, depicted in figure 1, we assume free electron with energy of E from x₀ < 0 region enters the planar superlattice at the angle of θ₀. In this region (x₀ < 0), the final wavefunction will be a combination of incident and reflected waves so that at the incident point of the superlattice (x = 0) the wavefunction component ψ_A (x0 = 0) is re-written as:

$$\begin{eqnarray}&&{\psi }_{A}\left(0\right)={\psi }_{i}\left(E,{k}_{y}\right)+{\psi }_{r}\left(E,{k}_{y}\right)=\left(1+r\right){\psi }_{i}\left(E,{k}_{y}\right),\end{eqnarray} \tag{ 11 }$$

ψ_i (E, k_y ) is the electron incident wave packet at x = 0, and ψ_r (E,k_y ) = rψ_i (E,k_y ) where r is reflection coefficient. From equations (5a ) and (6a ), it is easy to obtain φ₁(0) = (1+r)ψ_i (E,k_y ) and φ₂(0) = (1−r)ψ_i (E,k_y ). Using equation (8) ψ_B (0) is obtained as:

$$\begin{eqnarray}&&\begin{array}{l}{\psi }_{B}\left(0\right)=\displaystyle \frac{i}{{\tau }_{0}}{P}_{0}\,\sin \,{\theta }_{0}{\varphi }_{1}\left(0\right)+{P}_{0}\,\cos \,{\theta }_{0}{\varphi }_{2}\left(0\right)\\ \,=\,{P}_{0}\left[\left(\cos \,{\theta }_{0}+\displaystyle \frac{i}{{\tau }_{0}}\,\sin \,{\theta }_{0}\right)-r\left(\cos \,{\theta }_{0}-\displaystyle \frac{i}{{\tau }_{0}}\,\sin \,{\theta }_{0}\right)\right]{\psi }_{i}\left(E,{k}_{y}\right).\end{array}\end{eqnarray} \tag{ 12 }$$

And finally:

$$\begin{eqnarray}&&\left(\begin{array}{c}{\psi }_{A}\left(0\right)\\ {\psi }_{B}\left(0\right)\end{array}\right)=\left(\begin{array}{c}1+r\\ {P}_{0}\left[\left(\cos \,{\theta }_{0}+\displaystyle \frac{i}{{\tau }_{0}}\,\sin \,{\theta }_{0}\right)-r\left(\cos ,\,,{\theta }_{0},-,\displaystyle \frac{i}{{\tau }_{0}},\,,\sin ,\,,{\theta }_{0}\right)\right]\end{array}\right){\psi }_{i}\left(E,{k}_{y}\right).\end{eqnarray} \tag{ 13 }$$

By defining t as total transmission coefficient and by using equation (7) in final point of superlattice (x = x_e ) we have:

$$\begin{eqnarray}&&\left(\begin{array}{c}{\psi }_{A}\left({x}_{e}\right)\\ {\psi }_{B}\left({x}_{e}\right)\end{array}\right)=\left(\begin{array}{c}t\\ t{P}_{e}\left(\cos \,{\theta }_{e}+\displaystyle \frac{i}{{\tau }_{e}}\,\sin \,{\theta }_{e}\right)\end{array}\right){\psi }_{i}\left(E,{k}_{y}\right).\end{eqnarray} \tag{ 14 }$$

θ_e is exit angle at the end of the superlattice. X is called transfer matrix that connects wavefunction at incident point to that of the end of superlattice as:

$$\begin{eqnarray}&&\left(\begin{array}{c}{\psi }_{A}\left({x}_{e}\right)\\ {\psi }_{B}\left({x}_{e}\right)\end{array}\right)={\boldsymbol{X}}\left(\begin{array}{c}{\psi }_{A}\left(0\right)\\ {\psi }_{B}\left(0\right)\end{array}\right).\end{eqnarray} \tag{ 15 }$$

Since the whole structure of the superlattice is composed of N ribbons with width of w_j , then one can construct X as multiplying matrices M_j as:

$$\begin{eqnarray}{\bf{X}}=\left(\begin{array}{cc}{X}_{11} & {X}_{12}\\ {X}_{21} & {X}_{22}\end{array}\right)=\displaystyle \prod _{j=1}^{N}{M}_{j}\left({w}_{j},E,{k}_{y}\right).\end{eqnarray} \tag{ 16 }$$

Considering detM_j = 1, it is very easy to reach detX = 1. Therefore, using equations (13)–(15) reflection and transmission coefficient (r, t) for whole structure can be calculated. In particular, for τ_A,B = ±1, r and t are obtained as:

$$\begin{eqnarray}&&r\left(E,{k}_{y}\right)=\displaystyle \frac{{X}_{22}{P}_{0}{e}^{\pm i{\theta }_{0}}-{X}_{11}{P}_{e}{e}^{\pm i{\theta }_{e}}-{X}_{12}{P}_{0}{P}_{e}{e}^{\pm i\left({\theta }_{e}+{\theta }_{0}\right)}+{X}_{21}}{{X}_{22}{P}_{0}{e}^{\mp i{\theta }_{0}}+{X}_{11}{P}_{e}{e}^{\pm i{\theta }_{e}}-{X}_{12}{P}_{0}{P}_{e}{e}^{\pm i\left({\theta }_{e}-{\theta }_{0}\right)}-{X}_{21}},\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&t\left(E,{k}_{y}\right)=\displaystyle \frac{2{P}_{0}cos{\theta }_{0}}{{X}_{22}{P}_{0}{e}^{\mp i{\theta }_{0}}+{X}_{11}{P}_{e}{e}^{\pm i{\theta }_{e}}-{X}_{12}{P}_{0}{P}_{e}{e}^{\pm i\left({\theta }_{e}-{\theta }_{0}\right)}-{X}_{21}}.\end{eqnarray} \tag{ 18 }$$

In order to determine wavefunction ψ_A.B (x) at arbitrary point of x between x_j−1 and x_j (x_j−1 < x < x_j ) in jth layer, we must obtain connection matrix that connects wavefunction ψ_A.B (x₀) at initial point of x₀ to ψ_A.B (x) as:

$$\begin{eqnarray}&&\left(\begin{array}{c}{\psi }_{A}\left(x\right)\\ {\psi }_{B}\left(x\right)\end{array}\right)=Q\left({\rm{\Delta }}{x}_{j},E,{k}_{y}\right)\left(\begin{array}{c}{\psi }_{A}\left({x}_{0}\right)\\ {\psi }_{B}\left({x}_{0}\right)\end{array}\right),\end{eqnarray} \tag{ 19 }$$

where Δx_j = x−x_j−1 and matrix Q, connection matrix can be written as:

$$\begin{eqnarray}&&Q\left({\rm{\Delta }}{x}_{j},E,{k}_{y}\right)={M}_{j}\left({\rm{\Delta }}{x}_{j},E,{k}_{y}\right)\displaystyle \prod _{i=1}^{j-1}{M}_{i}\left({w}_{i},E,{k}_{y}\right).\end{eqnarray} \tag{ 20 }$$

In fact, matrix Q is composed of transfer matrix of (j − 1) ribbons or layers with width of w_i and jth layer with width of Δx_j .

Also, we can determine wavefunction of each sublattices as:

$$\begin{eqnarray}&&{\psi }_{A}\left(x\right)={\psi }_{i}\left(E,{k}_{y}\right)\left[\left(1+r\right){Q}_{11}+{P}_{0}\left[\left(\cos \,{\theta }_{0}+\displaystyle \frac{i}{{\tau }_{0}}\,\sin \,{\theta }_{0}\right)-r\left(\cos \,{\theta }_{0}-\displaystyle \frac{i}{{\tau }_{0}}\,\sin \,{\theta }_{0}\right)\right]{Q}_{12}\right],\end{eqnarray} \tag{ 21a }$$

$$\begin{eqnarray}&&\left\{{\psi }_{B}\left(x\right)={\psi }_{i}\left(E,{k}_{y}\right)\left[\left(1+r\right){Q}_{21}+{P}_{0}\left[\left(\cos \,{\theta }_{0}+\displaystyle \frac{i}{{\tau }_{0}}\,\sin \,{\theta }_{0}\right)-r\left(\cos \,{\theta }_{0}-\displaystyle \frac{i}{{\tau }_{0}}\,\sin \,{\theta }_{0}\right)\right]{Q}_{22}\right],\right.\end{eqnarray} \tag{ 21b }$$

As Q_mn 's are the elements of the matrix Q.

To complete our theoretical approach, we also work on y direction of wavefunction that is obtained by multiplying $Y\left(y,{k}_{y}\right)={e}^{i{k}_{y}{\rm{\Delta }}{y}_{j}}\displaystyle {\prod }_{i=1}^{j-1}{e}^{i{k}_{y}{\rm{\Delta }}{y}_{i}}$ to ${\psi }_{A,B}\left(x\right).$ ${\rm{\Delta }}{y}_{i},$ ${\rm{\Delta }}{y}_{j}$ is defined according to figure 1 as ${\rm{\Delta }}{y}_{i}={w}_{i}\,\tan \,{\theta }_{i}$ and ${\rm{\Delta }}{y}_{j}={\rm{\Delta }}{x}_{j}\,\tan \,{\theta }_{j}.$ So final wavefunction reads:

$$\begin{eqnarray}&&\left\{{\psi }_{A}\left(x,y\right)={\psi }_{i}\left(E,{k}_{y}\right)Y\left(y,{k}_{y}\right)\left[\left(1+r\right){Q}_{11}+{P}_{0}\left[\left(\cos \,{\theta }_{0}+\displaystyle \frac{i}{{\tau }_{0}}\,\sin \,{\theta }_{0}\right)-r\left(\cos \,{\theta }_{0}-\displaystyle \frac{i}{{\tau }_{0}}\,\sin \,{\theta }_{0}\right)\right]{Q}_{12}\right]\right.,\end{eqnarray} \tag{ 22a }$$

$$\begin{eqnarray}&&\left\{{\psi }_{B}\left(x,y\right)={\psi }_{i}\left(E,{k}_{y}\right)Y\left(y,{k}_{y}\right)\left[\left(1+r\right){Q}_{21}+{P}_{0}\left[\left(\cos \,{\theta }_{0}+\displaystyle \frac{i}{{\tau }_{0}}\,\sin \,{\theta }_{0}\right)-r\left(\cos \,{\theta }_{0}-\displaystyle \frac{i}{{\tau }_{0}}\,\sin \,{\theta }_{0}\right)\right]{Q}_{22}\right].\right.\end{eqnarray} \tag{ 22b }$$

2.3. Dispersion relation

Finally, it is straightforward to get dispersion relation for TMDC-PSL very similar to the known approach followed in conventional periodic structure within Bloch theorem. From equation (10) and knowing that $\cos \left({k}_{x}{\rm{\Lambda }}\right)=Tr\left({M}_{A}{M}_{B}\right),$ we obtain dispersion relation as:

$$\begin{eqnarray}&&\cos \left({k}_{x}{\rm{\Lambda }}\right)=\,\cos \left({q}_{A}{w}_{A}\right)\cos \left({q}_{B}{w}_{B}\right)+\left[\displaystyle \frac{1}{cos{\theta }_{A}cos{\theta }_{B}}\left(sin{\theta }_{A}sin{\theta }_{B}-\displaystyle \frac{{P}_{B}}{{P}_{A}}\right)\right]\sin \left({q}_{A}{w}_{A}\right)\sin \left({q}_{B}{w}_{B}\right),\end{eqnarray} \tag{ 23 }$$

where Λ = w_A + w_B .

2.4. Conductance (σ), Fano factor (F)

From equation (18) and transmittance ($T=| t{| }^{2}$ ), conductance and Fano factor reads [51, 52]:

$$\begin{eqnarray}&&G=\displaystyle \frac{4{e}^{2}{L}_{y}}{2\pi h}\displaystyle {\int }_{-\infty }^{+\infty }T\left(E,{k}_{y}\right)d{k}_{y}.\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&F=\displaystyle \frac{\displaystyle {\int }_{-\infty }^{+\infty }T\left(E,{k}_{y}\right)\left[1-T\left(E,{k}_{y}\right)\right]d{k}_{y}}{\displaystyle {\int }_{-\infty }^{+\infty }T\left(E,{k}_{y}\right)d{k}_{y}}.\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&\sigma =G\times \displaystyle \frac{{L}_{x}}{{L}_{y}}.\end{eqnarray} \tag{ 26 }$$

Where all degeneracies are considered and L_x , L_y are superlattice lengths in x, y direction.

In this section, we like to investigate the influence of different parameters on band structure with transport properties in planner TMDC superlattice (TMDC-PSL) by using adopted transfer matrix method presented in previous section. The size of unit cell is Λ = w_A + w_B and w_A = N_A z_A , w_B = N_B z_B are width of A and B ribbons. z_A,B is the size of unit cell in x direction (in-plane growth direction) listed in table 1 and N_A,B is the number of hexagons. We take the length of structure in x direction as L_x = 51Λ + w_A and L_y in y direction (see figure 1). First, we investigate the electron band structure and transmitivity of the structure. Figure 2(a) depicts the band structure of TMDC-PSL containing alternating layers or ribbons of WSe₂ and MoS₂ where w_A = 5z_A and w_B = 5z_B , respectively, with(without) spin-orbit coupling, SOC, effect, shown by blue(red) curves. Energy gap without SOC has been obtained as 1.3 eV. With considering SOC a blue shift in uppest valence band is calculated as Δ = 168 meV where no shift in lowest conduction band is obtained (figure 2(a)). Accordingly, in figure 2(b) transmitivity is illustrated. As it is expected, for energies in gap region transmission coefficient reaches zero and has considerable magnitude for energies close to bands, so that for energies inside band, transmitivity reaches unit.

Zoom In Zoom Out Reset image size

In figures 2(c) and (d) investigations go for valley K'. Compared to figure 2(a) Obvious difference can be seen where uppest valence band is redshifted. The values of energy gaps and uppest valence band shift are also different in valleys' K and K' as indicated in figure [53].

Important feature from comparing figures 2(a) and (c) is attributed to SOC including. Without SOC, the value of energy gap is almost same for K and K' valleys as 1.63 eV. SOC changes energy gap of TMDC-PSL in a few hundreds of meV, that is tunable by means of width of layers (or ribbons) and types of TMDCs.

In figure 3 energy gap pattern versus the ribbons width ratio (or hexagon number ratio, N_A :N_B is illustrated. It is clear that by taking N_A = 0 the value of energy gap is that of WSe₂ and N_B = 0 gives the energy gap of MoS₂. By changing the ratio of N_B and N_A the value of E_g of superlattice is varied in both K and K' valleys, with SOC. Without SOC E_g has same value for K and K' as mentioned and discussed in figure 2. Taking SOC into account results in considerable difference in E_g value in both valleys. In fact, SOC removes degeneracy in the band structure of TMDC-PSL differently in K and K' valleys as it does in the case of TMDC monolayers. In order to get more insight into the effect of width ratio between neighbor ribbons or strips A, B (MoS₂, WSe₂) in unit cell, the value of E_g versus N_A :N_B is shown in figure 3(d). Opposite trends can be seen for changing the value of E_g in K and K' valleys with considering SOC, so that in the case of K valley, E_g increases with increasing the contribution of MoS₂ unit cell numbers against to that of K' valley. In order to compare with other works, we report the result of reference [47] that author's calculated a similar structure of TMDC-PSL within DFT approach. Reasonable consistency can be observed, in particular for $\tfrac{{N}_{A}}{{N}_{B}}=\tfrac{7}{3}$ [N_A :N_B ] = [7:3] consistency is excellent. It is worth to note that Eg is varied with the width of each ribbon because the contribution of materials (MoS2 or WSe2) is different in unit cell when we change the number of hexagons, NA and NB. Besides, SOC in each TMDC material has different magnitude. Therefore, we expect that changing the contribution of material in unit cell by changing NA and NB ratio leads to the different SOC effect in PSL.

Zoom In Zoom Out Reset image size

Figure 4 shows 3D view and counterplot of band structure for certain energy values in energy bands of figure 2(a), recognized as A, B, C, D points in K and K' valleys.

Zoom In Zoom Out Reset image size

Hereafter, we focus on K valley with SOC. In figures 5(a) and (b) we present band structure in k_y direction for different K_x 's that is a cross section of counterplots in figure 4. Parabolic band structure in k_y direction per k_x (that is quantized due to spatial confinement effect in x direction) reflects free carrier behavior in y direction in which there is no confinement effects. We take unit cell length as Λ = 8z_A + 2z_B and Λ = 3z_A + 7z_B . It is clear that by varying ribbons width ratio the energy gap magnitude is increased by increasing N_B . considerable shift to lower energies occurred in valence band. Also main energy gap is widened with increasing N_B . It is worth to mention that unit cell length Λ is equal to w_A + w_B and varying w_A or w_B portion by changing N_A or N_B leads to different structure of superlattice with different ribbon widths. In order to show the energy gap variation with k_y , we take into account the typical values of k_y as zero and 0.5 nm⁻¹. We found valence band red shift for higher k_y . in figures 5(d), (e), 3(e), (d)) we take N_A = N_B = N for k_y = 0.5 nm⁻¹ and a few tens of meV increase in energy gap is observed.

Zoom In Zoom Out Reset image size

Now, we turn to consider the transmission probability in TMDC-PSL. In figure 6(a) transmission probability map versus incident electron energy and transverse wavevector (k_y ) is depicted. It is clear that transmitivity is zero in energy band gap region. At allowed energy minibands transmission probability has nonzero value and in some minibands it reaches unit. Figure 6(b) illustrates the effect of λ, fraction of MoS₂ nanoribbon width to whole unit cell _Λ, on transmitivity. It can be seen that for bigger λ (very wide MoS₂ nanoribbons) transmitivity shows oscillatory trends while for smaller λ it has smooth behavior showing that wider MoS₂ nanoribbon gives rise to resonant states in particular for incoming electron energies lower that −1 eV. In the other words, oscillatory trend of transmitivity indicates quantum nature of transmission through the structure. If the MoS₂ nanoribbon becomes wider, number of resonance states generated for energies lower than −1 eV is increased. Similar behavior is observed for energies higher than about 1 eV. In figure 6(e) the influence of incoming electron incident angle (see geometry in figure 1) is depicted that shows incident angle has no obvious effect on transmitivity. Figure 6(c) is plotted for ${k}_{y}=0.25\left(\tfrac{2\pi }{{\rm{\Lambda }}}\right),$ however our study implies that for other values of k_y incident angle has no effect on transmission. Figure 6(d) demonstrates incremental effect of N, (N = N_A = N_B ), on oscillatory behavior of transmission at all allowed energies. It might be expected that larger number of N that subsequently increases the number of interfaces creates more maxima in transmission spectrum as it is well known in introductory quantum mechanics. Figures 6(e), (f), (d) illustrate a closer view to transmission probability. From these figures it can be understood that effective parameters of N, k_y , λ, and _θ0 how enhance or weaken the transmission probability. Furthermore, it can be seen that oscillatory behavior of T that indicates the presence of resonance energies can be engineered by these parameters.

Zoom In Zoom Out Reset image size

Now we illustrate the evolution of ψ_A and ψ_B , sublattice wavefunctions(or probability amplitude), in figure 7. We take ${\rm{\Lambda }}=5{z}_{A}+5{z}_{B}=3.2515\,nm,$ ${L}_{x}=51{\rm{\Lambda }}+{w}_{A}=162.91\,nm$ as whole system length, _θ0 = 0 and also we choose three energy values: one inside band gap (second column) and two in allowed neighbor higher band (third column) and lower (first column) band. Panels in each row also depicts probability amplitudes (ψ_A , ψ_B ) for real and imaginary parts of wavefunction in sublattices A, B with probability densities $| {\psi }_{A}{| }^{2},$ $| {\psi }_{B}{| }^{2}.$ From first column in figure 7, for E = +0.825 eV in allowed band higher than band gap (reminiscent of conduction band in two-band model of a typical semiconductor) wavefunction behaves as propagating wave or Bloch wave in entire system. From dispersion relation, equation (23), for E = +0.825 eV, kx = 0.3287 nm⁻¹ then attributed wavelength is obtained as λ_B = 19.1153nm, so it is expected 8.52 periods $\left(\tfrac{{{\rm{L}}}_{x}}{{\lambda }_{B}}\right)$ are governed by propagating wave that is confirmed by figures 7(a), (b) (determined by the number of maxima or minima). Second column indicates wavefunctions (probability amplitude) and probability density for an energy inside band gap that is decaying wave which is caused by this fact that there is no electronic state in the band. Panels g-i show mentioned parameters behavior for an energy in bands immediately lower than band gap as E = −0.648 eV that shows again a propagating wave feature similar to panels a–c.

Zoom In Zoom Out Reset image size

According to reports on fabricating planar TMDC superlattice [3, 4], ribbons have inherent disorders in width, hence, in order to study electronic and transport properties more precisely, disorder effects are considered. To model these effects, we take whole system length as ${L}_{x}=51{\rm{\Lambda }}+{w}_{A}$ and disorder length as three values of ${R}_{1},$ ${R}_{2},$ ${R}_{3}$ $\left({w}_{A\left(B\right)}=5{z}_{A\left(B\right)}+{R}_{1,2,3}\right)$ which are random numbers in three intervals of $\tfrac{1}{2}\left(-{z}_{A,B},{z}_{A,B}\right),$ $\left(-{z}_{A,B},{z}_{A,B}\right),$ $\tfrac{3}{2}\left(-{z}_{A,B},{z}_{A,B}\right),$ respectively for ${\theta }_{0}=0.$ Figure 8 shows the electronic transmittivity through the structure. It is clear that main gap (around $K=0$) is independent of different disorders, but, higher and lower neighbor band gaps (that are determined by minima at two sides of main minimum ) are sensitive to disorder as shown in figures 8(b) and (c). Disorder effect in upper band gap is more remarkable. In figure 9 electronic transmittivity spectrum under different incident angle is depicted. It is obvious that main gap is very weakly dependent on incident angle but two other gaps experience shift with increasing incident angle. Therefore, main gap is almost independent of incident angle and disorder.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this part we investigate the effect of defect on transmission by considering defected ribbon in TMDC-PSL as ${(AB)}^{25}D{(BA)}^{25}.$ "D" denotes the defected ribbon, which is a wider ribbon of $WS{e}_{2}.$ In figure 10 transmitivity for 3 different widths of defect is illustrated for incident angle ${\theta }_{0}=0.$ Similar to disorder effect, it can be seen that main gap is not influenced by including defect but two aside gaps are shiftted. Also defect modes are appeared inside gaps determined by sharp peaks. In addition, shift in upper gap becomes remarkable with increasing defect ribbon width (see figures 10(b), (c)). Figure 10(d) demonstrates the behavior of all defect modes inside energy gaps by changing the width of defect ribbon from ${w}_{D}=1{z}_{A}$ to ${w}_{D}=15{z}_{A}.$

Zoom In Zoom Out Reset image size

It is obvious that all defect modes below zero energy experience blue shift and those of higher than zero energy are red shifted. From figure, when width of defect ribbon is $5{z}_{A},$ that is the width of ribbon A in superlattice, defect modes disappear because system is normal ${(AB)}^{25}A{(BA)}^{25}$ superlattice. Figure 10(e) presents electronic transmitivity within structure under different incident angles. Obviously, for particular width of defect ribbon the location of defect modes in all energy gaps are independent of the angles, see figures 10(f), (g), (h) that show this finding in closer view.

Figure 11 illustrates the evolution of both the probability densities and amplitudes of ${\psi }_{A,B}(x)$ inside TMDC- PSL with a defect as ${(AB)}^{25}D{(BA)}^{25}$ at defect mode energies (points A and B in figure 11(a)). In figures 11(b), (c), (d) we take ${E}_{A}$ (energy at point A) that is −1.22 eV and incident angle ${\theta }_{0}=0.$ It is clear that carrier's probability density is highly localized inside defect ribbon. To complete description of wavefunction evolution at energy ${E}_{A},$ it is interesting to see that incident wave is propagating at the entrance of structure (x = 0) and starts to decay because ${E}_{A}$ is inside the lowest gap (figure 11(a), gap A), and localize inside defect and continues decaying. Similar behavior is seen for ${\psi }_{B}(x)$ (figure 11(b)). Figure 11(d) shows localization more clearly in which probability density reaches peak in defect ribbon. Figures 11(e), (f), (g) depicts similar trends for defect mode of B in the above band gap in figure 11(a), and again localization feature is pronounced in defect region.

Zoom In Zoom Out Reset image size

3.1. Conductance (σ), Fano factor (F)

Now, we consider the transmission of an electron moving inside TMDC-PSL e. g. ${(AB)}^{25}A{(BA)}^{25}$ for different layers width ratio, $\lambda ,$ under special incident angle of ${\theta }_{0}=0$ and we focus on Fano factor that is very important parameter in addressing shot noise and conductivity in transport properties of nanostructures. Figure 12 illustrates transmission, conductivity and Fano factor versus energy with different $\lambda$ for finite periodic structure. Varying $\lambda$ has considerable effect on transmission probability that is depicted in figure 12(a). Two typical values of $\lambda$ have been considered that change the transmission probability spectrum by shifting it. The amount of shift is much more pronounced for lower subband (lower than main gap). Similar trend can be seen in conductivity spectrum, see figure 12(b). As it is expected, conductivity reaches zero in gap region for both values of $\lambda ,$ again shift of conductivity spectrum is more considerable in energies below zero. Corresponding Fano factor in the band gap tends to integer 1, figure 12(c). The effect of changing $\lambda$ has still significant effect on lower subband energies. Outside the band gap, as the energy increases (absolute value) conductivity increases linearly with oscillations caused by propagating modes in finite structure. Also we notice that the Fano factor changes from 1 to smaller values gradually, indicating that transport becomes ballistic. In order to get deeper insight, figures 12(d), (e), (f) depict transmitivity, conductivity and Fano factor versus $\lambda$ and certain values of energies in upper and lower subbands. Apparently for pure $Mo{S}_{2}$ (B ribbon) that is $\lambda =1,$ conductivity exhibits a factor of minimum $4{e}^{2}/\pi h$ at energies higher and lower than gap region. Furthermore, for $\lambda =1$ Fano factor has minimum value for these energies, even below $1/3.$

Zoom In Zoom Out Reset image size

In figure 13 the effect of incident angle is investigated. To give an overview in figures 13(a), (b), (c), transmission probability, conductivity and Fano factor spectrum for two typical values of incident angle, are shown. The main gap and other lower gaps are independent from incident angle and transmitivity variation with changing incident angle is not remarkable. From figures 13(b), (c), incident angle doesn't shift the conductivity and Fano factor spectrum but changes their values outside gap region. Furthermore, in figure 13(c) at energy around −1.25 eV Fano factor is above $1/3$ for ${\theta }_{0}={30}^{\circ }$ and less than $1/3$ for ${\theta }_{0}={5}^{\circ }$ showing that Fano factor can be adjusted by incident angle. In figures 13(d), (e), (f) for particular energies of lower and upper subbands (than main bandgap), incident angle is swept. Interestingly, it is found that incident angle has no considerable impact on transmission probability but changes conductivity and Fano factor more remarkably, which has decremental and incremental effect, respectively, for all values of energies (not shown here).

Zoom In Zoom Out Reset image size

Finally, we turn our attention to the effect of unit cell size, ${\rm{\Lambda }},$ on transport properties of TMDC-PSL. As it was mentioned before, increasing unit cell size results in narrowing main bandgap. Hence, from figure 14(a), it is obvious that transmitivity and main energy gap are varied remarkably . Similarly, considerable effects can be seen in conductivity and Fano factor spectrum. By decreasing unit cell size, conductivity is weakened, and subsequently Fano factor is enhanced. Also, decreasing unit cell leads to shifting conductivity and Fano factor spectrum and increasing the number of picks due to enhancing propagation modes of superlattice. In figures 14(d), (e), (f), the effect of varying unit cell size has been considered to show its overall effect, for typical energies in upper and lower energy subbands. Opposite behavior of conductivity and Fano factor is still seen.

Zoom In Zoom Out Reset image size

We found that without spin orbit coupling, the value of energy gap is almost same for $K$ and ${K}^{{\prime} }$ valleys . SOC changes energy gap of TMDC-PSL in a few hundreds of meV, that is tunable by the width of ribbons and kind of materials. Furthermore, SOC removes degeneracy in the band structure of TMDC-PSL and shifts valence bands differently in k and k' valleys. Opposite trends were seen for behavior of ${E}_{g}$ in $K$ and ${K}^{{\prime} }$ valleys with considering SOC. For $K$ valley, ${E}_{g}$ increases with increasing the portion of $Mo{S}_{2}$ in unit cell against to that occurred in ${K}^{{\prime} }$ valley. Also, valence band experiences red shift for higher value of ${k}_{y}.$ It is also shown that for bigger $\lambda$ (wide $Mo{S}_{2}$ nanoribbons) transmitivity shows oscillatory trends while for smaller $\lambda$ it has smooth behavior showing that wider $Mo{S}_{2}$ nanoribbon gives rise to resonant states. Parameters $N,$ ${k}_{y},$ $\lambda ,$ and ${\theta }_{0}$ enhance or weaken the transmission probability. We further found that main gap is almost independent of incident angle and disorder and two neighbor gaps experience shift with defect modes that are appeared inside these gaps. Finally, we have seen that for pure MoS₂ ($\lambda =1$) conductivity exhibits a factor of minimum $4{e}^{2}/\pi h$ at energies higher and lower than gap region. Furthermore, for $\lambda =1$ Fano factor has minimum value for these energies, even below $1/3.$ Fano factor above $1/3$ is seen for certain value of energies and incident angle. Decreasing unit cell (number of hexagones) leads to conductivity and Fano factor spectrum shift with increasing the number of peaks due to enhancing propagation modes of superlattice. While our theoretical study is simple compared to DFT approaches, presents an efficient way to control the key parameters effects on conductivity and band structure of TMDC-PSL and paves a way towards pertinent experiments.

The data generated and/or analysed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request.