Flat bands in the Weaire–Thorpe model and silicene

After the physics of graphene was initiated, originally by a theoretical prediction for a massless Dirac fermion by Wallace [1] back in the 1950s, and then by a recent experimental realization [2], interests have extended to wider class of system. Unusual behaviours of the massless Dirac fermions have since been experimentally confirmed and theoretically analysed [3]. History seems to repeat itself, for silicene: the material was theoretically predicted in the early 1990s by one of the present authors and Takeda [4], and the material began to be synthesised recently [5–7], after a long latent period. Successful synthesis of silicene on Ag has triggered serious attention in both physics and materials science. In graphene, the carriers, being in an atomically flat system, are π-electrons arising from the sp² hybridization of carbon orbitals. In silicene, by contrast, the honeycomb lattice is buckled, so that an sp³ character of Si is involved, which is crucial as pointed out previously in [4]. Thus silicene is not just a Si analogue of graphene, but distinct multi-orbital characters should appear in its electronic structure.

Hence silicene has a larger degree of freedom than graphene as a target material for applications and also theoretical considerations. In the present paper, we focus on this multi-orbital feature of silicene, where we opt for a simplified model, namely we propose an extension of the Weaire–Thorpe (WT) model, which was originally conceived for three-dimensional (3D) silicon with sp³ orbits in a tight-binding model on a diamond lattice [8]. In the original WT model for four sp³ orbitals on each atom, along with inter-atom hopping in 3D, the band structure is systematically examined, and two singular dispersions arise: one is a massless Dirac cone, and the other is a dispersionless (flat) band [9]. Massless Dirac cones in three dimensions are topologically protected. In two dimensions they can arise with supplemental symmetry protections, as in graphene and silicene [10]. On the other hand, the flat bands are not protected by symmetry in realistic materials. Still, however, the flat bands can have some topological/geometrical origins reflecting the multi-orbital character of a given material as we describe in the present paper. For silicene, here we re-formulate the WT model in terms of the overlapping molecular-orbital (MO) method [11] discussed by Hatsugai and Maruyama, which contains WT and enables us to describe flat bands away from the electron–hole symmetric point. The overlapping MO formulation indeed enables us to pinpoint, algebraically, an important difference: while in 3D the dipersive bands with cones are sandwiched by doubly-degenerate flat bands, in 2D the dipersive bands with cones are sandwiched by triply-degenerate and non-degenerate (nearly) flat bands, which is consistent with the original band calculation by Takeda and Shiraishi [4]. Thus there emerges an explanation of why the whole band structure of silicene comprises dispersive bands with Dirac cones along with nearly flat (narrow) bands.

In the present paper, we first start with an overlapping molecular orbital theory [11]. Applying this to the WT model enables us to generically treat the flat bands away from the electron–hole symmetric point in multi-orbital models while the usual flat-band theories [12–14] focus on those at the electron–hole symmetric point. We then give a picture for silicene in this formalism, and compare the result with the band structure of Takeda and Shiraishi. We finally discuss that the nearly flat bands imply a large density of states (DOS), which would offer an interesting possibility for inducing instabilities into symmetry-broken states, such as ferromagnetism or structural changes. For the band engineering, the bonds perpendicular to the atomic plane are suggested to be crucial.

Let us start with describing a class of lattice model Hamiltonians in terms of molecular orbitals proposed by Hatsugai and Maruyama [11], which is explained here to make the present paper self-contained. Consider fermions on an N-site lattice with the creation operator, $c_{i}^{\dagger }$, for the sites $i=1,\cdots ,N$ with $\{c_{i}^{\dagger },{{c}_{j}}\}={{\delta }_{ij}}$. After shifting the origin of energy by μ, let us consider the case in which the Hamiltonian is given as a sum of overlapping molecular orbitals $m(=1,...,M)$ as

$$\begin{eqnarray*}&&H-\mu \mathcal{N}=\mathop{\sum }\limits_{m=1}^{M}{{\mathcal{E}}_{m}}C_{m}^{\dagger }{{C}_{m}},\end{eqnarray*}$$

where $\mathcal{N}=\sum _{i=1}^{N}c_{i}^{\dagger }{{c}_{i}}$ is the number operator of the fermions. The coefficient ${{\mathcal{E}}_{m}}\in \mathbb{R}$ is the energy level of the molecular orbital m, while $C_{m}^{\dagger }$ creates the molecular orbital as

$$\begin{eqnarray*}&&C_{m}^{\dagger }=\mathop{\sum }\limits_{i=1}^{N}c_{i}^{\dagger }{{\psi }_{i,m}}={{c}^{\dagger }}{{\psi }_{m}},{{c}^{\dagger }}=(c_{1}^{\dagger },\cdots ,c_{N}^{\dagger }),{{\psi }_{m}}=\left[ \begin{array}{ccccccccccccccc} {{\psi }_{1,m}} \\ \vdots \\ {{\psi }_{N,m}} \\ \end{array} \right],\end{eqnarray*}$$

where ${{\psi }_{m}}$ is the wavefunction of the molecular orbital m. We do not require translational symmetry in the system or in the molecular orbitals as is the case with the WT model.

