Green’s function approach to edge states in transition metal dichalcogenides
Mojtaba Farmanbar, Taher Amlaki, and Geert Brocks
Faculty of Science and Technology and MESA+ Institute for Nanotechnology, University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands
(Received 21 March 2016; revised manuscript received 4 May 2016; published 31 May 2016) The semiconducting two-dimensional transition metal dichalcogenides MX2 show an abundance of one-dimensional metallic edges and grain boundaries. Standard techniques for calculating edge states typically model nanoribbons, and require the use of supercells. In this paper, we formulate a Green’s function technique for calculating edge states of (semi-)infinite two-dimensional systems with a single well-defined edge or grain boundary. We express Green’s functions in terms of Bloch matrices, constructed from the solutions of a quadratic eigenvalue equation. The technique can be applied to any localized basis representation of the Hamiltonian. Here, we use it to calculate edge states of MX2monolayers by means of tight-binding models. Aside from the basic zigzag and armchair edges, we study edges with a more general orientation, structurally modifed edges, and grain boundaries. A simple three-band model captures an important part of the edge electronic structures. An 11-band model comprising all valence orbitals of the M and X atoms is required to obtain all edge states with energies in the MX2band gap. Here, states of odd symmetry with respect to a mirror plane through the layer of M atoms have a dangling-bond character, and tend to pin the Fermi level.
DOI:10.1103/PhysRevB.93.205444 I. INTRODUCTION
In the wake of graphene, a whole new class has emerged of materials that are essentially two dimensional (2D) in nature [1,2]. The subset of materials with honeycomblike structures alone contains metals such as graphene, insulators such as boron nitride (h-BN), and semiconductors such as the transition metal dichalcogenides (MX2; M= Mo,W, X =
S,Se,Te) [3,4]. It becomes more and more feasible to grow nanostructures and in-plane heterostructures of 2D materials in a controlled way [5–11]. As a result, the electronic structure of edges and grain boundaries attracts increasing attention.
Graphene edges, for instance, are predicted to have re-markable one-dimensional electronic and magnetic properties [12,13]. The edges and grain boundaries of MX2 sheets are
generally metallic [10,14–20]. As the bulk 2D MX2materials
are semiconducting, the metallicity is truly localized at the edge or grain boundary, where one could see manifestations of the peculiar spectral and transport properties of one-dimensional (1D) metals [21]. In chemistry, MoS2edges have
been identified as sites that show a special catalytic activity [14,16,22]. Such experimental developments have motivated a large number of calculations on the electronic structure of edge states, both first-principles calculations, as well as model calculations.
The electronic structure of graphene edges in particular have been studied extensively. As graphene has a relatively simple electronic structure, some features of the edge states in graphene can be studied by analytical or simple numerical techniques [23–25]. The edges of MoS2have attracted renewed
attention recently [14–20,26–33]. Here, the complexity of the electronic structure requires more extensive models, even for relatively simple modeling at the tight-binding level [34,35], or at the continuum level [36,37]. A standard technique for calculating edge states uses supercells to model nanoribbons of a finite width. Drawbacks of this approach are that the electronic structure of the edge states is mixed with that of the bulklike interior of the nanoribbon, and that the two edges
of the nanoribbon can interact electronically. The ribbon has therefore to be sufficiently wide in order to electronically separate the two edges from one another and from the bulk. It may require the use of large supercells, which in particular if the materials are modeled from first principles leads to time-consuming calculations and complicates analysis of the results.
Green’s function techniques constitute an alternative ap-proach. They enable calculations on semi-infinite structures with a single well-defined edge, or on infinite structures containing a single grain boundary, and they generally do not require the use of large supercells. Green’s function techniques have been pioneered for calculations on surface states of three-dimensional (3D) materials [38–42]. Here, we formulate a special Green’s function technique for calculating edge and grain boundary states. This technique is inspired by the Green’s function formalism that has been introduced for calculating electronic transport through nanostructures [43–47].
We express both edge Green’s functions and bulk Green’s functions in terms of Bloch matrices, which are constructed from the eigenvalues and eigenfunctions of a quadratic eigen-value problem [43–48]. Structural and chemical modifications at the edges and grain boundaries are then tackled by connecting the modified edges to semi-infinite bulk structures. The method allows for a clean separation of edge and bulk properties at a moderate computational cost.
Our Green’s function approach requires a representation of the Hamiltonian on a localized basis, such as atomic orbitals or Wannier functions, or a real-space grid. It can be applied to tight binding as well as first-principles models. In this paper, we illustrate its use on tight-binding models for MX2
monolayers, MoS2 in particular. We apply the approach to
the three-band model for the simplified electronic structure of MoS2, developed by Mattheis [49] and Liu et al. [34]. We study
the electronic structures of the elementary (zigzag, armchair) MoS2 edges, and of edges of a more general orientation.
We illustrate the effect on the electronic structure of edge modifications and of the formation of grain boundaries. For
comparison, we also study edges within an 11-band model comprising all MoS2valence orbitals, developed by Cappelluti et al. [50].
The paper is organized as follows. In Sec.II A, we formulate the technique of calculating Green’s functions using Bloch matrices. How to apply this technique to MX2edges and grain
boundaries is described in Sec.II B. The tight-binding models are discussed in Sec.II Cand AppendixesAandB. We discuss results obtained with these models for MX2edges and grain
boundaries in Secs.III AandIII B.
II. THEORY A. Green’s functions
We divide the 2D layers into 1D strips (see Fig.1) [40]. Assuming translational symmetry with period a along the strip, the Hamiltonian matrix can be labeled by a Bloch wave number −π
a < k π
a. For clarity of notation, we suppress the label k in the following. The thickness of the strips is chosen such that a direct interaction exists between neighboring strips only. We label the strips by an index i, and divide the Hamiltonian matrix into blocks Hi,j, with an i,j = −∞, . . . ,∞ for an infinite system. Having only nearest-neighbor interactions between strips means that the Hamiltonian matrix is block tridiagonal, i.e., Hi,j = 0; j = i,i ± 1. For a unit cell in the strip containing N orbitals, all these matrix blocks are N× N. The columns of the retarded Green’s function matrix blocks
Gi,j obey
− Hi,i−1Gi−1,j+ (E+I− Hi,i)Gi,j − Hi,i+1Gi+1,j = Iδij, (1) where I is the N× N identity matrix, and E+= E + iη with
η the usual infinitesimal positive real number. Note that we assume a representation based upon an orthogonal basis set. We are foremost interested in the layer-resolved density of states, given by the usual expression ni = −π−1ImTrGi,i. Obviously, besides on the Bloch wave number k, the Green’s function matrix also depends on the energy E. Again, for ease of notation we often omit both these labels in the following.
FIG. 1. A 2D layer is divided in 1D strips. Translation symmetry along the strips gives the Bloch wave number k. Translation symmetry between the strips results in identical on-strip Hamiltonian matrix blocks H(k) and hopping matrix blocks B(k).
In an infinite system with translational symmetry between the layers, the strips are identical, and Eq. (1) becomes
− BGi−1,j + (E+I− H)Gi,j − B†Gi+1,j = Iδij, (2) with B= Hi,i−1, H= Hi,i, and B†= Hi,i+1 for i= −∞, . . . ,∞ (see Fig.1). For a left semi-infinite system with
i= −∞ to i = 0, the equations remain the same, except for
the i= 0 equation, which becomes − BGL
−1,j+ (E+I− H)GL0,j = Iδ0j. (3)
In surface science, the matrix block gL= GL0,0 is called the surface Green’s function. In the context of 2D structures, we call it the edge Green’s function.
In a similar way, for a right semi-infinite system with i= 0, . . . ,∞, only the equation for i = 0 is different from the bulk equation (2),
(E+I− H)GR0,j− B†GR1,j = Iδ0j, (4)
with the matrix block gR = GR0,0 the edge Green’s function. Note that if the 2D layer has no twofold symmetry, such as inversion, a mirror plane perpendicular to the 2D layer, or a twofold rotation axis perpendicular to the 2D layer, then
gR = gL, i.e., a right edge is different from a left edge. The density of states of the edge strip of a left/right semi-infinite system is given by
n(S)L/R =
N α=1
nα; nα= −π−1Im[gL/R]α,α, (5) where nα is the density of states projected on a single orbital
α.
Partitioning the infinite system into a right and a left semi-infinite half, the on-strip matrix blocks of the Green’s function of the infinite system can be expressed as
Gi,i =
g−1L − B†gRB −1
. (6)
1. Eigenmodes and Bloch matrices
We will express the Green’s function matrices of (semi-)infinite systems in terms of Bloch matrices [43]. Equation (2) for i = j is the same as the Schr¨odinger equation in tight-binding representation
− Bci−1+ (EI − H)ci− B†ci+1= 0, (7) where ci is the N -dimensional vector of orbital coefficients of the wave function on strip i. With translational symmetry between the strips, the elementary solution is a Bloch wave, and ci+1= λci, with λ a complex constant. For an infinite system this holds for all i, and for a left (right) semi-infinite system for i < 0 (i > 0). Using this relation as an ansatz in Eq. (7) gives a quadratic eigenvalue equation in λ of dimension
N:
A(λ)u= 0; A(λ) = −B + λ(EI − H) − λ2B†. (8) As A(λ)= λ2A†(1/λ∗), if λ is a root of det A(λ), then so
is 1/λ∗. Numerically, the solutions are usually found by solving an equivalent linear generalized eigenvalue equation
of dimension 2N [48]: I 0 EI− H −B − λ 0 I B† 0 z= 0, (9) which resembles the eigenvalue equation for the transfer matrix [40,41,51–54].
The maximally 2N solutions of this equation can be divided into two classes, i.e., right-going modes and left-going modes, labeled, respectively, by+ and − superscripts in the following. Right-going modes are either evanescent waves that decay to the right|λ+n| < 1 or waves that travel to the right, meaning |λ+n| = 1 and a positive group velocity. Left-going modes either decay to the left |λ−n| > 1 or travel to the left |λ−
n| = 1 with negative group velocity. With unthe eigenvector belonging to the eigenvalue λnin Eq. (9), the group velocity is given by [45]
vn= − 2a
Im[λnu†nB†un]. (10) Only for traveling waves is the group velocity nonzero.
We divide the eigenvectors into a set of N+ right-going modes u+nand a set of N−left-going modes u−n. The evanescent waves always come in pairs of a right-going and a left-going mode, i.e., if λ+n is the eigenvalue of a right-going mode, then
λ−n = 1/(λ+n)∗gives a left-going solution. Traveling waves do not necessarily come in such pairs, and the numbers of and left-going traveling waves may be different. Neither right-going modes nor left-right-going modes necessarily form a complete set in N -dimensional space, nor are they an orthogonal set in general.
One can use these two sets of modes to form the two N×
N±matrices
U±= (u±1,u±1, . . . ,u±N±), (11)
and construct the two Bloch matrices
F±= U±± U±, (12) where ± are the N±× N± diagonal matrices with the eigenvalues λ±n on the diagonal, and U± are the N±× N (pseudo)inverses of U± [55]. The Bloch matrices have the convenient property
(F±)p= U±(±)p U± (13) for any integer p, as U±U±= IM±, the M±× M± identity matrix, with M± N±. It follows that
ci = (F+)ic0++ (F−)ic−0 (14) satisfies the tight-binding equation (7), where c±0 set the boundary conditions in strip number 0. We assume that all relevant solutions can be expressed this way.
2. Green’s functions in terms of Bloch matrices
The general expression of Eq. (14) also applies to the columns of the Green’s function matrices [compare Eqs. (2) and (7)]. The boundary conditions require that a retarded Green’s function comprises traveling waves moving outwards from its point source and/or evanescent waves decaying away from the source. For left and right semi-infinite systems, this
gives
GLi,0 = (F−)
i
gL, i <0 (15)
GRi,0= (F+)igR, i >0. (16) Using Eq. (15) in Eqs. (2) and (3) then leads to
[−B(F−)−1+ E+I− H]g
L= I,
[−B(F−)−1+ E+I− H − B†F−](F−)−1gL= 0. (17) Solving these two equations gives the edge Green’s function of a left semi-infinite system as
gL = [B†F−]−1, (18) where the inversion should be treated as a pseudoinversion if
B†is singular [55]. Using Eq. (16) in Eqs. (2) and (4) gives the edge Green’s function of a right semi-infinite system
gR = [B(F+)−1]−1. (19)
Finally, using Eq. (6) gives the on-strip Green’s function matrix block of an infinite system
Gi,i = [B†F−− B†F+]−1. (20)
3. Ideal edge states
One cannot have traveling Bloch waves for energies in the band gap of a semiconductor. In semi-infinite systems one can, however, have solutions in the form of evanescent states that originate from the edge of the system. One can find the energies of these edge states from the edge Green’s functions [Eqs. (18) and (19)], which have isolated poles at these energies. Obviously in numerical calculations one has to work at complex energies E+ iη to avoid these poles, but η can be chosen small.
Alternatively, edge states can be obtained from the solutions of the eigenvalue problem [Eq. (8)], solved at real energies E [40,53,56,57]. As edge states should decay away from the edge, only the u−n modes can contribute to an edge state for a left semi-infinite system, and only the u+n modes for a right semi-infinite system. An edge state of a left semi-infinite system has amplitude zero beyond the edge of that system, i.e., c1= F−c−0 = 0 [Eq. (14)]. This means rank(F−) < N−,
the number of left-going solutions of Eq. (8). Because rank(F−)= rank(U−) [see Eq. (12) with rank(−)= N−and rank(U−)= rank( U−)], a necessary and sufficient condition for the existence of an edge state is that the eigenmodes u−n are linearly dependent. A similar reasoning holds for the edge states of a right semi-infinite system. The number of edge states at a particular energy E and wave number k of a left/right semi-infinite system is then given by NedgeL/R= N−/+− rank(U−/+).
4. Ideal grain boundaries
A model for an ideal grain boundary is shown in Fig.2. Space is divided into two parts with B the hopping matrix block connecting the left and right halves. We assume that the on-strip Hamiltonian matrix blocks of all the strips in the left half are given by HL right up to the boundary, and that the hopping matrix blocks between all nearest-neighbor strips in the left part are given by BL. The corresponding matrix blocks
FIG. 2. Grain boundary between a left semi-infinite system with on-strip Hamiltonian matrix blocks HL(k) and hopping matrix blocks
BL(k) and a right semi-infinite system with matrix blocks HR(k) and
BR(k). The coupling between left and right parts is given by the
hopping matrix block B(k).
for the right half are HR and BR, respectively. The Green’s function matrix blocks gI
Land g I
R pertaining to the two strips just left and right of the grain boundary interface can be derived like Eq. (6): gIL= g−1L − B†gRB −1 , (21) gIR = g−1R − BgLB† −1 , (22)
where gLand gR are the edge Green functions of the left and right semi-infinite systems, respectively. With Eqs. (18) and (19), one can express the interface Green’s functions in terms of Bloch matrices
gIL = [B†LF−L− B(BR(F+R)−1)−1B†]−1, (23)
gIR = [BR(F+R)−1− B(BLF−L)−1B†]−1. (24)
5. Modified edges
So far, we have assumed that all layers are identical right up to an edge or grain boundary. The formation of an edge or boundary often involves an electronic and a structural reconstruction, which makes the Hamiltonian matrix blocks of the strips adjacent to an edge or grain boundary different from those of the bulk strips.
We illustrate this on a left semi-infinite system with an edge layer (i= 0) that is different from the bulk. The tight-binding equations for the two layers closest to the edge are
− Bc−2+ (EI − H)c−1− B†c0= 0,
−Bc
−1+ (EI − H)c0= 0, (25)
where H= H is the onsite Hamiltonian of the edge layer, and
B= B is the coupling of this layer to the rest of the system
(see Fig.3). Using c−2= (F−)−1c−1[see Eqs. (14) and (15)] transforms Eq. (25) into
EI− H − B(F−)−1 −B† −B EI− H c−1 c0 = 0 0 . (26)
FIG. 3. Left semi-infinite system with a modified edge. H(k) and B(k) are the matrix blocks of the unperturbed system; H(k) and B(k) are the matrix blocks of the perturbed edge strip.
In the terminology used in nonequilibrium Green’s function (NEGF) transport calculations, the termL(E)≡ B(F−)−1is called the self-energy of the left lead [45]. The Green’s function matrix blocks pertaining to the two layers closest to the edge are then given by
GL= EI− H − B(F−)−1 −B† −B EI− H −1 . (27) The Green’s function gives the density of states in the usual way [cf. Eq. (5)]. The density of states is zero for energies inside the band gap, except for isolated poles at particular energies, which represent the edge states of the modified edges. This formalism can be adapted in an obvious way to model modified edges of right semi-infinite systems, or grain boundaries where the strips left and right of the interface are modified.
6. Charge-neutrality level
Experimentally, the Fermi level in MX2compounds is often
determined by unintentional doping due to defects [58–60]. In addition, as MX2compounds are polar materials, an internal
electric field is created if a crystallite is terminated by polar edges [61,62]. Such an electric field can cause a long-range charge transfer between different edges, even if the bulk material does not contain any impurities. The position of the Fermi level can then become dependent on the size and the shape of the sample, which makes the intrinsic Fermi level an ill-defined quantity.
Each edge or grain boundary has a well-defined energy level at which that edge or grain boundary is electrically neutral, the charge-neutrality level (CNL). We use the CNL as a reference point in the following. The charge-neutrality level ECNL is
defined as the energy at which
N(ECNL)= Nstrip, (28)
where Nstrip is the number of electrons that makes the strip
neutral, and N (E) is the electron counting function
N(E)= E
−∞
FIG. 4. Top view of the MX2 lattice with M and X atoms in purple and yellow, respectively. The lattice vectors a1 and a2 are indicated by arrows. The M and X zigzag edges run parallel to a2 and are the edges of right and left semi-infinite systems, respectively. The armchair edge runs parallel to 2a1+ a2.
with n(E) the k-integrated density of states
n(E)= a 2π π a −π a n(E,k)dk. (30)
The k-resolved density of states n(E,k) can be obtained from Eq. (5) for edge strips, and a similar expression for a bulk strip.
B. M X2edges
The hexagonal lattice of MX2is shown in Fig.4. We specify
an edge starting from a supercell spanned by vectors T1 and T2. This supercell is used to define a semi-infinite system,
choosing one of the vectors as the translation vector parallel to the edge.
1. Zigzag and armchair edges
Similar to graphene, the basic-type edges of the MX2lattice
are the zigzag and armchair edges as defined in Fig.4. A zigzag edge is defined by T1= a1and T2= a2, with T2as the vector
parallel to the edge. The matrix blocks discussed in Sec.II A
become
H= H0,0+ H0,1ei2π k+ H0,−1e−i2πk, (31) B= H−1,0+ H−1,−1e−i2πk, (32) where Hp,q denotes the real-space Hamiltonian matrix block that describes the interaction between atoms in the unit cell situated at the origin and atoms in the unit cell situated at
pa1+ qa2. The matrix elements of Hp,qdepend on the specific tight-binding model that is used to represent the electronic structure of MX2. They are given in the Appendixes.
Note that by solving the quadratic eigenvalue equations (8) and (9), one has simultaneous access to the Green’s functions of the edges of a right and a left semi-infinite system via Eqs. (18) and (19). Unlike graphene, the MX2 monolayer
lacks inversion symmetry, which means that the zigzag edge termination of a right semi-infinite system is different from that of a left semi-infinite system (see Fig.4). The zigzag edge of a right semi-infinite system is terminated by metal atoms.
We call this the M edge, consistent with previous studies on MoS2, where it is called the Mo edge [14–16]. The zigzag
edge of a left semi-infinite system is then called the X edge, as it is terminated by chalcogen atoms. In previous studies on MoS2, it has been called the S edge.
An edge in armchair orientation can be constructed from
T1= 2a1+ a2 and T2= a2 with T1 the translation vector
parallel to the edge (see Fig.4). The corresponding supercell is twice as large as the unit cell, so the dimension of the matrix blocks defining the edge are twice the dimension of the blocks defining the zigzag edge
H= H0,0 H1,1+ H−1,0e−i2πk H1,0ei2π k+ H−1,−1 H0,0 , (33) B= H0,−1 H1,0+ H−1,−1e−i2πk 0 H0,−1 . (34)
Note that the termination at the edge is controlled by the contents of the cell used to define the primitive vectors a1,a2.
In particular, the cell defined in Fig.4does not lead to a pristine armchair edge, but one with additional M or X atoms attached to the edge. It is straightforward to remove these atoms and by applying the technique outlined in Sec.II A 5 obtain the electronic structure of a pristine armchair edge.
2. General edges
Edges with a somewhat more general orientation are defined by the supercell
T1 = m(a1+ a2)+ n(2a1+ a2); T2= a2, (35)
with T1the translation vector parallel to the edge, defined as m
zigzag vectors plus n armchair vectors (see Fig.5). The angle with the direction of a2is given by
θ = arccos 1 2m √ m2+ 3mn + 3n2 . (36)
Because of the symmetry of the lattice, one only has to cover the 30◦angle between a zigzag orientation m= 1, n = 0, θ = 60◦, and an armchair orientation m= 0, n = 1, θ = 90◦. Left and right edges (M-type and X-type edges) are then obtained by using T2= a2and T2= −a2, respectively.
A series of edge orientations is obtained by setting n= 1 and varying m, where the translation vector along the edge
FIG. 5. Left: generalized zigzag edge parallel to T1= 3(a1+ a2)+ (2a1+ a2). Right: generalized armchair edge parallel to T1= (a1+ a2)+ 2(2a1+ a2).
is the sum of m zigzag vectors and one armchair vector. We call this a generalized zigzag edge (see Fig. 5). The supercell defined by the lattice vectors [Eq. (35)] contains
m+ 2 unit cells. For instance, values m = 1, 2 give angles θ =
79.1◦,70.9◦, respectively. The construction of the Hamiltonian matrix blocks defining the edge is straightforward. As an example, the Hamiltonian matrix blocks for m= 2, n = 1 are
H= ⎛ ⎜ ⎜ ⎜ ⎝ H0,0 H1,1 0 H−1,0e−i2πk H−1,−1 H0,0 H1,1 0 0 H−1,−1 H0,0 H1,1 H1,0ei2π k 0 H−1,−1 H0,0 ⎞ ⎟ ⎟ ⎟ ⎠, (37) B= ⎛ ⎜ ⎜ ⎜ ⎝ H0−1 H1,0 0 H−1,−1e−i2πk 0 H0,−1 H1,0 0 0 0 H0,−1 H1,0 0 0 0 H0,−1 ⎞ ⎟ ⎟ ⎟ ⎠. (38)
To generate edges with an orientation closer to the armchair edge (θ= 90◦), one can use the series with m= 1 and vary
n. The translation vector parallel to the edge is then the sum of one zigzag vector and n armchair vectors, which we call a generalized armchair edge (see Fig.5). The supercell then contains 2n+ 1 unit cells. As an example, m = 1, n = 3 gives
θ= 85.3◦. The edge termination can be controlled by adding and removing atoms at the edge, and apply the technique outlined in Sec.II A 5.
C. Tight-binding models
In most of the calculations shown here, we use tight-binding models with spin degeneracy, where spin-orbit coupling (SOC) is neglected. Including SOC changes the edge bands very little for MX2compounds containing light elements. For instance,
SOC gives a50 meV splitting in the edge bands of MoS2,
which is a small number relative to the dispersion of those bands. SOC gives a larger effect for heavy elements. The SOC-induced splitting in the edge band of WTe2 is 190 meV.
The effects of SOC are demonstrated in AppendixC. In the following, we proceed with Hamiltonians without SOC terms. The main contributions to the valence and conduction bands of MX2around the band gap come from the valence d shell of
the M atom and the p shell of the X atoms. A minimal basis set then comprises 11 orbitals: 5 metal d orbitals, and 6 p orbitals of the two chalcogen atoms [50]. Monolayer MX2 has D3h
point-group symmetry. The trigonal prismatic coordination of the metal atom splits the d states into three groups: A1{d3z2−r2},
E{dxy,dx2−y2}, and E
{dxz,dyz}. The six p orbitals of the two chalcogen atoms split into the groups A1{p−z =
pz(X1)− pz(X2)}, E {p+ x = px(X1)+ px(X2),p+y = py(X1)+ py(X2)}, A 2{pz+= pz(X1)+ pz(X2)}, and E{p−x = px(X1)− px(X2),p−y = py(X1)− py(X2)}. Mirror
symmetry in the plane of the metal atoms σh allows for hybridization between orbitals that are even with respect to
σh, i.e., A
1and E
orbitals, or between orbitals that are odd, i.e., A2and E orbitals.
The set of orbitals with even symmetry thus comprises the six orbitals d3z2−r2, dxy, dx2−y2, p−z, p+x, and py+, and the set
with odd symmetry the five orbitals dxz, dyz, pz+, p−x, and
p−y. As the even/odd symmetry is conserved for the edges and grain boundaries considered in this paper, all corresponding Hamiltonian matrices are blocked, and the even/odd solutions can be obtained separately. The matrix blocks Hp,qrequired for constructing the Hamiltonian matrices of Sec.II Bare given in AppendixB. The values of the tight-binding parameters have been obtained by fitting the bulk band structure to bands obtained from density functional theory (DFT) calculations with the generalized gradient approximation (GGA/PBE) functional [63]. For the even states we use the parameters given by Rostami et al. [64], and for the odd states we use the parameters given in AppendixB.
The 11-band model can be simplified further. From early theoretical studies and recent first-principles calculations, one observes that the MX2 bands at the top of the valence band
and at the bottom of the conduction band are dominated by the metal d orbitals, in particular those with even symmetry, i.e., d3z2−r2, dxy, and dx2−y2 [49,50,65]. Contributions to the
bands around the gap from the dxz and dyz orbitals and from the chalcogen p orbitals are much smaller. Matheiss has constructed an effective tight-binding model for MX2,
where only the metal sites are taken into account explicitly [49]. These sites form a two-dimensional triangular lattice, and the presence of the X atoms lowers the symmetry of this lattice from D6h to D3h. The metal orbitals with even
symmetry d3z2−r2, dxy, and dx2−y2 are used to construct an
effective three-band model. The matrix blocks Hp,q of this model are given in AppendixA. We use the parameters given by Liu et al. [34], which have been obtained by fitting the bulk bands to GGA/PBE results.
III. RESULTS A. Three-band model
1. Zigzag and armchair edges
The eigenvalues and eigenvectors of Eq. (8) are used to construct the Bloch matrices [Eq. (12)]. The bulk density of states n(k,E) is shown in Fig. 6(a), as calculated from the Green’s function matrix of a bulk strip [Eqs. (2) and (20)], in the zigzag edge orientation. In the three-band model, the lowest band is the valence band, and the overlapping two upper bands form the conduction band. Note that the zero of energy is chosen at the top of the valence band. The same Bloch matrices give access to the Green’s function matrices of the zigzag edge strips [Eqs. (18) and (19)]. Here, gR(k,E) and gL(k,E) are the energy and k-resolved Green’s function of the M edge and the X edge, respectively (see Fig.4). Again using MoS2 as an example, the associated densities of state nR(k,E) and nL(k,E) of the Mo edge, respectively the S edge, are shown in Figs.6(b)and6(c).
The Mo edge has a prominent edge state inside the bulk band gap dispersing upwards from k= 0 to12. The S edge has an edge state that starts in the bulk conduction band at k= 0, and disperses downward to k=12. In surface physics, such states are termed Shockley states [66]. These edge states are also found in the more detailed 11-band model with a similar
(a)
(b)
(c)
(d)
E (
e
V
)
-1 0 1 2 3 4 E (eV) 0 2 4 6 8 DoS (eV -1 ) -1 0 1 2 3 4 E (eV) 0 1 2 3 4 5 -1 0 1 2 3 4 E (eV) 0 1 2 3 4 -1 0 1 2 3 4 E (eV) 0 1 2 3 4 5 6(e)
(f)
(g)
(h)
k
k
k
k
FIG. 6. (a) k-resolved density of states (DOS) per MoS2unit of bulk MoS2in the three-band tight-binding model, with the k vector parallel to the zigzag edge, and the zero of energy at the top of the valence band. The DOS is plotted on a logarithmic scale (right-hand side n denotes amplitude 10−n) using a broadening parameter η= 0.05 eV [Eq. (1)]. (b), (c), (d) k-resolved DOSs of the Mo edge, the S edge, and the armchair edge, respectively. (e), (f), (g), (h) k-integrated DOS per MoS2unit of bulk MoS2, the Mo edge, the S edge, and the armchair edge, all plotted on a linear scale (η= 0.05 eV). The red solid lines give the counting function and the green dashed lines indicate the charge-neutrality level (CNL) [Eqs. (28)–(30)].
dispersion (see Sec. III B). The latter model shows a richer edge state structure, originating from bands that are omitted in the simple three-band model. Nevertheless, the three-band model finds prominent Mo-edge and S-edge states with the correct dispersion. Edge states are also found at energies close to∼3 eV, which is in the hybridization gap between the two conduction bands.
The structure of a bulk strip in armchair orientation has a mirror plane perpendicular to the MX2 plane and along the
armchair orientation (see Fig.4). It follows that the densities of states of right and left edges, nR(k,E) and nL(k,E), are identical for the armchair orientation. The k-resolved density of states of the MoS2 armchair edge [Eqs. (33) and (34)]
is shown in Fig. 6(d). There are two clear edge states with energies in the band gap. One edge state is just below the conduction band and roughly follows the dispersion of the conduction band edge, whereas the other one is positioned at ∼0.5 eV above the valence band, following the dispersion of valence band edge.
The k-integrated densities of states [Eq. (30)] of a bulk MoS2strip in zigzag orientation, the Mo edge, the S edge, and
the armchair edge, are given in Figs.6(e)–6(h), respectively. These densities of states show the van Hove singularities at the band edges that are typical of 1D structures. With the zero of energy at the top of the valence band, the bulk valence band lies in the energy range−0.5–0.0 eV, and the two conduction bands in the range 1.6–3.5 eV [see Fig.6(e)]. The Mo edge shows edge bands in the energy range 0.1–1.4 eV, and around 3 eV.
The S edge has an edge band starting at 0.5 eV, which merges with the conduction band at higher energies, and additional edge states with energies 2.7–3.0 eV. The armchair edge has two edge bands in the gap in the energy ranges 0.3–0.6 and 1.4–2.0 eV, respectively, and an edge state around 3.1 eV.
Also shown in Figs. 6(e)–6(h) are the electron counting function (red curves) [Eq. (29)] and the charge-neutrality level (CNL; green lines) [Eq. (28)]. Obviously, for the bulk the CNL is at the top of the valence band [see Fig.6(e)]. The position of the CNL at the Mo edge corresponds to a2
3 filled edge band in
the gap [see Fig.6(f)]. One might interpret this as effectively two-thirds of the bonds being broken for a Mo atom at the Mo edge. Likewise, the position of the CNL at the S edge corresponds to a 13 filled edge band in the gap [see Fig.6(g)]. This correlates with one-third of the bonds being broken for a Mo atom at the S edge. The supercell of the armchair edge has two Mo atoms along the edge, where the local environments of one is similar to that of a Mo atom at an Mo-edge, and the local environment of the other is similar to that of a Mo atom at the S edge (see Fig.4). Summing the23and13edge state occupations at the Mo and S edges, one predicts that the armchair edge has one completely filled edge band. The calculated CNL shows that this is indeed the case; the armchair edge has one fully occupied edge band, and one empty edge band [see Fig.6(h)]. The CNLs of all three basic edges are quite close [cf. Figs.6(f)–6(h)]. This is an artifact of the three-band model, and of disregarding the odd bands in particular. The latter play an important role in setting the CNLs in MX2edges, as
-1 0 1 2 3 4 E (eV) 0 1 2 3 4 5 DoS (eV -1 ) -1 0 1 2 3 4 E (eV) 0 1 2 3 4 5 (a) (b) (c) (d) k k E (eV)
FIG. 7. (a), (b) k-resolved DOSs of the W edge and the Te edge within the three-band tight-binding model of H−WTe2. (c), (d) k-integrated DOSs of the W edge and the Te edge. The same conventions and parameters are used as in Fig.6.
we will discuss in Sec.III B. In polar lattices such as MX2
the CNL does not fix the intrinsic Fermi level, unlike in a nonpolar lattice such as graphene [61,62]. In a finite-sized
MX2 sample, electrons can be redistributed among all the
edges in the sample, driven by the internal electric field set up by the intrinsic polarization of the material.
So far, we have focused on edges of MoS2, but all MX2
compounds with a similar structure (D3hpoint group, trigonal
prismatic coordination of M atoms by X atoms) have similar edge structures. Figure7shows the densities of states of the W and Te zigzag edges of H−WTe2 as an example. These
densities of states are very similar to those of the corresponding Mo and S edges of MoS2[see Figs.6(b),6(c),6(f), and6(g)].
The band gap of bulk WTe2 is somewhat smaller than that
of MoS2, and the band widths are somewhat larger. WTe2
has a band gap of 1.1 eV (GGA/PBE), a valence band in the range−0.7–0.0 eV, and conduction bands in the range 1.1–3.7 eV. The W edge has a state in the gap in the energy range 0.0–1.4 eV, and the Te edge has a state in the gap in the range 0.5–2.0 eV. As is the case for the bulk bands, the band widths of these edge states are somewhat larger than the corresponding states in MoS2.
2. General edges and grain boundaries
The electronic structure of the more general edges defined in Sec.II B 2can be calculated along the same lines. Figures8(a)
and 8(b) give the densities of states of generalized zigzag edges of MoS2with translation vectors along the edge defined
by m= 1 and 3, respectively, and n = 1 [see Eq. (35) and Fig.5]. The generalized zigzag edge has a rich structure of edge states within the bulk band gap. The peak at 0.4 eV in
0 2 4 E (eV) 0 1 2 3 4 5 0 2 4 E (eV) 0 1 2 3 4 5 -1 0 1 2 3 4 E (eV) 0 1 2 3 4 5 DoS (eV -1) -1 0 1 2 3 4 E (eV) 0 1 2 3 4 5 DoS (eV -1) (a) (b) (c) (d)
FIG. 8. k-integrated DOSs of the generalized zigzag edges defined by Eq. (35) with (a) m= 1, n = 1 and (b) m = 3, n = 1, and of the generalized armchair edges with (c) m= 1, n = 2 and (d)
m= 1, n = 3.
Figs.8(a)and8(b)results from the armchair part of this edge [compare to Fig. 6(h)]. The band in the range 0.5–1 eV in
8(a), and the structure of bands in the range 0.8–1.9 eV in8(b)
originates mainly from the zigzag parts.
The counting system for the three-band model, as outlined in the previous section, gives a filling23 per Mo atom at a Mo edge and a filling 1 per two Mo atoms at an armchair edge. For a generalized zigzag edge one would then predict the CNL to correspond to 2m/3+ n filled edge states. This would mean that the CNL is inside an edge band, unless m is a multiple of three. This is confirmed by a calculation of the CNL according to Eq. (28). In Fig.8(a), where m= 1, the position of the CNL corresponds to a23filled edge band, whereas in Fig.8(b), where
m= 3, the CNL is in a gap between two edge states. The gaps
between the edge states of the general edges are quite small, however. In a sample that involves long-range charge transfer between the edges, as discussed above, these semiconducting edges easily become doped.
Figures 8(c) and 8(d) show the densities of states of generalized armchair edges of MoS2, with translation vectors
along the edge with m= 1 and n = 2 and 3, respectively [see Eq. (35) and Fig. 5]. Again, this general edge has a rich structure of states within the bulk band gap. In particular, the peaks around 0.4 and 1.7 eV have a strong armchair character. The counting model gives 1/3+ n filled edge states, so it predicts that the CNL always lies inside an edge band. Calculations of the CNL according to Eq. (28) confirm this, as shown in Figs. 8(c) and 8(d). In conclusion, within the three-band tight-binding model, pristine charge-neutral edges are metallic for most edge orientations.
A similar analysis can be applied to the states found at grain boundaries. We illustrate the use of Eq. (21) for calculating
-1 0 1 2 3 4 E (eV) 0 1 2 3 DoS (eV -1 ) -1 0 1 2 3 4 E (eV) 0 1 2 3 k k E (eV) (a) (b) (c) (d)
FIG. 9. (a) k-resolved DOS of a model grain boundary between a Mo edge and a S edge, with a coupling between the edges
α= 0.2 [Eq. (32)]. (b) k-resolved DOS with α= 0.8. (c), (d) The corresponding k-integrated DOSs.
the states at grain boundaries using a simple grain boundary in MoS2. It consists of a left semi-infinite part, terminated by
a S edge, connected to a right semi-infinite part, terminated by a Mo edge. As hopping matrix defining the coupling between left and right parts (Fig.2), we choose a scaled version of the bulk hopping matrix [Eq. (32)]:
B= αB; 0 α 1. (39)
Obviously, α= 0 gives two uncoupled S and Mo edges with corresponding edge states [Figs. 6(b),6(c), 6(f), and 6(g)], whereas α= 1 gives the bulk electronic structure [Figs.6(a)
and 6(e)] without any edge states. The values 0 < α < 1 then represent a simple model for a weak link with zigzag orientation.
Figure9(a)gives the k-resolved density of states for α= 0.2, which represents a relatively weak coupling between the left and right parts. In the band gap one observes the two edge bands that are typical of the Mo edge and the S edge [see Figs.6(b)and6(c)]. Due to the coupling between the S and Mo edges at the grain boundary, the two bands interact, which results in an avoided crossing between the two edge states at k≈ 0.35 and E ≈ 1 eV. The avoided crossing creates an energy gap in the range 0.7–1.2 eV, which is clearly visible in the k-integrated density of states for α= 0.2, shown in Fig.9(c). The counting model predicts23, respectively13, filling per Mo atom of an edge state at the Mo edge and the S edge, if the grain boundary is charge neutral. This implies that the lower edge band is fully occupied, whereas the upper edge band is empty. A calculation of the CNL according to Eq. (28) confirms this [see Fig.9(c)].
Figure 9(b) gives the k-resolved density of states for a strong coupling between the left and right parts, α= 0.8. The
interaction between the two edge bands is now much stronger than for the case discussed in the previous paragraph, which results in a much larger gap in the range 0.2–1.6 eV [see Fig.9(d)]. The two edge bands are pushed toward the valence band and the conduction band, respectively, and they more closely follow the dispersions of the valence and conduction band edges. The occupancy of a charge-neutral grain boundary is the same as before, i.e., a fully occupied lower edge band and an empty upper edge band.
As a final point in this section we illustrate the technique introduced in Sec.II A 5 for handling edges that involve an electronic or a structural reconstruction. For the zigzag edges, we have found a 23, respectively 13, occupied edge band for the Mo edge and the S edge. It suggests that a reconstruction that triples the translational period along the edge may lead to fully occupied edge states for both edges. As a proof of principle, we test a very simple reconstruction, where one in three Mo atoms at the Mo edge or the S edge has a different onsite energy. The Green’s function matrix is calculated from Eq. (25). The densities of states of the modified edge layers are shown in Figs.10(a)and10(c)for the Mo edge and10(b)
and10(d)for the S edge. Obviously increasing the periodicity from 1× to 3× folds the edge band, so one observes three edge states instead of one. Decreasing the onsite energy of every third Mo atom at the edge by 1 eV introduces energy gaps between the edge states.
By calculating the CNL, we observe that at the Mo edge the two lowest edge bands, starting at 0.4 and 0.7 eV in Figs.10(b)
and10(d), are filled, whereas the third band, starting at 1.1 eV, is empty. The lowest edge band of the S edge at 0.0 eV [Figs.10(a)and10(c)] is filled, whereas the two upper bands, starting at 0.9 and 1.6 eV, respectively, are empty.
FIG. 10. (a), (c) k-resolved and k-integrated DOS of a Mo edge with a 3× reconstruction, where the onsite energy of every third Mo atom is 1 eV lower than that of the other two; (b), (d) the same information for the S edge.
B. 11-band tight-binding model
As discussed in Sec.II C, the 11-band tight-binding model of MX2uses a basis set composed of the 5 valence d orbitals
of the M atom and the 6 valence p orbitals of the two
X atoms. Because mirror symmetry σh in the MX2 plane
holds both for monolayers as well as for edges, states that are even or odd with respect to this mirror symmetry are treated separately. The even symmetry set comprises the six orbitals d3z2−r2, dxy, dx2−y2, pz−, p+x, and p+y, and the odd
symmetry set the five orbitals dxz, dyz, p+z, p−x, and p−y, leading to six- and five-dimensional Hamiltonian matrices, respectively.
Again, we use MoS2as an example of a MX2monolayer.
The procedures for obtaining the Green’s functions and the densities of states are the same as those described in the prev-ious section. In the following, we will only discuss the results obtained for the basic edges, the zigzag and the armchair edges, and compare those to the results obtained with the three-band model.
1. Even bands
The k-resolved density of states of the even states of a bulk strip in the zigzag edge orientation is shown in Fig.11(a)in an energy region around the band gap, and the corresponding
k-integrated density of states is given in Fig.11(e). As before, the zero of energy is positioned at the top of the valence band. Qualitatively, these densities of states are similar to those of the three-band model [compare to Figs.6(a)and6(e)].
Quantitatively, the 11-band model gives a highest valence band that is wider by∼0.3 eV, whereas the two lowest conduction bands are narrower by∼0.5 eV in total.
Figures 11(b) and 11(f) show the resolved and k-integrated densities of states of the Mo edge. There are two prominent edge bands with energies in the MoS2band gap. One
edge band emerges from the bulk valence band at k≈ 0.3, and disperses upward with increasing k to 1 eV. A second band disperses downward from the conduction band to 1.6 eV. The first edge band is also found in the three-band model [compare Figs.6(b)and6(f)]. In the three-band model. it is found at a slightly higher energy, such that it is completely isolated from the bulk valence band. The second edge band is absent in the three-band model. We will discuss the character of these bands in more detail below. Like in the three-band model, there are also edge states at other energies, for instance, just below the highest valence band at∼−1.1 eV, and in the hybridization gaps of the conduction bands.
The k-resolved and k-integrated densities of states of the S edge are given in Figs.11(c)and11(g). At energies in the MoS2band gap there is one prominent edge state, which starts
at k= 0 in the conduction band and disperses downward with increasing k to 0.7 eV. The three-band model shows the same edge state with more or less the same dispersion [compare to Figs. 6(c) and 6(g)]. In the 11-band model this state is somewhat more prominently isolated from the conduction band. Unlike for the Mo edge, the 11-band model does not give additional even edge states for the S edge in the band gap, as compared to the 3-band model.
-1 0 1 2 3 E (eV) 0 1 2 3 4 DoS (eV -1 ) -1 0 1 2 3 E (eV) 0 1 2 3 4 -1 0 1 2 3 E (eV) 0 1 2 3 -1 0 1 2 3 E (eV) 0 1 2 3 (a) (b) (e) (f) (g) (h) E (e V ) k k k k (c) (d)
FIG. 11. As Fig.6but for the states within the 11-band model with even symmetry with respect to a mirror plane through the plane of the Mo atoms. (a), (b), (c), (d) k-resolved DOS per MoS2 unit of bulk MoS2, the Mo edge, the S edge, and the armchair edge, respectively. (e), (f), (g), (h) The corresponding k-integrated DOSs of bulk MoS2, the Mo edge, the S edge, and the armchair edge. The same conventions and parameters are used as in Fig.6.
Finally, Figs.11(d)and11(h)show the resolved and k-integrated densities of states of the armchair edge. There are three edge bands with energies in the MoS2 band gap. One
band is situated at−0.1–0.5 eV and a second edge band lies at 0.9–1.5 eV. These two edge bands roughly correspond to the ones that are found in the three-band model of the armchair edge [see Figs.6(d)and6(h)]. In the 11-band model these two edge bands are found at a slightly lower energy in the gap. In addition, the present model finds a third edge band in the gap, just below the conduction band at 1.7–1.9 eV.
In conclusion, although there are quantitative differences between the electronic structures found with the 11-band model and the 3-band model, qualitatively they give similar results concerning the prominent edge states of even symmetry found in the MoS2band gap. For the Mo edge and the armchair
edge, the 11-band model gives an additional edge state, as compared to the 3-band model, with energies just below the conduction band.
Figures 11(e)–11(h) also show the CNLs. One needs of course all the bands to calculate the CNLs, including the odd bands to be discussed in the next section. A comparison to Fig.6 reveals that the odd bands are in fact instrumental in fixing the CNLs. For instance, the CNLs of the Mo edge and
the S edge differ by 1.4 eV [Figs.11(f)and11(g)], whereas the corresponding CNLs in Figs.6(f)and6(g)differ by 0.3 eV only. The 11-band model gives one completely occupied even edge band for the Mo edge one, one empty one [see Fig.11(f)]. For the S edge it gives one completely empty even edge band [see Fig.11(g)]. This is unlike the the three-band model, where we found partially filled even edge bands, both for the Mo edge and the S edge. As we will see in the next section, edge states of odd symmetry, which are absent in the three-band model, pin the CNLs of the Mo edge and the S edge, and make these edges metallic. The CNL at the armchair edge is in the gap between the two edge bands [see Fig.11(h)], as it also is in the three-band model [see Fig.6(h)].
The orbital character of the edge states can be analyzed using the projected density of states, according to Eq. (5). The projections on the Mo d3z2−r2 orbital, the dxyand dx2−y2
orbitals, and the p orbitals of the S atoms are given in Figs.12(a)–12(c), respectively, for the Mo edge. The results indicate that both edge states in the MoS2 band gap have a
mixed Mo and S character. The dominant Mo contributions clearly come from the dxyand dx2−y2orbitals. The upper edge
state has S p character mixed in, in particular for k < 0.4, whereas the lower edge state has some S p character mixed
FIG. 12. Projected densities of states (PDOS) calculated according to Eq. (5). (a)–(c) For the Mo edge, projected on the Mo d3z2−r2orbital,
the dxy and dx2−y2 orbitals, and the p orbitals of the S atoms, respectively. (d)–(f) The same information for the S edge. The PDOSs are
FIG. 13. As Fig.11, but for the states with odd symmetry with respect to a mirror plane through the plane of the Mo atoms.
in for k > 0.4. There is little d3z2−r2character mixed in these
edge bands, except at the band edges.
The projected densities of states of the S edge, projected on the same orbitals, are given in Figs.12(d)–12(f). Also here, the edge state in the MoS2 band gap has a mixed Mo and S
character. As for the Mo contribution, for k= 0 the dominant character is Mo dxy and dx2−y2. That changes somewhat for
larger k; at k= 0.5 the dominant character is Mo d3z2−r2. The
contribution of the S orbitals is fairly constant throughout the whole edge band.
2. Odd bands
The k-resolved and k-integrated bulk densities of states
n(k,E) corresponding to the odd states are given in Figs.13(a)
and 13(e). The band gap between the odd states is 3.3 eV, which is significantly larger than the gap between the even states. The top of the highest valence band of the odd states is approximately 0.8 eV below the top of the highest valence band of the even states, whereas the bottom of the lowest conduction band of the odd states is approximately 0.6 eV above that of the even states.
The k-resolved densities of state of the odd bands are shown in Figs.13(b)and13(c)for the Mo edge and of the S edge, respectively, and the corresponding k-integrated densities of states are shown in Figs.13(f)and13(g). Both edges have a prominent edge band with moderate dispersion inside the gap between the odd states. The edge band of the Mo edge lies close to the conduction band between 1.3 and 1.8 eV, whereas the edge band of the S edge is close to the top of the valence band between−0.5 and 0.2 eV.
If we also take the bulk band structure of the even bands into consideration [Figs. 11(a)and11(d)], then most of the
edge band at the S edge [Figs.13(c)and13(g)] overlaps with the valence band. It is still a true edge state, though, because the interaction between the odd edge state and the even bulk bands is symmetry forbidden. The edge band at the Mo edge [Figs.13(b)and13(f)] partially overlaps with the conduction band of the even states. Also, this is a true edge state over the whole 1D Brillouin zone, as interaction with the bulk states is symmetry forbidden.
Figures13(d)and13(h)give the k-resolved and k-integrated densities of states corresponding to the armchair edge. These results show two edge bands in the gap region with a modest to small dispersion. The lowest of these edge bands has a dispersion of 0.4 eV, and taking also the even bulk bands into consideration, it lies in the MoS2 valence band. The
upper edge band is near dispersionless; it is found in the gap at 1.6 eV.
The calculated CNLs are also given in Figs.13(e)–13(h). From these figures it becomes clear that the odd edge states are instrumental in controlling the CNLs of the edges. The CNL of the Mo edge is found in the edge band that is close to the MoS2 conduction band [see Fig.13(f)], and the CNL
of the S edge lies in the edge band close to the MoS2valence
band [see Fig.13(g)]. Within the 11-band model both these edges are metallic in their charge-neutral state. The CNL of the armchair edge lies between the two edge bands [Fig.13(h)], which means that the neutral armchair edge is semiconducting. The orbital character of the odd symmetry edge states is identified using the projected densities of state. Figures14(a)
and14(b) show the projections on the d and p orbitals of the density of states of the Mo edge, and Figs. 14(c) and
14(d)show the projections for the S edge. The edge state at the Mo edge has somewhat of a mixed character, but it is
FIG. 14. As Fig.12, but for the states with odd symmetry with respect to a mirror plane through the plane of the Mo atoms.
dominated by Mo d orbitals. The edge state at the S edge has a clear S p character. More than the even edge states discussed in the previous section, these two odd edge states have a “dangling-bond” character, i.e., they correspond to bulk bonds that are broken at the edge atoms.
IV. SUMMARY AND CONCLUSIONS
Growth of 2D semiconductors, such as the transition metal dichalcogenides MX2, commonly results in sheets that contain
quasi-1D metallic structures, which are located at the edges of 2D crystallites or at the boundaries between 2D crystal grains. Such edges and grain boundaries form a practical realization of 1D metals. Metallic edges have also been identified as the active sites in MX2 catalysts. A standard technique for
modeling edge or grain boundary states starts from nanoribbon structures, which involves using large supercells if one aims at obtaining converged results that enable the separation between edge and bulk properties. Moreover, as a nanoribbon has two edges, electron transfer between the edges can occur in order to equilibrate the system, which complicates identifying the properties of single edges.
In this paper, we formulate a Green’s function technique for calculating edge states of (semi-)infinite 2D systems with a single well-defined edge or grain boundary. We express bulk, edge, and boundary Green’s functions in terms of Bloch matrices. The latter are constructed from the solutions of a quadratic eigenvalue equation, which gives the traveling and evanescent Bloch states of the system. Electronic and structural reconstructions of edges or grain boundaries are easily incorporated in the technique. The formalism can be implemented in any localized basis set representation of the Hamiltonian.
Here, we use it to calculate edge and grain boundary states of MX2 monolayers by means of tight-binding models. A
simple three-band model is employed to explore the electronic structures of the basic pristine zigzag and armchair edges in MX2. Within this model, the zigzag edges are metallic
in their charge-neutral state, whereas the armchair edge is semiconducting. The three-band model is also used to analyze the electronic structures of MoS2 edges with a more general
orientation between zigzag and armchair; these edges are generally metallic. The same model is applied to obtain electronic structures of grain boundaries and of structurally modifed edges. The inclusion of spin-orbit coupling leads to a splitting of the edge states, which ranges from50 meV in MoS2to190 meV in WTe2. This moderate splitting does not
alter the basic character of the edge states.
The three-band model captures only part of MX2 edge
electronic structure. A more complete picture of the rich electronic structure of MX2edges is obtained from an 11-band
tight-binding model comprising the d metal valence orbitals and the p valence orbitals of the chalcogen atoms. Edge states in the band gap with even symmetry (with respect to a mirror through the plane of M atoms) are qualitatively similar to those found in the three-band model. As in the latter model, charge-neutral zigzag edges are metallic and the armchair edge is semiconducting in the 11-band model. Unlike the three-band model, however, the charge-neutrality level is fixed by edge states of odd symmetry in the band gap. The odd states have a character that reflects the dangling bonds on the edge atoms, unlike the even states, which have more of a mixed M−d and
X−p character. The odd edge states are likely to be found near
the Fermi level under experimental conditions, and can not be discarded when modeling MX2edges.
ACKNOWLEDGMENT
This work is part of the research program of the Foundation for Fundamental Research on Matter (FOM), which is part of the Netherlands Organization for Scientific Research (NWO).
APPENDIX A: THREE-BAND TIGHT-BINDING MODEL
Following Mattheis, the three-band tight-binding model discussed in Sec.II Conly explicitly uses the sublattice of the
Matoms and the d orbitals of A1and Esymmetry, i.e., the orbitals that are symmetric with respect to σh, mirror symmetry in the MoS2plane [49]. We number these orbitals as
d3z2−r2 = |0 , dxy= |1 , dx2−y2= |2 . (A1)
The matrix blocks Hp,q discussed in Sec.II Bare then 3× 3 matrices, where Hp,q denotes the real-space Hamiltonian matrix block pertaining to the interaction between atoms in the unit cell at the origin and atoms in the unit cell at pa1+ qa2.
Using only nearest-neighbor interactions, one has|p|,|q| 1. Defining the tight-binding parameters i = i|h|i , i = 0,1,2, and tij = i|h|j , i,j = 0,1,2, and making use of the rotation properties of the d orbitals, one derives the matrix
blocks H0,0= ⎡ ⎢ ⎣ 0 0 0 0 1 0 0 0 2 ⎤ ⎥ ⎦; H1,0= ⎡ ⎢ ⎢ ⎣ t00 √3 2 t02+ 1 2t01 − 1 2t02+ √ 3 2 t01 √ 3 2 t02− 1 2t01 1 4t11+ 3 4t22 √ 3 4 (t11− t22)− t12 −1 2t02− √ 3 2 t01 √ 3 4 (t11− t22)+ t12 1 4t22+ 3 4t11 ⎤ ⎥ ⎥ ⎦; (A2) H0,1= ⎡ ⎢ ⎣ t00 −t01 t02 t01 t11 −t12 t02 t12 t22 ⎤ ⎥ ⎦; H1,1= ⎡ ⎢ ⎢ ⎣ t00 − √ 3 2 t02− 1 2t01 − 1 2t02+ √ 3 2 t01 −√3 2 t02+ 1 2t01 1 4t11+ 3 4t22 − √ 3 4 (t11− t22)+ t12 −1 2t02− √ 3 2 t01 − √ 3 4 (t11− t22)− t12 1 4t22+ 3 4t11 ⎤ ⎥ ⎥ ⎦. (A3)
The hopping associated with matrix elements t00,t02,t22,t11
is even with respect to inversion and the one associated with
t01,t12is odd. It is then easy to see that
H−p,−q= (Hp,q)T. (A4) The values of the tight-binding parameters are taken from Liu et al. [34], which have been obtained by fitting the tight-binding band structure of a MX2sheet to a DFT GGA/PBE
band structure [63], using the highest valence band and the two lowest conduction bands in the fit. GGA/PBE gives an optimized lattice parameter of 3.19 ˚A for MoS2, which is 1%–
2% larger than the reported experimental values [67–69]. With this lattice parameter, the calculated band gap is 1.63 eV, which is smaller than the experimental optical band gap of 1.85 eV [70].
APPENDIX B: 11-BAND TIGHT-BINDING MODEL
The 11-band tight-binding model of MX2uses a basis set
composed of all 5 d orbitals of the M atom, and the 6 p orbitals of the two X atoms. Following the approach of Cappelluti
et al. [50,71], we include next-nearest-neighbor interactions, and use the Slater-Koster two-center approximation for the hopping matrix elements [72]. The real-space matrix blocks
Hp,qdiscussed in Sec.II Bthen have the form
H0,0 = d t0,0dp t0,0dpT p , H1,0= t1,0dd 0 0 t1,0pp , H0,1= t0,1dd 0 t0,1dp T t0,1pp , H1,1= t1,1dd 0 (t1,1dp) T t1,1pp , (B1) with the remaining blocks constructed according to Eq. (A4). The bulk Hamiltonian for an infinite layer with two-dimensional translation symmetry is written in this notation as
H(k1,k2)= H0,0+ A + A†;
A= H1,0ei2π k1+ H0,1ei2π k2+ H1,1ei2π (k1+k2). (B2)
The parameters in the model can be found by fitting the tight-binding band structure obtained with this bulk Hamiltonian to a band structure obtained from a DFT calculation. The model turns out to be too restrictive to obtain a satisfactory fit for all 11 bands with one parameter set. A good fit can, however, be obtained if we divide the bands into a set of even symmetry and a set of odd symmetry, and use different parameters for the two sets. Mirror symmetry σhin the MX2plane holds for
mono-layers, as well as for edges, so states that are even or odd with respect to σhcan be treated separately. In AppendixB 1we give the expressions for the matrix blocksa and tp,qab in Eq. (B1) for the even states and in AppendixB 2for the odd states.
1. Even states
The even states are composed of orbitals of E and A1 symmetry: d3z2−r2 = |0 , dxy= |1 , dx2−y2= |2 , 1 √ 2[px(X1)+ px(X2)]= |3 , 1 √ 2[py(X1)+ py(X2)]= |4 , 1 √ 2[pz(X1)− pz(X2)]= |5 . (B3) The matrix blocksaand tp,qab in Eq. (B1) are then all 3× 3:
d = ⎡ ⎢ ⎣ 3z2−r2 0 0 0 xy 0 0 0 xy ⎤ ⎥ ⎦, p= ⎡ ⎢ ⎣ px+ Vppπ 0 0 0 px+ Vppπ 0 0 0 pz − Vppσ ⎤ ⎥ ⎦, (B4)
where aare the onsite orbital energies and Vabαare the Slater-Koster two-center integrals [72]. Similarly,
t0,0dp = √ 2 7√7 ⎡ ⎢ ⎣ 9Vpdπ− √ 3Vpdσ 3 √ 3Vpdπ− Vpdσ 12Vpdπ+ √ 3Vpdσ −Vpdπ− 3 √ 3Vpdσ −5 √ 3Vpdπ− 3Vpdσ −6Vpdπ+ 3 √ 3Vpdσ −5√3Vpdπ− 3Vpdσ 9Vpdπ − √ 3Vpdσ −2 √ 3Vpdπ+ 3Vpdσ ⎤ ⎥ ⎦, (B5)
t1,0dd = ⎡ ⎢ ⎢ ⎣ 1 4(3Vddδ+ Vddσ) 3 8(−Vddδ+ Vddσ) √ 3 8 (−Vddδ+ Vddσ) 3 8(−Vddδ+ Vddσ) 1 16(3Vddδ+ 4Vppπ+ 9Vddσ) √ 3 16(Vddδ− 4Vddπ + 3Vddσ) √ 3 8 (−Vddδ+ Vddσ) √ 3 16(Vddδ− 4Vddπ+ 3Vddσ) 1 16(Vddδ+ 12Vppπ+ 3Vddσ) ⎤ ⎥ ⎥ ⎦, (B6) t1,0pp = ⎡ ⎢ ⎣ 1 4(3Vppπ+ Vppσ) √ 3 4 (Vppπ− Vppσ) 0 √ 3 4 (Vppπ− Vppσ) 1 4(Vppπ+ 3Vppσ) 0 0 0 Vppπ ⎤ ⎥ ⎦, (B7) t0,1dd = ⎡ ⎢ ⎢ ⎣ 1 4(3Vddδ+ Vddσ) 3 8(Vddδ− Vddσ) √ 3 8 (−Vddδ+ Vddσ) 3 8(Vddδ− Vddσ) 1 16(3Vddδ+ 4Vppπ+ 9Vddσ) √ 3 16(−Vddδ+ 4Vddπ− 3Vddσ) √ 3 8 (−Vddδ+ Vddσ) √ 3 16(−Vddδ+ 4Vddπ− 3Vddσ) 1 16(Vddδ+ 12Vppπ+ 3Vddσ) ⎤ ⎥ ⎥ ⎦, (B8) t0,1pp = ⎡ ⎢ ⎣ 1 4(3Vppπ+ Vppσ) √ 3 4 (−Vppπ+ Vppσ) 0 √ 3 4 (−Vppπ+ Vppσ) 1 4(Vppπ+ 3Vppσ) 0 0 0 Vppπ ⎤ ⎥ ⎦, (B9) t0,1dp = 2√2 7√7 ⎡ ⎢ ⎣ 0 −3√3Vpdπ+ Vpdσ 6Vpdπ+12 √ 3Vpdσ Vpdπ 0 0 0 −3Vpdπ− 2 √ 3Vpdσ 2 √ 3Vpdπ− 3Vpdσ ⎤ ⎥ ⎦, t1,1dd = ⎡ ⎢ ⎣ 1 4(3Vddδ+ Vddσ) 0 √ 3 4 (Vddδ− Vddσ) 0 Vddπ 0 √ 3 4 (Vddδ− Vddσ) 0 1 4(Vddδ+ 3Vddσ) ⎤ ⎥ ⎦, t1,1 pp = ⎡ ⎢ ⎣ Vppσ 0 0 0 Vppπ 0 0 0 Vppπ ⎤ ⎥ ⎦, (B10) t1,1dp = √ 2 7√7 ⎡ ⎢ ⎣ −9Vpdπ+ √ 3Vpdσ 3 √ 3Vpdπ− Vpdσ 12Vpdπ+ √ 3Vpdσ −Vpdπ− 3 √ 3Vpdσ 5 √ 3Vpdπ + 3Vpdσ 6Vpdπ− 3 √ 3Vpdσ 5√3Vpdπ+ 3Vpdσ 9Vpdπ− √ 3Vpdσ −2 √ 3Vpdπ+ 3Vpdσ ⎤ ⎥ ⎦. (B11)
We use values of the parameters as obtained by Rostami et al. from fitting the even tight-binding bands to DFT GGA/PBE bands [64]. A lattice parameter of 3.16 ˚A for MoS2 has
been used in these calculations, which is 1% smaller than the optimized GGA/PBE value, but close to the reported experimental values. It results in a calculated band gap of 1.76 eV, which is slightly smaller than the experimental optical band gap of 1.85 eV [70]. We have shifted the tight-binding
bands such that the zero of energy coincides with the top of the valence band.
2. Odd states
The odd states are composed of orbitals of E and A1 symmetry: dxz= |1 , dyz= |2 , 1 √ 2[px(X1)− px(X2)]= |3 , √1 2 py(X1)− py(X2) = |4 , √1 2[pz(X1)+ pz(X2)]= |5 . (B12) The matrix blocksain Eq. (B1) are 2× 2 if a = d and 3 × 3 if a = p, whereas the matrix blocks tp,qab are 2× 2 if ab = dd,3 × 3 if ab= pp, and 2 × 3 if ab = dp: d = xz 0 0 xz , p= ⎡ ⎢ ⎣ px− Vppπ 0 0 0 p x− V ppπ 0 0 0 p z + V ppσ ⎤ ⎥ ⎦, (B13) t0,0dp = √ 2 7√7 √ 3Vpdπ + 9Vpdσ −6Vpdπ + 3√3Vpdσ −√3Vpdπ − 9Vpdσ −6V pdπ+ 3 √ 3Vpdσ 5 √ 3Vpdπ + 3Vpdσ −Vpdπ − 3 √ 3Vpdσ , (B14)
