Group theoretical and topological analysis of the quantum spin Hall effect in silicene

Silicene is a monolayer of silicon atoms arranged in a buckled honeycomb lattice (see figure 1 for a schematic). In contrast to graphene, the two basis atoms of the unit cell (called A and B) are separated perpendicular to the atomic plane at a distance 2l with l = 0.23 Å [14]. As there is no translation symmetry in the out-of-plane direction, the material is quasi-two-dimensional. The buckling is quantified by an angle θ ⩾ 90° as shown in figure 1.

Zoom In Zoom Out Reset image size

In figure 2, we provide an illustration of the lattice of silicene—in real and reciprocal space. Basis vectors defining the unit cell are given by

$$\begin{equation} \vec{a}_1=a \left(\frac{\sqrt{3}}{2},\frac{3}{2},0\right), \quad \vec{a}_2=a \left(-\frac{\sqrt{3}}{2},\frac{3}{2},0\right) \end{equation} \tag{ 1 }$$

in real space and by

$$\begin{equation} \vec{b}_1=\frac{2\pi}{a} \left(\frac{1}{\sqrt{3}},\frac{1}{3}\right), \quad \vec{b}_2=\frac{2\pi}{a} \left(-\frac{1}{\sqrt{3}},\frac{1}{3}\right) \end{equation} \tag{ 2 }$$

in reciprocal space. Here, a is the distance between two neighboring silicon atoms. In real space without buckling, the lattice exhibits D_6h symmetry (i.e. the graphene case). All the symmetry operations sketched in figure 2 are present, as well as bulk inversion i and reflection at the atomic plane itself, σ_h. With the buckling present (silicene), symmetries C₆, C₂, C₂', σ'' and σ_h are broken and we are left with point group D_3d.

Zoom In Zoom Out Reset image size

Let us now go to reciprocal space. For symmorphic groups, the Γ point always has the same symmetry as the real lattice. However, the notation in use of symmetry axes with a single prime referred to the fact, that these axes cross the corners of the underlying hexagon (see the caption of figure 2). Under Fourier transformation, the hexagonal lattice is rotated by 90° with respect to the axes of a fixed coordinate system. The positions of the C₂- and σ-axes stay the same then, while the corners of the hexagon come to rest on the previously double-primed symmetry axes now. To have a corresponding group theoretical notation in reciprocal space as well, we rename the symmetry axes. Fourier transformation of the lattice can thus effectively be considered as an interchange of single and double-primed operations (see figure 2).

In particular, the D_3d-symmetry at the Γ-point is equivalent to operations C₃, C₂', σ'', i and S₆ in reciprocal space, which corresponds to the point group D₃ with the additional symmetry classes {σ'',i,S₆}. At the K points of the BZ, the symmetry of the group of the wave vector is further reduced to the point group D₃.

The full knowledge of the lattice symmetries makes it in principle possible to construct a low-energy Hamiltonian by expansion around high-symmetry points. This can be done with powerful and well-established approaches, for instance, the invariant expansion [18–20]. In undoped silicene, the Fermi level is located at the K-points of the BZ. Hence, a low-energy Hamiltonian constructed by a symmetry analysis near the K-points will capture the essential physical properties of the system.

We perform an invariant expansion around the K-points of silicene, which were identified to exhibit the symmetry point group D₃. The π-orbital wavefunction transforms like the two-dimensional irreducible representations (IR) Γ₃ of the group D₃, while the spin part is represented by the IR Γ₄ of the double group. Therefore, our total wavefunction is of the form of the product Γ₃ × Γ*₄ = Γ₄ + Γ₅ and the starting point for the derivation of the Hamiltonian. A detailed presentation of this expansion is given in appendix A.

Consequently, the low-energy Hamiltonian of silicene near the K-points is found to be (in lowest orders of $\vec {k}$ and E_z)