Then we have a simple theorem that $H-\mu \mathcal{N}$ has $(N-M)$-fold degenerate zero-energy states when $N-M\gt 0$, where M is the total number of molecular orbitals in the whole real space. The number of zero energy states $N-M$ can be macroscopic. Since the theorem is general, we can apply it to periodic systems as well. Then we can characterize the wavefunctions as Bloch states, and N may be regarded as the number of energy bands when we express the Hamiltonian in the Bloch basis. The total number of molecular orbitals, M, in this Bloch basis corresponds to the total number of terms in the Hamiltonian in the momentum representation, which varies from one model to another, as we shall see.

Note here that if there are zero-energy states in the present Hamiltonian, they should be located at $E=\mu$ for the original Hamiltonian H. Thus we can describe non-zero energy states algebraically by suitably choosing μ. This is trivial but useful as we demonstrate in sections 3 and 4.

While we can normalize the molecular orbitals as $\psi _{m}^{\dagger }\;{{\psi }_{m}}=1$, they overlap with each other in general, which implies that the anticommutation relation $\{{{C}_{m}},C_{{{m}^{\prime }}}^{\dagger }\}$ may not be simple. We further decompose the Hamiltonian as

$$\begin{eqnarray*}&&H-\mu \mathcal{N}={{{\boldsymbol{c}} }^{\dagger }}h{\boldsymbol{c}} ,\qquad h=\mathop{\sum }\limits_{m=1}^{M}{{\mathcal{E}}_{m}}{{P}_{m}},\end{eqnarray*}$$

where ${{P}_{m}}={{\psi }_{m}}\psi _{m}^{\dagger }$ is a projection operator with $P_{m}^{2}={{P}_{m}}$. For non-orthogonal MOs we have ${{P}_{m}}{{P}_{{{m}^{\prime }}}}\ne 0$ when the different MOs m and ${{m}^{\prime }}$ have a non zero overlap. Since P_ms project into a one-dimensional space, the dimension of the Hamiltonian h is at most M. We still express the Hamiltonian as an N×N matrix, but this should be redundant as the $(N-M)$-dimensional subspace has to be null, with h having $N-M$ zero eigenvalues. We can explicitly show this by writing the N×N matrix h as (see footnote 4 of [11])

where $\mathcal{E}=\ {\rm diag}({{\mathcal{E}}_{1}},\cdots ,{{\mathcal{E}}_{M}})$ is an M-dimensional diagonal matrix and Ψ(Ψ₁, ⋯, Ψ_M) is an N×M matrix composed of the molecular orbitals (ψ' s) as columns. Then, by simple algebra, the secular equation for h becomes

$$\begin{eqnarray*}&&\;{{{\rm det} }_{N}}\;(\lambda {{E}_{N}}-h)={{\lambda }^{N}}\;{{{\rm det} }_{N}}\;({{E}_{N}}-{{\lambda }^{-1}}\Psi \mathcal{E}{{\Psi }^{\dagger }})={{\lambda }^{N-M}}\;{{{\rm det} }_{M}}\;(\lambda {{E}_{M}}-\mathcal{E}{{\Psi }^{\dagger }}\Psi )=0,\end{eqnarray*}$$

where E_n is an n-dimensional unit matrix. Then one can see that h has $(N-M)$-fold degenerate zero eigenvalues $\lambda =0$. These zero-energy states are topological in that they are stable against continuous deformations of the parameters such as ${{\mathcal{E}}_{m}}$ and ${{\psi }_{im}}$. These zero states are stable as far as the number of the molecular orbitals is fixed. It is a finite dimensional analogue of Atiyah and Singer's index theorem [12].

Following this idea, we state the theorem in a slightly extended form, which we shall use later in the present paper. For the projection ${{P}_{m}}=P_{m}^{2}=P_{m}^{\dagger }$, let us define its dimension, ${\rm dim}\;{{P}_{m}}={\rm Tr}{{P}_{m}}={\rm rank}\;{{P}_{m}}$. Since the projection operator has eigenvalues 0 or 1, the dimension is a number of nonzero eigenvalues of P_m, which also coincides with the rank of the matrix represerntation of P_m. Then the number of zero modes, Z, should satisfy a condition,

$$\begin{eqnarray*}&&Z\geqslant N-\mathop{\sum }\limits_{m}{\rm dim}\;{{P}_{m}}.\end{eqnarray*}$$

We note here that the flat bands at the zero energy of the chiral symmetric models [9, 12–14] can be discussed with the argument presented above, where a square of the Hamiltonian is considered. Since the Hamiltonian in the chiral class is written as $\left[ \begin{array}{ccccccccccccccc} O & D \\ {{D}^{\dagger }} & O \\ \end{array} \right]$ in a suitable basis, its square is $\left[ \begin{array}{ccccccccccccccc} D{{D}^{\dagger }} & O \\ O & {{D}^{\dagger }}D \\ \end{array} \right]$. When D is an N×M matrix ($N\gt M$), one may identify $\Psi =D$ and $\mathcal{E}={{E}_{M}}$. The $(1,1)$ block of this squared Hamiltonian, $D{{D}^{\dagger }}$, corresponds to the present Hamiltonian. The theorem here then guarantees the existence of zero-energy states of the chiral symmetric Hamiltonian with degeneracy $N-M$. We can also see that the overlapping molecular orbitals in real space are discussed in the context of the rigorous treatment of ferromagnetism on the Hubbard model [13, 15, 16]. Non particle–hole symmetric flat bands on a special shape of the lattice (partial line graphs) are discussed as well [17]. The relation to the present analysis could therefore be an interesting problem, and should be discussed in the future.

Weaire and Thorpe considered a simple but multi-orbital tight-binding model for sp³ electrons on a diamond lattice [8], where the original motivation was to treat amorphous silicone. Let us reproduce the model here for later reference. We start with sp³-hybridized orbitals on a single tetrahedron. The local Hamiltonian for the tetrahedron reads

$$\begin{eqnarray*}&&{{H}_{sp3}}={{\epsilon }_{s}}c_{s}^{\dagger }{{c}_{s}}+{{\epsilon }_{p}}(c_{{{p}_{x}}}^{\dagger }{{c}_{{{p}_{x}}}}+c_{{{p}_{y}}}^{\dagger }{{c}_{{{p}_{y}}}}+c_{{{p}_{z}}}^{\dagger }{{c}_{{{p}_{z}}}})\equiv {{{\boldsymbol{c}} }^{\dagger }}{{h}_{s{{p}^{3}}}}{\boldsymbol{c}} ,\end{eqnarray*}$$

where ${{c}_{i}}\ \ (i=s,{{p}_{x}},{{p}_{y}},{{p}_{z}})$ is an annihilation operator of the bond orbitals with energy levels ${{\epsilon }_{s}}$ for the s orbital and ${{\epsilon }_{p}}$ for the p orbitals, and

$$\begin{eqnarray*}{{h}_{s{{p}^{3}}}}={{\epsilon }_{p}}{{E}_{4}}+{{V}_{1}}\left[ \begin{array}{ccccccccccccccc} 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 \\ \end{array} \right],{\boldsymbol{c}} =\left[ \begin{array}{ccccccccccccccc} {{c}_{0}} \\ {{c}_{1}} \\ {{c}_{2}} \\ {{c}_{3}} \\ \end{array} \right]=\frac{1}{2}\left[ \begin{array}{ccccccccccccccc} 1 & 1 & 1 & 1 \\ 1 & 1 & -1 & -1 \\ 1 & -1 & 1 & -1 \\ 1 & -1 & -1 & 1 \\ \end{array} \right]\left[ \begin{array}{ccccccccccccccc} {{c}_{s}} \\ {{c}_{{{p}_{x}}}} \\ {{c}_{{{p}_{y}}}} \\ {{c}_{{{p}_{z}}}} \\ \end{array} \right].\end{eqnarray*}$$

Here E₄ is the 4 × 4 unit matrix, and ${{V}_{1}}=\frac{1}{4}({{\epsilon }_{s}}-{{\epsilon }_{p}})$ is proportional to the s-p level offset. With this bond basis, the WT model for silicon atoms on a diamond lattice considers only the hopping, denoted by V₂, between bond-sharing orbitals.

For the diamond lattice (see figure 1(a)), the Hamiltonian in the Bloch picture of the WT model (H_k in the appendix B of [8], but note that here we take the origin of energy at ${{\epsilon }_{p}}$, so that the energy in [8] is shifted by $-{{V}_{1}}$ from ours) is written as

$$\begin{eqnarray*} H_{\rm WT}(\bm{k}) &=& \left[ \begin{array}{cccccccc} V_1 & V_1 & V_1 & V_1 & V_2 & 0& 0& 0 \\ V_1 & V_1 & V_1 & V_1 & 0& V_2 & 0& 0 \\ V_1 & V_1 & V_1 & V_1 & 0&0 &V_2 & 0 \\ V_1 & V_1 & V_1 & V_1 & 0& 0& 0& V_2 \\ V_2 & 0 & 0& 0 & V_1 & V_1 {\rm e}^{{\rm i}(k_y+k_z)} &V_1 {\rm e}^{{\rm i}(k_z+k_x)} & V_1{\rm e}^{{\rm i}(k_x+k_y)} \\ 0& V_2 & 0& 0 &V_1 {\rm e}^{{\rm -i}(k_y+k_z)} & V_1 &V_1 {\rm e}^{{\rm i}(k_x-k_y)} & V_1{\rm e}^{{\rm -i}(k_z-k_x)} \\ 0& 0& V_2 & 0 &V_1 {\rm e}^{{\rm -i}(k_z+k_x)} &V_1 {\rm e}^{{\rm -i}(k_x-k_y)} & V_1 & V_1{\rm e}^{{\rm i}(k_y-k_z)} \\ 0 & 0& 0& V_2 &V_1 {\rm e}^{{\rm -i}(k_x+k_y)} &V_1 {\rm e}^{{\rm i}(k_z-k_x)} &V_1 {\rm e}^{{\rm -i}(k_y-k_z)} & V_1 \end{array} \right] \\ & =& \left[ \begin{array}{cc} {H_V(0)} & {V_2 E_4}\\ {V_2 E_4} & {H_V(\bm{k})}, \end{array} \right] \\ H_V(\bm{k})&=& 4V_1 \psi_{\bm{k}} \psi_{\bm{k}} ^\dagger, \\ \psi_{\bm{k}} &=& \frac 1 2 \left[ \begin{array}{c} {{\rm e}^{{\rm i}(k_x+k_y+k_z)}}\\ {{\rm e}^{{\rm i} k_x}}\\ {{\rm e}^{{\rm i} k_y}}\\ {{\rm e}^{{\rm i} k_z}} \end{array} \right] = \frac 1 2 \left[ \begin{array}{c} {1}\\ {{\rm e}^{{\rm -i} k_1}}\\ {{\rm e}^{{\rm -i} k_2}}\\ {{\rm e}^{{\rm -i} k_3}} \end{array} \right] {\rm e}^{{\rm i} (k_x+k_y+k_z)}, \end{eqnarray*}$$

where ${{k}_{1}}={{k}_{y}}+{{k}_{z}}$, ${{k}_{2}}={{k}_{z}}+{{k}_{x}}$, ${{k}_{3}}={{k}_{x}}+{{k}_{y}}$, and ${{\psi }_{{\boldsymbol{k}} }}\psi _{{\boldsymbol{k}} }^{\dagger }\propto {{H}_{V}}$ is a projection onto the space spanned by ${{\psi }_{{\boldsymbol{k}} }}$. We can then apply the discussion in section 2 to note that the 4 × 4 matrix ${{H}_{V}}({\boldsymbol{k}} )$ has at most one nonzero energy that corresponds to a localized molecular orbital at each tetrahedron, while there are three zero-energy flat bands (in the present choice of the origin of energy).

Zoom In Zoom Out Reset image size

With this observation, we can introduce two representations for the 8 × 8 Hamiltonian as ${{H}_{{\rm WT}}}({\boldsymbol{k}} )\pm {{V}_{2}}{{E}_{8}}$, where E₈ is an 8 × 8 unit matrix. Although one might think the choice of the origin of the energy to be irrelevant, the purpose here is that we want to deal with flat bands that have nonzero energies. To do that, we can shift the origin of the energy to apply an algebraic argument. Now, if we take the plus sign, we have

$$\begin{eqnarray*}\begin{array}{rcl} {{H}_{{\rm WT}}}({\boldsymbol{k}} )+{{V}_{2}}{{E}_{8}} & = & \left[ \begin{array}{ccccccccccccccc} {{H}_{V}}(0) & 0 \\ 0 & {{H}_{V}}({\boldsymbol{k}} ) \\ \end{array} \right]+{{V}_{2}}\left[ \begin{array}{ccccccccccccccc} {{E}_{4}} & {{E}_{4}} \\ {{E}_{4}} & {{E}_{4}} \\ \end{array} \right] \\ {} & = & 4{{V}_{1}}{{P}_{1}}+4{{V}_{1}}{{P}_{2}}+2{{V}_{2}}{{P}_{3+}}, \\ \end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}\begin{array}{rcl} {{P}_{1}} & = & \left[ \begin{array}{ccccccccccccccc} {{\psi }_{0}}\psi _{0}^{\dagger } & O \\ O & O \\ \end{array} \right]={{\Psi }_{1}}\Psi _{1}^{\dagger },{{P}_{2}}=\left[ \begin{array}{ccccccccccccccc} O & O \\ O & {{\psi }_{{\boldsymbol{k}} }}\psi _{{\boldsymbol{k}} }^{\dagger } \\ \end{array} \right]={{\Psi }_{2}}\Psi _{2}^{\dagger },{{P}_{3+}}=\frac{1}{2}\left[ \begin{array}{ccccccccccccccc} {{E}_{4}} & {{E}_{4}} \\ {{E}_{4}} & {{E}_{4}} \\ \end{array} \right]={{\Psi }_{3+}}\Psi _{3+}^{\dagger }, \\ {{\Psi }_{1}} & = & \frac{1}{2}\left[ \begin{array}{ccccccccccccccc} 1 \\ 1 \\ 1 \\ 1 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{array} \right],\quad {{\Psi }_{2}}=\frac{1}{2}\left[ \begin{array}{ccccccccccccccc} 0 \\ 0 \\ 0 \\ 0 \\ 1 \\ {{{\rm e}}^{-{\rm i}{{k}_{1}}}} \\ {{{\rm e}}^{-{\rm i}{{k}_{2}}}} \\ {{{\rm e}}^{-{\rm i}{{k}_{3}}}} \\ \end{array} \right],\quad {{\Psi }_{3+}}=\frac{1}{\sqrt{2}}\left[ \begin{array}{ccccccccccccccc} {{E}_{4}} \\ {{E}_{4}} \\ \end{array} \right], \\ \end{array}\end{eqnarray*}$$