$$\begin{eqnarray} &&\fl \mathcal{H}^{K(K')} = a_1\mathbb{I}+ a_2 \tau_z \sigma_z s_z + a_3 \sigma_z s_0 E_z + a_4 \tau_z \sigma_0 s_z E_z + a_5 ( \tau_z \sigma_x s_y - \sigma_y s_x) E_z+ a_7 (\tau_z \sigma_x k_x +\sigma_y k_y) s_0 \nonumber \\ &&+ a_{10} E_z [\sigma_x(s_x k_y + s_y k_x)+\tau_z\sigma_y(s_x k_x -s_y k_y)] + a_{11} E_z \sigma_0 (s_x k_y - s_y k_x) \nonumber \\ &&+ a_{13}\sigma_z (s_x k_y - s_y k_x)+ a_{16}E_z(\sigma_x k_x +\tau_z\sigma_y k_y)s_z \end{eqnarray} \tag{ 3 }$$

presented in the basis (ψ_Aβ,ψ_Aα,ψ_Bβ,ψ_Bα)^T, where α = |↑〉, β = |↓〉. Here, τ_z = ± 1 distinguishes the inequivalent valleys K and K'. The Pauli matrices σ act on sublattices A and B, while s are the corresponding matrices in spin space. Additional constraints due to time-reversal symmetry (TRS) have already been taken into account, as explained in appendix A.

In equation (3), the terms proportional to a₂, a₅ and a₇ are well known from the graphene literature to describe its fundamental properties near the K-points [4]. Further corrections proportional to a₁₀ and a₁₁ can as well be found in graphene [21]. Beyond this, silicene exhibits specific terms proportional to a₃ and a₁₃ as reported in [12, 13]. Additionally, we find an electric field induced SOI-term proportional to a₄, that has not been reported for silicene before⁵.

Moreover, the Hamiltonian exhibits higher order corrections to the linear dispersion (proportional to a₁₆). We conclude, that the low-energy Hamiltonian of silicene near the K-points contains interesting SOI terms that are absent in graphene because of the difference in the symmetry of the two honeycomb lattices.

To supplement the symmetry analysis, a corresponding tight-binding model is discussed next. This model enables us to construct a valid Hamiltonian for the full BZ which is crucial for the subsequent topological analysis.

Let us start with a brief discussion of the internal spin–orbit coupling and subsequently introduce the other important terms of the tight-binding model. In real space, silicene is described in terms of a lattice model including the (standard) spin–orbit coupling term [4, 12]

$$\begin{equation} H_{\rm so}=\frac{\hbar}{4m_0c^2}(\vec{\nabla}V\times\vec{p})\cdot \vec{s}=-\frac{\hbar}{4m_0c^2}(\skew5\vec{F}\times\vec{p})\cdot \vec{s}, \end{equation} \tag{ 4 }$$

where $\skew5\vec{F}$ is the force stemming from the electric potential V , $\skew3\vec {p}$ is the momentum and $\vec {s}$ the spin of the electron. The (internal) electric force in silicene is provided by the crystal field in in-plane and out-of-plane direction. Since the momentum operator is oriented along nearest neighbor or next-nearest neighbor bonds, the silicene lattice forbids terms involving a crystal in-plane force coupling to a nearest neighbor bond. This consideration leads to the following spin–orbit Hamiltonian [12]:

$$\begin{eqnarray} &&H_{\rm so}=\imath \lambda_{\rm so} \sum_{\langle\langle i,j\rangle\rangle; \alpha\beta} v_{ij} c_{i\alpha}^\dagger s_z^{\alpha\beta} c_{j\beta}-\imath \frac{2}{3}\lambda_{r,2} \sum_{\langle\langle i,j\rangle\rangle; \alpha\beta} \mu_{i} c_{i\alpha}^\dagger (\vec{s}\times \skew5\hat{d}_{ij})_z^{\alpha\beta} c_{j\beta}\end{eqnarray} \tag{ 5 }$$

with, so far, undefined parameters λ_so and λ_r,2. In the latter equation, α and β are spin quantum numbers; the indexes i,j label the atomic site/orbital. Here and in the following, 〈i,j〉 denotes nearest neighbors and 〈〈i,j〉〉 next-nearest neighbors. The neighboring sites are each time connected by the vector $\vec {d}_{ij}$ with its corresponding unit vector $\skew5\hat {d}_{ij}$. The sign v_ij = ± 1 refers to the next-nearest neighbor hopping being anticlockwise or clockwise with respect to the positive z-axis. Furthermore, μ_i = ± 1 is introduced to distinguish between the A(B) site.

Additionally, there is a regular nearest-neighbor hopping term and on-site energies of the form

$$\begin{equation} H_0=\sum_{i\alpha} \epsilon_i c_{i\alpha}^\dagger c_{i\alpha}- \sum_{\langle i,j\rangle; \alpha\beta} t_{ij} c_{i\alpha}^\dagger c_{j\beta}, \end{equation} \tag{ 6 }$$

where _i is the on-site energy of the atomic orbital and t_ij the hopping parameter for hopping between the orbitals i and j of neighboring atomic sites. When an external electric field E_z is applied perpendicularly to the atomic plane of silicene, a staggered sublattice potential of the form

$$\begin{equation} H_{E}=\lambda_{\rm e} \sum_{i\alpha} E_z^i \mu_i c_{i\alpha}^\dagger c_{i\alpha} \end{equation} \tag{ 7 }$$

is generated [14] with a, so far, undefined parameter λ_e. Interestingly, Eⁱ_z allows for on-site transitions between the p_z and s orbitals [23].

To describe a full next-nearest neighbor tight-binding model, we also introduce spin–orbit terms involving external electric forces in equation (4), which we may index as H_R here due to their resemblance to Rashba terms. The simplest ones are

$$\begin{eqnarray} &&\fl H_{\rm R} =\imath \lambda_{r,1} \sum_{\langle i,j\rangle; \alpha\beta} c_{i\alpha}^\dagger (\vec{s}\times\hat{d}_{ij})_z^{\alpha\beta} c_{j\beta} E_z^j +\imath \lambda_{{\rm e},2} \sum_{\langle\langle i,j\rangle\rangle; \alpha\beta} v_{ij} c_{i\alpha}^\dagger s_z^{\alpha\beta} c_{j\beta} \mu_{i} E_z^j \nonumber \\ &&+ \imath \lambda_{r,3} \sum_{\langle\langle i,j\rangle\rangle; \alpha\beta} c_{i\alpha}^\dagger (\vec{s}\times\hat{d}_{ij})_z^{\alpha\beta} c_{j\beta} E_z^j. \end{eqnarray} \tag{ 8 }$$

Interestingly, the second term proportional to λ_e,2 is a multiplicative combination of the spin–orbit term and the staggered sublattice potential. This term corresponds to the additional, electric field induced spin–orbit term that we have found in the invariant expansion model (proportional to a₄). The full tight-binding Hamiltonian is then given by

$$\begin{equation} H=H_0+H_{\rm so}+H_E+H_{\rm R}. \end{equation} \tag{ 9 }$$

The terms given above will reproduce our Hamiltonian derived by symmetry analysis, equation (3), when being expanded around the K-points.

To be able to estimate the coupling constants introduced above, we now briefly discuss and extend a tight-binding model presented in [12] including π and σ-bands in silicene using the basis

$$\begin{equation} \{| p_z^A \rangle, |p_z^B \rangle, |p_y^A \rangle, |p_x^A \rangle, |s^A\rangle, |p_y^B \rangle, |p_x^B \rangle, |s^B\rangle \}\times\{\uparrow, \downarrow \}. \end{equation} \tag{ 10 }$$

Nearest-neighbor hopping matrix elements between given orbitals can then be expressed by parameters V₁, V₂ and V₃ as functions of Slater bond parameters V_ppπ, V_ppσ, V_spσ and curvature angle θ with

$$\begin{eqnarray*} &\displaystyle V_1=\frac{3}{4} \sin^2 \theta (V_{{\rm pp}\pi}-V_{{\rm pp}\sigma}),&\\ &\displaystyle V_2=\frac{3}{2} \sin\theta V_{{\rm sp}\sigma}, &\\ &\displaystyle V_3=\frac{3}{2} \sin\theta \cos\theta (V_{{\rm pp}\pi}-V_{{\rm pp}\sigma})& \end{eqnarray*}$$

(see [12] for details of the modeling). Beyond the analysis done in [12], we consider the influence of an electric field E_z applied perpendicularly to the atomic plane. A staggered sublattice-Hamiltonian is then induced that takes in the basis (10) the following form:

$$\begin{equation} H_E= e E_z \left( \begin{array}{@{}cccccccc@{}} 0 & 0 & 0 & 0 & z_0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & z_0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ z_0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & z_0 & 0 & 0 & 0 & 0 & 0 & 0 \end{array}\right) \times \mathbb{I}_{(2\times2)}\end{equation} \tag{ 11 }$$

as only transitions between p_z and s orbitals of the same site are allowed [23]. In this equation, e is the charge of an electron and z₀ the Stark element weighting the transition. Spin–orbit coupling is now included in the same way as in [12] by the term $H_{\mathrm { SO}}=\Delta \vec {L}\cdot \vec {s}$ where $\vec {L}$ is the angular momentum vector, $\vec {s}$ the spin operator and Δ the coupling constant. The tight-binding model of π and σ bands allows us to estimate some of the parameters introduced in the previous section as a function of the parameters V₁₋₃ (that are in principle known) as well as the Stark element z₀. To do so, we now apply a unitary transformation U to the Hamiltonian that corresponds to the following change of basis:

$$\begin{eqnarray} &&\fl U^{\mathrm{T}} \{ p_z^A, p_z^B, p_y^A, p_x^A, s^A,p_y^B, p_x^B, s^B \}\times\{\uparrow, \downarrow \}=\{ \phi_1, \phi_4, \phi_2, \phi_5, \phi_3, \phi_6, \phi_7, \phi_8 \}\times\{\uparrow, \downarrow \} \end{eqnarray} \tag{ 12 }$$

with mixed orbitals, for example, $\phi _1=u_{11} p_z^A + u_{21} s^A+ u_{31}(\frac {1}{\sqrt {2}}(p_x^B-\imath p_y^B))$. Then, the transformed Hamiltonian

$$\begin{equation} H'= U^\dagger (H_{\rm SO}+H_E) U \end{equation} \tag{ 13 }$$

can be analyzed perturbatively on the first (4 × 4) block corresponding to the basis {ϕ₁↑,ϕ₁↓,ϕ₄↑,ϕ₄↑}, which is expected to describe energy eigenstates near the Fermi energy [12]. Listing only terms including the external electric field, we find in first and second order perturbation theory the following three terms:

$$\begin{eqnarray*} &&H_{}^{(1,2)}= \lambda_{\rm e} \sigma_z s_0+ \lambda_{{\rm e},2} \sigma_0 s_z+\lambda_{r,1} (\sigma_x s_y + \sigma_y s_x) \end{eqnarray*}$$

with parameters

$$\begin{eqnarray} &&\displaystyle \lambda_{\rm e}= 2 e E_z z_0 u_{11} u_{21} \approx \frac{2 e E_z z_0 V_3}{V_2},\nonumber\\ &&\displaystyle \lambda_{{\rm e},2} \approx -\frac{\Delta d e E_{z} z_{0} V_{3} }{2\,V_{2}^3},\\ &&\displaystyle \lambda_{r,1} \approx \frac{\Delta e E_{z} z_{0}}{2 V_{2}}.\nonumber \end{eqnarray} \tag{ 14 }$$

In the last step, the expressions were approximated under the assumption of low buckling, meaning θ → 90° and V₃ → 0, where only the term of lowest order in V₃ was kept. The term proportional to λ_e represents the on-site hopping occurring only in a buckled structure. With the one proportional to λ_e,2, an electric field induced term of first order in SOI Δ is generated that is related to a next-nearest neighbor hopping. It depends linearly on V₃, so this term is specific to the buckled silicene structure as well. Finally, the term proportional to λ_r,1 is the well-known Rashba spin–orbit coupling reported already in [4, 24]. All terms given above agree nicely with our previous group theoretical analysis.

We now estimate the size of the spin–orbit terms coupling to an external electric field in silicene. The bond parameters V₁₋₃, the energy d, the angle θ = 101.7°, the lattice constant a, as well as the SOI Δ are adapted from [12]. We then approximate the Stark-element z₀ = 〈ϕ_n,0,0|z|ϕ_n,1,0〉 as transition matrix elements between atomic wave functions, where in silicene we have n = 3, since the outer shells are provided by the 3s and 3p orbitals. The corresponding wave functions were chosen as the wave functions of the hydrogen atom with the same quantum numbers. We then find $z_0=3\sqrt {6}\times a_{\mathrm { B}}/Z_{\mathrm { eff}}$ in silicene, where a_B is the Bohr radius and an outer electron experiences an effective atomic charge of Z_eff ≈ 4.29 [25] due to screening. Hence, we estimate z₀ = 0.906 Å. For typical values of E_z = 50 eV /300 nm, we calculate

$$\begin{eqnarray} &\displaystyle \lambda_{\rm e} \approx 8.5\,{\rm meV}, \nonumber\\ & \displaystyle\lambda_{{\rm e},2} \approx 12.8\, \mu {\rm eV},\\ & \displaystyle\lambda_{r,1} \approx 22.7\, \mu {\rm eV}. \nonumber \end{eqnarray} \tag{ 15 }$$

Thus, the additional term proportional to λ_e,2 is small compared to the term proportional to λ_e, but similar in magnitude as the Kane–Mele–Rashba term proportional to λ_r,1.

Let us briefly give an idea of how the topology of silicene may depend on an external electric field, orienting ourselves along the lines of [14]. We consider again the effective Hamiltonian of equation (3), derived by analytical expansion around the high-symmetry K-points. Without electric fields, the energy spectrum at K exhibits a gap of size |2a₂|. Interestingly, in the presence of SOI, the bulk gap can be closed by an increasing electric field perpendicular to the plane. The gap closes approximately at the very K-points, only if the corresponding parameter λ_r,1 is small compared to the transfer energy t, so λ_r,1 ≪ t. The low-energy Hamiltonian allows to give an analytical expression for the gap-closing critical field E^c_z [14]. At this point we may expect a topological phase transition from a QSH to a trivial insulating phase. This is indicated by the fact, that a gap is reestablished at the K-points, if the electric field is further increased, exceeding E^c_z. However, a rigorous exploration of prospective topological quantum phase transitions in the silicene parameter space requires the calculation of a topological bulk invariant.

Therefore, we now turn to the general calculation of the topological $\mathbb Z_2$-invariant defining the QSH phase [5]. While the $\mathbb Z_2$-invariant is of course a gauge invariant quantity by definition, the original literature [5, 16] did not provide a constructive recipe for its numerical calculation. This is so because a calculation following the original definition requires a macroscopic gauge, though giving a gauge invariant result. By macroscopic gauge, we mean that the phase relation between Bloch functions at remote points in the BZ has to be fixed. However, when the band structure is calculated numerically, such phase relations are typically not accessible thus preventing the direct calculation of the invariant. This problem has only been resolved rather recently [17, 26, 27] by a more constructive recipe for the direct calculation of the $\mathbb Z_2$-invariant. Here, we follow the method introduced by Prodan [17] which makes use of the elegant and manifestly gauge invariant formulation of the adiabatic connection [3, 28–30] originally introduced in a seminal work by Kato back in 1950 [28]. Recently, the manifestly gauge invariant formulation of topological invariants has been generalized to other topological band structures [3]. The crucial step of this approach consists in going from the Bloch functions of the occupied bands to the projection operator P(k) onto the occupied states. As already pointed out by Kato in [28], the advantage of this construction is that the projection operator is obviously basis (gauge) independent.

The adiabatic connection is defined as

$$\begin{equation} \mathcal A_\mu(k) = -\left[(\partial_\mu P(k)),P(k)\right], \end{equation} \tag{ 16 }$$

where $\partial _\mu =\frac {\partial }{\partial k^\mu }$ is the derivative in momentum space. Note that this connection is manifestly independent of the basis choice within the occupied bands in contrast to the more familiar Berry connection. Therefore, the adiabatic connection can be calculated numerically in a straight forward way which is in general not possible for the Berry connection. In [17], the $\mathbb Z_2$ invariant of the QSH state is constructed in a similar way to [16] but using the adiabatic connection instead of the Berry connection. We refer the reader to [17] for this both accessible and fairly self contained explicit construction and review only the resulting expression for the $\mathbb Z_2$-invariant Ξ for a rectangular lattice with lattice constants a_x, a_y, for completeness:

$$\begin{equation} \Xi= \alpha \frac{\det(U_{(0,0),(0,b_y)})\det(U_{(b_x, 0),(b_x, b_y)}) {\mathrm{Pf}}(\theta_{(0,0)}) {\mathrm{Pf}}(\theta_{(b_x,0)})}{{\mathrm{Pf}}(\theta_{(0,b_y)}) {\mathrm{Pf}}(\theta_{(b_x,b_y)}) \sqrt{\det(U_{(0, b_y),(0,-b_y)}) } \sqrt{\det(U_{(b_x, b_y),(b_x,-b_y)}) } }. \end{equation} \tag{ 17 }$$

Here, $b_x=\frac {\pi }{a_x},~b_y=\frac {\pi }{a_y}$ are the boundaries of the BZ, θ_(kx,ky) is the matrix of the TRS operation which is anti-symmetric at the time-reversal invariant momenta, and Pf denotes the Pfaffian. We note that Ξ = −1 defines the non-trivial QSH phase whereas Ξ = 1 for a trivial insulator. The parameter α will be explained in detail below and

$$\begin{equation} U_{(k^{(\mathrm{i})},k^{(\mathrm{f})})}=\mathcal T \,\mathrm{e}^{-\int_{k^{(\mathrm{i})}}^{k^{(\mathrm{f})}} \mathcal A} \end{equation} \tag{ 18 }$$

describes the adiabatic evolution along the straight line from k⁽ⁱ⁾ to k^(f) with the path ordering operator $\mathcal T$. Unitaries generated by the adiabatic connection such as equation (18) can be conveniently calculated numerically as products of projection operators [3, 17, 30]. Explicitly,

$$\begin{equation} U_{(k^{(\mathrm{i})},k^{(\mathrm{f})})}=\lim_{n\rightarrow \infty}\prod_{j=0}^n P(k_j),\quad k_j= k^{(\mathrm{i})} + j~\frac{k^{(\mathrm{f})}-k^{(\mathrm{i})}}{n}. \end{equation} \tag{ 19 }$$

The path ordering appearing in equation (18) then amounts to the ordering of the product in equation (19) from the right to the left with increasing j. The construction resulting in equation (19) is similar to a Trotter decomposition for a time evolution operator and is correct to leading order in $\frac {1}{n}$. The gauge invariance of equation (17) might not be obvious at first glance due to the basis dependence of the Pfaffians. However, the combinations

$$\begin{equation} \xi_{k_x}=\frac{\det(U_{(k_x,0),(k_x,b_y)}){\mathrm{Pf}}(\theta_{(k_x,0)})}{{\mathrm{Pf}}(\theta_{(k_x,b_y)})}=\pm \sqrt{\det(U_{(k_x, b_y),(k_x,-b_y)})} \end{equation} \tag{ 20 }$$

with k_x = 0, b_x appearing in equation (17) are indeed gauge invariant up to the choice of the branch of the multivalued square-root as is obvious from the right-hand side of equation (20). That is the point where α comes into play. When calculating

$$\begin{eqnarray*} &&\Xi=\alpha\frac{\xi_{0}\xi_{b_x}}{\sqrt{\det(U_{(0, b_y),(0,-b _y)})}\sqrt{\det(U_{(b_x, b_y),(b_x,-b_y)}}} \end{eqnarray*}$$

it has to be assured that the branch choices of the square-root at k_x = 0 and b_x are continuously connected to each other. This can be done by continuously interpolating the phase factor $\det (U_{(k_x, b_y),(k_x,-b_y)})$ between k_x = 0 and b_x [17]. If this phase has a winding such that the standard branch cut between Riemann sheets (line (− ∞,0] in the complex plane) is crossed n times during this interpolation, the naive result for Ξ is corrected by the factor α = (− 1)ⁿ.

While this procedure might be numerically challenging close to critical points or for large super-cells in disordered systems, it is basically a straightforward recipe. Note that the arbitrary basis choice for the representation matrix θ_(kx,ky) does not require any knowledge about the relative phase between the basis vectors at different points in k-space. Furthermore, the mentioned phase interpolation does not require any numerically inaccessible information either, since the phase factor $\det (U_{(k_x, b_y),(k_x,-b_y)})$ at every k_x is a gauge invariant quantity.

On the one hand, the procedure for the calculation of the $\mathbb Z_2$-invariant just described is general but requires a rectangular lattice whereas silicene crystalizes in a honeycomb lattice. On the other hand, the result for the topological invariant cannot depend on the choice of the unit cell. Therefore, we will now introduce a rectangular super cell which immediately allows us to apply the above method. To this end, the common basis vectors given in equation (1) are combined to a new set of basis vectors

$$\begin{equation} \vec{a}'_1=\vec{a}_1+\vec{a}_2= \left(0,3 a,0\right), \quad \vec{a}'_2=\vec{a}_1-\vec{a}_2= \left(\sqrt{3} a,0,0\right) \end{equation} \tag{ 21 }$$

spanning a rectangular lattice with four atoms per site A, B, A' and B', as shown in figure 3. The BZ is again rectangular with basis vectors

$$\begin{equation} \vec{b}'_1= \left(0,\frac{2 \pi}{3 a},0\right), \quad \vec{b}' _2= \left(\frac{2 \pi}{\sqrt{3} a},0,0\right). \end{equation} \tag{ 22 }$$

Zoom In Zoom Out Reset image size

The Bloch Hamiltonian of size (8 × 8) associated with our tight-binding model consists of the following contributions:

$$\begin{equation} H=H_{\mathrm{t}}+H_{\rm so}+H_{\mathrm{e}}+H_{\rm e,2}+H_{\lambda_{r,1}}+H_{\lambda_{r,2}}. \end{equation} \tag{ 23 }$$

The form of the k-dependent matrices is listed in equations (B.1)–(B.6) in appendix B. In the extended unit cell, the K-points (former corner points of the BZ) are mapped onto points $(\pm \frac {2 \pi }{3 \sqrt {3} a},0)$ inside the rectangle forming the new BZ.

For the Hamiltonian we chose the basis {|ψ_A〉,|ψ_B〉,|ψ_A'〉,|ψ_B'〉} × {↑,↓}, where the wavefunctions are of the form [31]

$$\begin{eqnarray*} &&|\psi_X \rangle =\frac{1}{\sqrt{N}}\sum_j \,\mathrm{e}^{\imath \vec{k}\cdot \vec{R}_j} |\phi_j^X \rangle \end{eqnarray*}$$

with X∈{A,B,A',B'}. Here, ϕ are atomic wavefunctions, in particular $\langle \vec {r}|\phi _j^X\rangle = \phi (\vec {r}-\vec {R}_j^X)$. $\vec {R}_j$ denote lattice vectors in real space, connecting the reference points of unit cells, while $\vec {R}_j^X$ are the positions of atoms labeled X, in the unit cell j, relative to the reference point. The sum runs over all cells of the crystal.

Having constructed a rectangular lattice for silicene with four atoms per unit cell (see figure 3), we now investigate its topological phase diagram by direct evaluation of equation (17). The unitary adiabatic time evolutions (see equation (18)) appearing in equation (17) are evaluated numerically using equation (19). The numerical calculations are done with a k-space discretization mesh of n = 200 (see equation (19)) and 200 steps for the interpolation determining the phase winding α (see discussion below equation (20)). Close to the critical points, we have increased both discretizations from 200 up to 1000 steps to reach convergence. To ensure a regular evolution of the phase factor, we set the numerical threshold for the absolute value of $\det (U)$ to be greater than 0.7.

For notational simplicity, the electric field E_z is absorbed in the parameters λ_i, that is λ_iE_z → λ_i. All parameters are measured in units of the hopping matrix element t, which is known to be of the order of 1.6 eV in silicene [13]. We keep in mind that λ_e, λ_e,2 and λ_r,1 depend on the external electric field. We first benchmark the general method in a parameter regime where we have a good intuition for the expected QSH physics from previous literature [4, 14]. Let us start with just one non-vanishing parameter t, all other parameters being zero. This choice corresponds to the well known gapless Dirac cones at the K-points. Upon switching on λ_so [4], the system exhibits a bulk gap of size |2λ_so| and we reproduce Ξ = −1 [5]. We now add a term proportional to λ_e depending on the external field. Upon increasing λ_e the bulk gap closes at the K-points for a critical value as predicted before. For even stronger fields a trivial gap opens, i.e. Ξ = + 1. Similarly, the gap closes for terms proportional to λ_e,2. Next, we add the Rashba-term proportional to λ_r,1 and the intrinsic SOI-term proportional to λ_r,2. The first term λ_r,1 tends to establish three additional minima of the energy gap, that are shifted away from the K-point. On the other hand, λ_r,2, primarily providing gapless states at the corners of the BZ, causes a modification of the band slope. From an analysis of the overall band structure, one finds, that the interplay of both terms λ_r,1 and λ_r,2 may possibly close the band gap away from the K-points (see figure 4).

Zoom In Zoom Out Reset image size

After these general considerations, we would now like to focus on the interesting parameter regime identified in [15], i.e. we consider λ_e, λ_r,1, λ_r,2 as free parameters that are measured in units of t. All other parameters are set to zero. Unfortunately, the analysis in [15] was not suitable to predict the correct phase diagram for the $\mathbb Z_2$ invariant. However, our more rigorous analysis confirms the general conjecture that in this regime, a combination of λ_r,1 and λ_r,2 is able to drive a topological quantum phase transition away from the K-points which induces a non-trivial QSH phase. If we only switch on λ_e, a trivial gap of size |2λ_e| opens. By tuning λ_r,1 and λ_r,2, the bulk gap can be closed away from the K-points, and a gap characterized by Ξ = −1, i.e. a QSH state emerges. We present the phase diagram in the identical parameter regime as presented in [15] in figure 5. Our direct calculation of the $\mathbb Z_2$-invariant disagrees with [15] both qualitatively (absence of the QSHE2) and quantitatively (significantly different phase boundary for the QSH phase). However, we would again like to point out that we can definitely confirm the phenomenology of a QSH phase in the absence of the Kane–Mele term λ_so. Instead, the extended QSH region shown in figure 5 is only driven by the two silicene specific (Rashba) SOI terms λ_r,1 and λ_r,2 which is conceptually very interesting. The existence of gapless edge states was additionally verified by the simulation of zigzag-edge silicene nanoribbons, as shown in figure 6. However, we emphasize, that our topological analysis is based on bulk properties and does not depend on specific forms of the boundaries.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Lastly, we found the phase diagram in figure 5 to be robust against small perturbations proportional to the term λ_so. For larger Kane–Mele parameters $\lambda _{\mathrm { so}} > \lambda _{\mathrm { e}}/(3 \sqrt {3})$, another non-trivial regime (Ξ = −1) enters the phase diagram. This regime occurs independently of λ_r,2 and can be expelled again by increasing λ_r,1. Essentially, the previously described quantum phase transitions are not affected.

Let us, for completeness, shortly review the general theory of the invariant expansion.

In this approach, the Hamiltonian is constructed to be invariant under all symmetry operations of the point group of the underlying lattice. When spin is included, it is composed block-wise according to IR of the corresponding double group. Each block $\mathcal {H}^{\alpha \beta }(\vec {\mathcal {K}})$, that may depend on some general tensor components $\vec {\mathcal {K}}$, has to reflect the symmetry properties of all the IRs Γ_κ, that are contained in the product

$$\begin{eqnarray*} &&\Gamma_{\alpha}\times \Gamma_{\beta}^*=\sum_{\kappa} n_{\kappa} \Gamma_{\kappa}, \end{eqnarray*}$$

where the integer n_κ is the multiplicity. The most general ansatz is then

$$\begin{equation} \mathcal{H}^{\alpha \beta}(\vec{\mathcal{K}})=\sum_{\kappa, \lambda, \mu} a_{\lambda\mu}^{\alpha\beta, \kappa} \sum_{l} X_{l}^{(\kappa, \lambda)} \mathcal{K}_l^{(\kappa, \mu)\ *}.\end{equation} \tag{ A.1 }$$

Here, κ runs over all IR Γ_κ included in the product and a^αβ,κ_λμ are material-specific constants. The symmetrized tensor operator components $\mathcal {K}_l^{(\kappa , \mu )\ *}$ are parameters like the wave vector $\vec {k}$, as well as the external magnetic $\vec {\mathcal {B}}$ or electric field $\vec {\mathcal {E}}$. Furthermore, X^(κ,λ)_l are the symmetrized matrices we wish to find to construct the Hamiltonian. Evidently, λ and μ are indices numbering the different possible matrix- and tensor-components. These indices run from unity to limits that depend on the order of the expansion. One set of the κ, λ and μ is called an invariant.

The irreducible tensor components belonging to an IR are composed analogously to its eigenfunctions. With the help of projection operators [32], that contain the matrix representations of each IR, we can combine tensor components at any order and project them onto the appropriate IR of the point group. To construct the basis matrices X^(κ,λ)_l, we can use the Wigner–Eckart theorem [33]. Since any point group is a subgroup of the full rotation group $\mathcal {R}$, the eigenfunctions of any IR of the point group are also eigenfunctions of $\mathcal {R}$ and can be written in terms of spherical harmonics. Angular momentum quantum numbers are assigned to each IR, by defining the axial vector components R_x, R_y and R_z to be the real angular momentum eigenstates with quantum number l = 1. Each IR can now be classified with angular momentum quantum numbers by comparing its eigenfunctions to a table of real spherical harmonics. A symmetrized matrix X^k*_ν (that transforms like the IR Γ*_k of dimension l_k and numbering ν = 1,...,l_k) which is contained in the block Γ*_i × Γ_j is then derived by Clebsch–Gordan coefficients (CGC)

$$\begin{equation} (X^{k *}_{\nu})_{\lambda \mu}=\left( \begin{array}{@{}cc|c@{}} k & j & i\ l \\ \nu & \mu & \lambda \\ \end{array} \right)^*, \end{equation} \tag{ A.2 }$$

where the order of i, j and k is crucial. In equation (A.2), l is again the multiplicity, that shows how often the IR Γ_k is included in the product Γ*_i × Γ_j. For multiplicities larger than one, we find linear independent sets of basis matrices.

We now start to perform an invariant expansion of the group D₃, that we found to be the group of the wave vector at the K-points in silicene (see [34] for a general discussion of this point group). We are interested in the two-dimensional π-bands of silicene. The orbital part of the wavefunction corresponds to the two-dimensional IR Γ₃ of the group D₃. Since we wish to include spin in our symmetry analysis, we have to complement this IR by Γ₄, the appropriate double group IR of D₃. Consequently, we start from the following product of representations Γ₃ × Γ*₄ = Γ₄ + Γ₅. The resulting (4 × 4) Hamiltonian based on the invariant expansion can then be composed in blocks of four (2 × 2) matrices

$$\begin{equation} \mathcal{H}= \left(\begin{array}{@{}cc@{}} \mathcal{H}_{44} & \mathcal{H}_{45} \\ \mathcal{H}_{54} & \mathcal{H}_{55} \end{array}\right).\end{equation} \tag{ A.3 }$$

The construction by CGCs (as discussed in the previous section) implies that this π-band Hamiltonian is given in the basis (ψ_Aβ,ψ_Bα,−ψ_Aα,ψ_Bβ)^T, with ψ_A(B) ≡ R_x ∓ ıR_y and α = |↑〉, β = |↓〉. Each block of the Hamiltonian is thus decomposed into IRs. General tensor components of any desired order can be found with the help of projection operators from the point group character table. In table A.1, we list components of interest up to third order. Identifying combinations of general tensor components of $\vec {R}$ with angular momentum quantum numbers, we can assign such quantum numbers to the IRs itself. The basis matrices are then identified from CGCs according to equation (A.2) and listed in table A.2. Note that we have applied a unitary basis transformation to guarantee a more symmetric basis (ψ_Aβ,ψ_Bα,ψ_Aα,ψ_Bβ)^T. This transformation is done for comparison with the graphene case carefully analyzed in [21]. We use the symmetrized matrices and tensor components listed in tables A.1 and A.2 to expand the Hamiltonian up to first orders in $\vec {k}$ and the electric field E_z perpendicular to the plane. Note that since $\mathcal {H}$ is Hermitian, the coefficients of the diagonal blocks have to be real, while those of the off-diagonal blocks can in principle be imaginary. We write γ and ıγ' = γ to include both cases with real coefficients γ, γ'. Then, we obtain for the Hamiltonian

$$\begin{eqnarray}&&\fl \mathcal{H}_{44}= \alpha_1 \mathbb{I} + \alpha_2 \sigma_z E_z + \alpha_3 (\sigma_x k_x - \sigma_y k_y)+\alpha_4 E_z (\sigma_x k_y +\sigma_y k_x), \end{eqnarray} \tag{ A.4 }$$

$$\begin{eqnarray}&&\fl \mathcal{H}_{55}= \beta_1 \mathbb{I} + \beta_2 \sigma_x 1 +\beta_3 \sigma_z E_z + \beta_4 \sigma_y E_z, \end{eqnarray} \tag{ A.5 }$$

$$\begin{eqnarray} &&\fl \mathcal{H}_{45}= \gamma_1 (\mathbb{I} k_x -\imath \sigma_z k_y)+\imath \gamma_1' (\mathbb{I} k_x -\imath \sigma_z k_y)+ \gamma_2 (\sigma_x k_x + \sigma_y k_y)+ \imath \gamma_2' (\sigma_x k_x + \sigma_y k_y) \nonumber \\ &&+ \gamma_3 E_z (\mathbb{I} k_y +\imath \sigma_z k_x)+ \imath \gamma_3' E_z (\mathbb{I} k_y +\imath \sigma_z k_x)+ \gamma_4 E_z (\sigma_x k_y -\sigma_y k_x) \nonumber \\ &&+\imath \gamma_4' E_z (\sigma_x k_y -\sigma_y k_x). \end{eqnarray} \tag{ A.6 }$$

For a better physical interpretation, we once more change the basis by an unitary transformation and present the total Hamiltonian in the basis (ψ_Aβ,ψ_Aα,ψ_Bβ,ψ_Bα)^T. After this basis transformation, the Pauli matrices σ and s obtain the following physical meaning: σ acts on the space of sublattices A and B, and s on the spin space. In lowest order, the Hamiltonian in our new basis exhibits 16 terms labeled by coefficients a₁ to a₁₆:

$$\begin{eqnarray}&&\fl \mathcal{H}^K = a_1\mathbb{I}+ a_2 \sigma_z s_z + a_3 \sigma_z s_0 E_z +a_4 \sigma_0 s_z E_z + a_5 ( \sigma_x s_y - \sigma_y s_x) E_z+ a_{6} ( \sigma_x s_x + \sigma_y s_y)\nonumber\\ &&+a_7 ( \sigma_x k_x +\sigma_y k_y) s_0 + a_8 \sigma_0 (s_x k_x +s_y k_y) + a_9[\sigma_x(s_x k_x - s_y k_y)-\sigma_y(s_y k_x +s_x k_y)] \nonumber \\ &&+ a_{10} E_z [\sigma_x(s_x k_y + s_y k_x)+\sigma_y(s_x k_x -s_y k_y)]+ a_{11} E_z \sigma_0 (s_x k_y - s_y k_x)\nonumber \\ &&+ a_{12} E_z (\sigma_x k_y - \sigma_y k_x) s_0 + a_{13}\sigma_z (s_x k_y - s_y k_x) + a_{14} (\sigma_x k_y - \sigma_y k_x) s_z\nonumber \\ &&+ a_{15} E_z \sigma_z (s_x k_x+s_y k_y)+a_{16}E_z(\sigma_x k_x +\sigma_y k_y)s_z. \end{eqnarray} \tag{ A.7 }$$

Note, that the coefficients in equations (A.4)–(A.6) and in equation (A.7) are consistent, in the sense, that they are related by mutual linear combinations. Here, we would already like to point out the additional, interesting spin–orbit term proportional to the coefficient a₄, coupling directly the out-of-plane components of spin and electric field.

Table A.2. Symmetrized matrices for blocks arising from the product Γ₃ × Γ*₄ = Γ₄ + Γ₅ of the double group of D₃. Our choice of basis is (ψ_Aβ,ψ_Bα,ψ_Aα,ψ_Bβ)^T. σ_i denote the (2 × 2) Pauli matrices with i∈{x,y,z}.

$\mathcal {H}_{44}$	Γ4×Γ*4=Γ1+Γ2+Γ3	$\Gamma _1: \mathbb {I}$
		Γ2:σz
		Γ3:(σx,−σy)
$\mathcal {H}_{55}$	Γ5×Γ*5=2Γ1+2Γ2	$\Gamma _1: \mathbb {I}; \sigma _x$
		Γ2:σz;σy
$\mathcal {H}_{45}$	Γ4×Γ*5=2Γ3	$\Gamma _3: (\mathbb {I},-\imath \sigma _z);(\sigma _x, \sigma _y)$

The terms derived so far in the invariant expansion have not yet been checked for consistency with TRS which holds throughout this work since we only consider the influence of electric fields. Importantly, the point group D₃ + σ'' provides symmetry operations that map the inequivalent corner points K and K' of the Brillouin zone onto each other. For example, we can consider the reflection at a plane perpendicular to the atomic plane and including the y-axis (called R_y in [21]). This symmetry operation has the following impact on the Hamiltonian:

$$\begin{eqnarray*} &&\mathcal{D}(R_y) \mathcal{H}^K(\vec{\mathcal{K}}) \mathcal{D}(R_y)^{-1}=\mathcal{H}^{K'}(R_y^{-1}\vec{\mathcal{K}}), \end{eqnarray*}$$

where $\mathcal {D}(R_y)$ is the matrix representation of the operation R_y. Polar ($\vec {x}$) and axial ($\vec {R}$) tensor components will then transform like R⁻¹_y:(x,y,z) = (− x,y,z) and R⁻¹_y:(R_x,R_y,R_z) = (R_x,−R_y,−R_z). The sublattices will not be changed by R_y. In a compact notation, we can write R_y = σ₀s_x. With the help of R_y we find the form of the Hamiltonian at one of the two inequivalent K points. All terms, that change sign under R_y will thus be modified by the additional parameter τ_z = ± 1 to mark the difference between the valleys K and K'. The true time reversal operator equals $\mathcal {T}= - \sigma _0 \tau _x s_y \mathcal {C}$, where $\mathcal {C}$ is the complex conjugation operator [35] and τ_x a Pauli matrix corresponding to the valley isospin. Its impact on the Hamiltonian can be written as

$$\begin{eqnarray*} &&\mathcal{T} \mathcal{H}^K(\vec{\mathcal{K}}) \mathcal{T}^{-1}=\mathcal{H}^{K' *}(\xi\vec{\mathcal{K}}), \end{eqnarray*}$$

where ξ = ± 1 depending on whether the tensor component $\vec {\mathcal {K}}$ changes sign under time reversal or not. Hence, both operators R_y and $\mathcal {T}$ lead to transformations between the valleys. We now combine both of them to a (new) time-reversal operator within a single valley Θ (in the spirit of [21]): $\Theta (R_y)=\mathcal {T} \mathcal {D}(R_y)=\imath s_z \mathcal {C}$. This operator yields an additional symmetry constraint for all terms

$$\begin{equation} \Theta \mathcal{H}(\vec{\mathcal{K}}) \Theta^{-1}=s_z \mathcal{H}^{*}(\xi R_y^{-1}\vec{\mathcal{K}})s_z \end{equation} \tag{ A.8 }$$

which forces the coefficients a₆, a₈, a₉ and a₁₂ in equation (A.7) to vanish. Our result is the low-energy Hamiltonian of silicene near the K-points (in first orders in $\vec {k}$ and E_z) presented in the basis (ψ_Aβ,ψ_Aα,ψ_Bβ,ψ_Bα)^T. It is given by

$$\begin{eqnarray}&&\fl \mathcal{H}^{K(K')} =a_1\mathbb{I}+ a_2 \tau_z \sigma_z s_z + a_3 \sigma_z s_0 E_z + a_4 \tau_z \sigma_0 s_z E_z+ a_5 ( \tau_z \sigma_x s_y - \sigma_y s_x) E_z+ a_7 (\tau_z \sigma_x k_x +\sigma_y k_y) s_0 \nonumber \\ &&+ a_{10} E_z [\sigma_x(s_x k_y + s_y k_x)+\tau_z\sigma_y(s_x k_x -s_y k_y)]+ a_{11} E_z \sigma_0 (s_x k_y - s_y k_x) \nonumber \\ &&+ a_{13}\sigma_z (s_x k_y - s_y k_x)+ a_{16}E_z(\sigma_x k_x +\tau_z\sigma_y k_y)s_z. \end{eqnarray} \tag{ A.9 }$$