where P_is are projections with $P_{i}^{2}={{P}_{i}}$ and ${\rm dim}{{P}_{i}}={\rm Tr}\;{{\Psi }_{i}}\Psi _{i}^{\dagger }={\rm Tr}\;\Psi _{i}^{\dagger }{{\Psi }_{i}}=1,1,4$ respectively for $i=1,2,3+$, and $\Psi _{1}^{\dagger }{{\Psi }_{1}}=\Psi _{2}^{\dagger }{{\Psi }_{2}}=1,\Psi _{3+}^{\dagger }{{\Psi }_{3+}}={{E}_{4}}$. If we count the dimensions, ${{H}_{{\rm WT}}}({\boldsymbol{k}} )+{{V}_{2}}{{E}_{8}}$ has at most ${\rm dim}{{P}_{1}}+{\rm dim}{{P}_{2}}+{\rm dim}{{P}_{3+}}=1+1+4$ nonzero-energy bands, that is, there are $8-6=2$ zero-energy flat bands.

Similarly we have

$$\begin{eqnarray*}\begin{array}{rcl} {{H}_{{\rm WT}}}({\boldsymbol{k}} )-{{V}_{2}}{{E}_{8}} & = & \left[ \begin{array}{ccccccccccccccc} {{H}_{V}}(0) & 0 \\ 0 & {{H}_{V}}({\boldsymbol{k}} ) \\ \end{array} \right]-{{V}_{2}}\left[ \begin{array}{ccccccccccccccc} {{E}_{4}} & -{{E}_{4}} \\ -{{E}_{4}} & {{E}_{4}} \\ \end{array} \right] \\ {} & = & 4{{V}_{1}}{{P}_{1}}+4{{V}_{1}}{{P}_{2}}-2{{V}_{2}}{{P}_{3-}}, \\ \end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}&&{{P}_{3-}}={{\Psi }_{3-}}{{({{\Psi }_{3-}})}^{\dagger }},\end{eqnarray*}$$

$$\begin{eqnarray*}&&{{\Psi }_{3-}}=\frac{1}{\sqrt{2}}\left[ \begin{array}{ccccccccccccccc} -{{E}_{4}} \\ {{E}_{4}} \\ \end{array} \right],\end{eqnarray*}$$

where ${{P}_{3-}}$ is a projection similar to ${{P}_{3+}}$. Now we see that ${{H}_{{\rm WT}}}-{{V}_{2}}{{E}_{8}}$ has at most $1+1+4$ nonzero-energy bands, that is, there are again two zero-energy flat bands.

Hence we end up with two flat bands at each $E=\pm {{V}_{2}}$ in the WT Hamiltonian. To compare with [8], we need to shift the energy, after which the flat bands are at $-{{V}_{1}}\pm {{V}_{2}}$, since our choice of ${{H}_{{\rm WT}}}$ is shifted by $-{{V}_{1}}{{E}_{8}}$ from the Hamiltonian in [8]. An essential point is that we have succeeded in describing the flat bands at nonzero energies algebraically in terms of the overlapping molecular orbitals, exemplified here in the WT model.

Now we come to the original aim of describing silicene. We start with an observation that 2D silicene can also be captured in a manner similar to the WT model in 3D. As indicated in figure 1(b), we have three primitive vectors, ${{e}_{1}},{{e}_{2}},{{e}_{3}}$, from which we have three reciprocal vectors for the 3D diamond lattice. Corresponding 2D momentum components are given by $({{k}_{1}},{{k}_{2}})$ with ${{k}_{3}}=0$. The Hamiltonian can be obtained from that in 3D by cutting the bonds at the blue bonds in figure 1. Then the Hamiltonian, as a simple extension of the WT model for silicene, can be taken as

$$\begin{eqnarray*}H_{{\rm Silicene}}^{0}=\left[ \begin{array}{ccccccccccccccc} {{H}_{V}}(0) & {{V}_{2}}E_{4}^{C} \\ {{V}_{2}}E_{4}^{C} & {{H}_{V}}({\boldsymbol{k}} ) \\ \end{array} \right],\quad E_{4}^{C}=\left[ \begin{array}{ccccccccccccccc} {{E}_{3}} & {} \\ {} & 0 \\ \end{array} \right]={{E}_{4}}-\mathcal{E},\quad \mathcal{E}={\rm diag}(0,0,0,1).\end{eqnarray*}$$

It is important to note that the bonds normal to the two-dimensional plane, which are originally dangling bonds after the dissection that can be treated with hydrogen termination for instance, are different from the other bonds. A simple way of implementing this is to modify the Hamiltonian into

$$\begin{eqnarray*}H_{{\rm Silicene}}^{{{\epsilon }_{H}}}({\boldsymbol{k}} )=\left[ \begin{array}{ccccccccccccccc} {{H}_{V}}(0)-{{\epsilon }_{H}}\mathcal{E} & {{V}_{2}}E_{4}^{C} \\ {{V}_{2}}E_{4}^{C} & {{H}_{V}}({\boldsymbol{k}} )-{{\epsilon }_{H}}\mathcal{E} \\ \end{array} \right],\end{eqnarray*}$$

where the energy ${{\epsilon }_{H}}$ can be controlled by how the bond is chemically terminated.

Now, an interesting observation is that one can precisely apply the discussion given in section 2 to show that the model has flat bands. This follows from a simple observation that

$$\begin{eqnarray*}\begin{array}{rcl} H_{{\rm Silicene}}^{{{\epsilon }_{H}}}+{{V}_{2}}{{E}_{8}} & = & 4{{V}_{1}}{{P}_{1}}+4{{V}_{1}}{{P}_{2}}+{{V}_{2}}\left[ \begin{array}{ccccccccccccccc} E_{4}^{C} & E_{4}^{C} \\ E_{4}^{C} & E_{4}^{C} \\ \end{array} \right]+({{V}_{2}}-{{\epsilon }_{H}})\left[ \begin{array}{ccccccccccccccc} \mathcal{E} & O \\ O & \mathcal{E} \\ \end{array} \right] \\ {} & = & 4{{V}_{1}}{{P}_{1}}+4{{V}_{1}}{{P}_{2}}+2{{V}_{2}}P_{3+}^{C}+({{V}_{2}}-{{\epsilon }_{H}}){{P}_{5}}, \\ \end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}\begin{array}{rcl} H_{{\rm Silicene}}^{{{\epsilon }_{H}}}-{{V}_{2}}{{E}_{8}} & = & 4{{V}_{1}}{{P}_{1}}+4{{V}_{1}}{{P}_{2}}-{{V}_{2}}\left[ \begin{array}{ccccccccccccccc} E_{4}^{C} & -E_{4}^{C} \\ -E_{4}^{C} & E_{4}^{C} \\ \end{array} \right]-({{V}_{2}}+{{\epsilon }_{H}})\left[ \begin{array}{ccccccccccccccc} \mathcal{E} & O \\ O & \mathcal{E} \\ \end{array} \right] \\ {} & = & 4{{V}_{1}}{{P}_{1}}+4{{V}_{1}}{{P}_{2}}-2{{V}_{2}}P_{3-}^{C}-({{V}_{2}}+{{\epsilon }_{H}}){{P}_{5}}, \\ \end{array}\end{eqnarray*}$$

where

$$\begin{eqnarray*}\begin{array}{l} P_{3\pm }^{C}=\frac{1}{2}\left[ \begin{array}{ccccccccccccccc} E_{4}^{C} & \pm E_{4}^{C} \\ \pm E_{4}^{C} & E_{4}^{C} \\ \end{array} \right]=\Psi _{3\pm }^{C}{{(\Psi _{3\pm }^{C})}^{\dagger }}, \\ \Psi _{3\pm }^{C}=\frac{1}{\sqrt{2}}\left[ \begin{array}{ccccccccccccccc} \pm E_{4}^{C} \\ E_{4}^{C} \\ \end{array} \right], \\ {{(\Psi _{3\pm }^{C})}^{\dagger }}\Psi _{3\pm }^{C}=E_{4}^{C}={{(E_{4}^{C})}^{2}},\ {\rm Tr}E_{4}^{C}=3, \\ {{P}_{5}}={{\Psi }_{5}}\Psi _{5}^{\dagger }=\left[ \begin{array}{ccccccccccccccc} \mathcal{E} & O \\ O & \mathcal{E} \\ \end{array} \right]=P_{5}^{2}, \\ {{\Psi }_{5}}=({{\Psi }_{5,1}},{{\Psi }_{5,2}}),{{\Psi }_{5,1}}=\left[ \begin{array}{ccccccccccccccc} \mathcal{E} \\ O \\ \end{array} \right],{{\Psi }_{5,2}}=\left[ \begin{array}{ccccccccccccccc} O \\ \mathcal{E} \\ \end{array} \right], \\ {{\Psi }_{5,1}}\Psi _{5,1}^{\dagger }=\left[ \begin{array}{ccccccccccccccc} \mathcal{E} & O \\ O & O \\ \end{array} \right],\ {{\Psi }_{5,2}}\Psi _{5,2}^{\dagger }=\left[ \begin{array}{ccccccccccccccc} O & O \\ O & \mathcal{E} \\ \end{array} \right],\ \\ \Psi _{5,i}^{\dagger }{{\Psi }_{5,i}}={{\mathcal{E}}^{2}}=\mathcal{E},(i=1,2),\ \Psi _{5,1}^{\dagger }{{\Psi }_{5,2}}=\Psi _{5,2}^{\dagger }{{\Psi }_{5,1}}=O, \\ \Psi _{5}^{\dagger }{{\Psi }_{5}}=\left[ \begin{array}{ccccccccccccccc} \Psi _{5,1}^{\dagger }{{\Psi }_{5,1}} & \Psi _{5,1}^{\dagger }{{\Psi }_{5,2}} \\ \Psi _{5,2}^{\dagger }{{\Psi }_{5,1}} & \Psi _{5,2}^{\dagger }{{\Psi }_{5,2}} \\ \end{array} \right]=\left[ \begin{array}{ccccccccccccccc} \mathcal{E} & O \\ O & \mathcal{E} \\ \end{array} \right]. \\ \end{array}\end{eqnarray*}$$

The dimensions of the projections are evaluated as

$$\begin{eqnarray*}\begin{array}{rcl} {\rm dim}P_{3\pm }^{C} & = & {\rm Tr}\Psi _{3\pm }^{C}\Psi {{_{3\pm }^{C}}^{\dagger }}={\rm Tr}\Psi {{_{3\pm }^{C}}^{\dagger }}\Psi _{3\pm }^{C}={\rm Tr}E_{4}^{C}=3, \\ {\rm dim}{{P}_{5}} & = & {\rm Tr}{{P}_{5}}={\rm Tr}{{\Psi }_{5}}\Psi _{5}^{\dagger }={\rm Tr}\Psi _{5}^{\dagger }{{\Psi }_{5}}=2{\rm Tr}\mathcal{E}=2. \\ \end{array}\end{eqnarray*}$$

Since we have now expressed the Hamiltonian as a linear combination of projection operators, we can see that the wavefunctions associated with P₁ and P₂ are localized within each tetrahedron, the ones with P₅ are localized within the bonds perpendicular to the plane while the ones with $P_{3\pm }^{C}$ extend over the 2D plane. Then, counting the dimensions tells us that $H_{{\rm Silicene}}^{{{\epsilon }_{H}}}\pm {{V}_{2}}{{E}_{8}}$ has at most ${\rm dim}{{P}_{1}}+{\rm dim}{{P}_{2}}+{\rm dim}{{P}_{3\pm }}+{\rm dim}{{P}_{5}}=1+1+3+2=7$ nonzero-energy bands, that is, there is a 8−7 = 1 flat band at ±V₂ generically.

At a special point of ${{\epsilon }_{H}}=\pm {{V}_{2}}$, the coefficient of the projection P₅ for the expansion of $H_{{\rm Silicene}}^{{{\epsilon }_{H}}=\pm {{V}_{2}}}\pm {{V}_{2}}{{E}_{8}}$ happens to vanish. The additional two dimensional space of the Hamiltonian then becomes null, which makes the flat bands at energies V₂ three-fold degenerate.

As for the signs of the parameters, one has ${{V}_{1}}\lt 0$ since ${{\epsilon }_{s}}\lt {{\epsilon }_{p}}$, and ${{V}_{2}}\lt 0$ since the hopping gains energy. Further, it is natural to assume the bonds normal to the plane are close to being dangling bonds, that is ${{\epsilon }_{H}}\lt 0$ for free-standing silicene. Then we may consider the case where ${{\epsilon }_{H}}={{V}_{2}}$ belongs to a regime that contains free-standing silicene. Thus we can see, algebraically, that the difference in the structure of the Hamiltonian produces the following: while in 3D the dispersive bands with Dirac cones are sandwiched by doubly-degenerate flat bands, the band structure in 2D silicene has dispersive bands with cones sandwiched by triply-degenerate flat bands and a non-degenerate flat band.

Now let us demonstrate the above analytic formulation from numerical results for the band structure in figure 2. Panel (a) depicts the ideal case of ${{\epsilon }_{H}}={{V}_{2}}=-1.0$, where we can indeed see the triply-degenerate flat bands at $E=-{{V}_{2}}=1$ along with a non-degenerate flat band at $E={{V}_{2}}=-1$. The result confirms the analytic discussion above. Panel (b) depicts a general case of ${{\epsilon }_{H}}\lt {{V}_{2}}=-1$, where we can see that the two out of the three flat bands become somewhat dispersive around $E=-{{\epsilon }_{H}}$ while one flat band remains flat at $E=-{{V}_{2}}=1$. The two nearly-flat bands are derived from the dangling bonds perpendicular to the plane. Interestingly, if we compare the result for the dispersion of silicene obtained by Takeda and Shiraishi [4] as reproduced in figure 3 here, we can see that they roughly agree with each other in terms of the bands' widths, ordering and multiplicity. Specifically, we can make the following observation. In both the realistic band structure [4] and the present algebraic treatment for silicene, we have esentially two bands that contain Dirac dispersions. The bands are sandwiched by two narrow (ideally flat) bands, each of which touches the dispersive ones at the Γ point. Above these, there are two nearly-flat bands [4]. Thus the three narrow bands above the Dirac cone can be traced back to the three-fold degenerate flat bands in the idealized model. Of course, this is only a rough mapping, but these qualitative agreements in the band structure between the WT-like model and the realistic band calculation gives an insight into a physical origin of the electronic band structure of silicene. The existence of the (nearly) flat bands has an intuitive origin as well. If electrons hop like molecular orbitals, this imposes a strong constraint on the description in terms of the original electron. Namely, the degrees of freedom other than the molecular-orbital hopping remain more or less frozen, which form the nearly-flat (or less dispersive) bands.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Within the present simple model, one may expect that the electrons at the perpendicular bonds can be stabilized by hydrogen termination. It can be modeled as ${{\epsilon }_{H}}=-{{V}_{2}}\gt 0$. Then the flat bands originally situated at $-{{V}_{2}}\gt 0$ will move down to V₂, that is, below the energy region of the Dirac cones. Then the Fermi energy, which is originally situated at the Dirac point since the sp³-bands are half-filled in silicene, will move away from the original Dirac point.

4.1. Dirac cones and effects of buckling

We can even extend the present argument when there is some buckling in the structure. For this we start with an observation that existence of the Dirac cones in silicene is due to the time-reversal (T) and crystal symmetry (e.g. reflection, R, that exchanges bonds 1 and 2) [1, 18, 19]. This is readily seen in the present model, where the momentum dependence is expressed in terms of ${{\psi }_{{\boldsymbol{k}} }}$. At K and K $^{\prime }$ points, ${{\psi }_{{\boldsymbol{k}} }}$ is expressed and transformed as

$$\begin{eqnarray*}&&{{\psi }_{{\boldsymbol{k}} }}=\frac{1}{2}\left[ \begin{array}{ccccccccccccccc} 1 \\ \omega \\ {{\omega }^{2}} \\ 1 \\ \end{array} \right],\quad TR{{\psi }_{{\boldsymbol{k}} }}=\frac{1}{2}\;T\left[ \begin{array}{ccccccccccccccc} 1 \\ {{\omega }^{2}} \\ \omega \\ 1 \\ \end{array} \right]={{\psi }_{{\boldsymbol{k}} }},\end{eqnarray*}$$

with $\omega \ne 1$ being a nontrivial cubic root of unity, which is the origin of the degeneracy at K and K $^{\prime }$ (denoted by S points in [4]). Since the degeneracy is lifted at the general momentum, the energy gap is linearly vanishing (Dirac cones) around K and ${{K}^{\prime }}$. The property is stable as far as the symmetries T and R remain.

Then what are the effects of the buckling in the sp³ structure observed in silicene? Since the main effect of the buckling is a modification of the bond-to-bond angles among the sp³ orbitals, the local hamiltonian within a tetrahedron is roughly estimated to depend on the buckling angle θ as

$$\begin{eqnarray*}\begin{array}{rcl} {{h}_{{\rm local}}}(\theta ) & = & \mathop{\sum }\limits_{\langle i,j\rangle }{{V}_{ij}}c_{i}^{\dagger }{{c}_{j}}+{\rm h}.{\rm c}., \\ {{V}_{ij}} & = & \left\{ \begin{array}{ccccccccccccccc} {{V}_{1}} & (\langle i,j\rangle =\langle 01\rangle ,\langle 12\rangle ,\langle 20\rangle ) \\ V_{1}^{\prime } & (\langle i,j\rangle =\langle 03\rangle ,\langle 13\rangle ,\langle 23\rangle ) \\ \end{array} \right. \\ \frac{V_{1}^{\prime }}{{{V}_{1}}} & = & \frac{{\rm cos} \theta }{{\rm cos} {{\theta }_{0}}},{\rm cos} {{\theta }_{0}}=-\frac{1}{3}, \\ \end{array}\end{eqnarray*}$$

where c₃ is the annihilation operator for the bond that is perpendicular to the silicene plane. Since the two bonds coupling the neighboring tetrahedra are remain straight even with the buckling, the inter-site term V₂ is not modified. As far as the buckling is small, the effects can be considered by a Hamiltonian,

$$\begin{eqnarray*}\begin{array}{rcl} H_{{\rm Silicene}}^{{{\epsilon }_{H}},{\rm cos} \theta }({\boldsymbol{k}} ) & = & \left[ \begin{array}{ccccccccccccccc} H_{V}^{\theta }(0)-\epsilon _{H}^{\theta }\mathcal{E} & {{V}_{2}}E_{4}^{C} \\ {{V}_{2}}E_{4}^{C} & H_{V}^{\theta }({\boldsymbol{k}} )-\epsilon _{H}^{\theta }\mathcal{E} \\ \end{array} \right], \\ H_{V}^{\theta }({\boldsymbol{k}} ) & = & 4{{V}_{1}}\psi _{{\boldsymbol{k}} }^{\theta }{{(\psi _{{\boldsymbol{k}} }^{\theta })}^{\dagger }}, \\ \psi _{{\boldsymbol{k}} }^{\theta } & = & \ {\rm diag}\left( 1,1,1,\frac{{\rm cos} \theta }{{\rm cos} {{\theta }_{0}}} \right){{\psi }_{{\boldsymbol{k}} }}, \\ \epsilon _{H}^{\theta } & = & {{\left( \frac{{\rm cos} \theta }{{\rm cos} {{\theta }_{0}}} \right)}^{2}}{{\epsilon }_{H}}. \\ \end{array}\end{eqnarray*}$$

This implies the Dirac cones and the flat bands we have focused are stable against the buckling. This is also regarded as an aspect of the topological stability. When the buckling becomes large enough, however, sp²-character of the hybridized orbitals will generate direct hoppings among the out-of-plane bonds, which is not included in the present formalism.

In the present paper, after introducing a generic argument for flat bands, the Weaire–Thorpe model, originally conceived for 3D silicon, is extended to 2D silicene. A surprise revealed here is that the flat bands that arise in the WT model for the hybridized sp³ orbitals in 3D also appear, in ideal situations, in silicene, but with different degeneracies in the flat bands for an algebraic reason. In this picture, the band structure, including the flat ones, is crucially controlled by the out-of-atomic-plane orbits. We have further pointed out there are pairs of bands, each of which contains a Dirac cone.

The flat bands emerging in our treatment can be theoretically interesting and important as a multi-orbital effect⁵ . Finite samples will accommodate characteristic edge states. The flat bands in the idealized model, which should become dispersive in realistic situations, will still have a large density of states. If we can shift the chemical potential close to the flat bands, e.g., by chemical doping, interesting phenomena are expected. Among these are (i) structural instabilities, such as those observed experimentally, and (ii) flat-band ferromagnetism [20]. We believe that these will aid the synthesis and characterization of silicene.

This work was supported in part by JSPS grant numbers 26247064(YH,HA), 25610101(YH) and 25107005 from MEXT(YH,HA).