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

Silicene consists of a monolayer of silicon atoms in a buckled honeycomb structure. It was recently discovered that the symmetry of such a system allows for interesting Rashba spin-orbit effects. A perpendicular electric field is able to couple to the sublattice pseudospin, making it possible to electrically tune and close the band gap. Therefore, external electric fields may generate a topological phase transition from a topological insulator to a normal insulator (or semimetal) and vice versa. The contribution of the present article to the study of silicene is twofold: First, we perform a group theoretical analysis to systematically construct the Hamiltonian in the vicinity of the $K$ points of the Brillouin zone and discover a new, but symmetry allowed term. Subsequently, we identify a tight binding model that corresponds to the group theoretically derived Hamiltonian near the $K$ points. Second, we start from this tight binding model to analyze the topological phase diagram of silicene by an explicit calculation of the $Z_2$ topological invariant of the band structure. To this end, we calculate the $Z_2$ topological invariant of the honeycomb lattice in a manifestly gauge invariant way which allows us to include $S_z$ symmetry breaking terms -- like Rashba spin orbit interaction -- into the topological analysis. Interestingly, we find that the interplay of two Rashba terms can generate a non-trivial quantum spin Hall phase in silicene. This is in sharp contrast to the more extensively studied honeycomb system graphene where Rashba spin orbit interaction is known to compete with the quantum spin Hall effect in a potentially detrimental way.


Introduction
One of the main subjects of current interest in condensed matter physics is the search for materials that host topological insulator (TI) phases [1][2][3]. Two dimensional TIs exhibit the quantum spin Hall effect (QSHE) with gapless edge states and a finite energy gap in the bulk [4][5][6]. The first proposal of this state of matter was made by Kane and Mele [4] on the basis of graphene in the presence of spin-orbit interaction (SOI). However, the relevant SOI in graphene turns out to be rather small [7] such that the effect seems to be inaccessible in experiments. This situation is different in HgTe/CdTe quantum wells where the QSHE was also predicted theoretically [6] and experimentally seen soon after [8].
Recently, a single layer of silicon atoms-called silicene-has been synthesized exhibiting an analogous honeycomb structure as graphene [9][10][11]. Since silicon is heavier than carbon, the spin-orbit gap in silicene is much larger than in graphene. Therefore, if it was possible at some point to prepare clean silicene, it should be feasible to experimentally access the QSHE in this material. Similar to graphene, the unit cell of silicene contains two atoms which gives rise to two different sublattices. In contrast to graphene, however, the silicene sublattices are found to be arranged in a buckled structure pointing out-of-plane [12]. Due to the broken sublattice symmetry, the mobile electrons in silicene are therefore able to couple differently to an external electric field than the ones in graphene. This difference is the origin of new (Rashba) 4 spin-orbit coupling effects that allow for external tuning and closing of the band gap in silicene [13]. Consequently, an electrically induced topological quantum phase transition is possible. It is natural to ask whether this phase transition can in principle go both ways, i.e. whether the electric field can be used to destroy and generate the QSHE. [4,14] clearly show that a different potential on the two sublattices of a honeycomb lattice leads to a transition from a TI to a trivial insulating state. In [15], some indications have been presented that the interplay of two silicene specific (Rashba) spin-orbit terms can even induce the QSHE starting from a trivial insulating band structure in the absence of these terms.
The quantum spin Hall (QSH) phase is distinguished from a normal insulating phase by a bulk Z 2 topological invariant [5,16]. For a minimal model of the QSHE in graphene, this invariant has been analytically calculated in a seminal work by Kane and Mele [5]. However, the original formulation of the Z 2 invariant in terms of Bloch functions does not contain a constructive prescription as to its numerical evaluation. Subsequent work on the topological properties of the band structure of silicene was restricted to the absence of terms breaking the spin S z -conservation or to employing the bulk-boundary correspondence in silicene nanoribbons [14,15].
Evidently, a full topological analysis of silicene is missing and, as we show below, important to clearly identify phenomenological differences between graphene and silicene. In this work, we employ Prodan's method [17] to calculate the topological invariant without any further symmetry assumptions in a manifestly gauge invariant way to provide a conclusive analysis of the novel features of silicene regarding QSH physics. In particular, we establish that, in contrast to graphene, the QSHE can be generated by (Rashba) SOI in silicene.
The bulk of this paper consists of two parts which we keep fairly self contained to allow the reader to follow our analysisà la carte. In section 2, we analyze the symmetries of the lattice of silicene. This analysis allows us to mathematically construct the low-energy Hamiltonian (close to the K points of the Brillouin zone (BZ)) by means of the invariant expansion method with a particular focus on terms involving a perpendicular electric field. Thereby, we discover for silicene an additional, electric field induced spin-orbit term of the low-energy Hamiltonian. Furthermore, a tight-binding calculation is performed to verify the terms previously derived from symmetries and to estimate their magnitude. The reader who is more interested in QSH physics can directly go to section 3, where we study the topological properties of the band structure of silicene by explicitly calculating the Z 2 topological invariant in a manifestly gauge invariant way. In this section, the possibility of a topological phase transition induced by an external field is carefully examined which enables us to correct previously proposed phase diagrams. Finally, we conclude in section 4. Some technical details of the invariant expansion and the tight-binding model are presented in the appendix.

Identification of the lattice symmetry
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.
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 The buckling-angle θ is found to be 101.7 • in a silicene lattice model [12]. . The operations C n denote rotations by 2π/n around the z-axis perpendicular to the plane; C n and C n refer to rotations around the labeled axis, lying within the atomic plane. σ describes reflection planes spanned by the labeled axis and the z-axis. The primed operations cross the corners of the underlying hexagon, while the double-primed ones do not. A and B denote different sublattices, K and K inequivalent corner points of the first BZ. The red dots implicate, that for silicene in real space the atomic sites are shifted perpendicular to the plane. in real space and by 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 6 , C 2 , C 2 , σ and σ h are broken and we are left with point group D 3d .
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 2 -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 3 , C 2 , σ , i and S 6 in reciprocal space, which corresponds to the point group D 3 with the additional symmetry classes {σ , i, S 6 }. At the K points of the BZ, the symmetry of the group of the wave vector is further reduced to the point group D 3 .

Invariant expansion
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][19][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 3 . The π -orbital wavefunction transforms like the two-dimensional irreducible representations (IR) 3 of the group D 3 , while the spin part is represented by the IR 4 of the double group. Therefore, our total wavefunction is of the form of the product 3 × * 4 = 4 + 5 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 k and E z ) 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 timereversal symmetry (TRS) have already been taken into account, as explained in appendix A. In equation (3), the terms proportional to a 2 , a 5 and a 7 are well known from the graphene literature to describe its fundamental properties near the K -points [4]. Further corrections proportional to a 10 and a 11 can as well be found in graphene [21]. Beyond this, silicene exhibits specific terms proportional to a 3 and a 13 as reported in [12,13]. Additionally, we find an electric field induced SOI-term proportional to a 4 , that has not been reported for silicene before 5 .
Moreover, the Hamiltonian exhibits higher order corrections to the linear dispersion (proportional to a 16 ). 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.

Corresponding tight-binding model
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] where F is the force stemming from the electric potential V , p is the momentum and 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]: 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 d i j with its corresponding unit vectord i j . The sign v i j = ±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 where i is the on-site energy of the atomic orbital and t i j 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 is generated [14] with a, so far, undefined parameter λ e . Interestingly, E i 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 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 4 ). The full tight-binding Hamiltonian is then given by The terms given above will reproduce our Hamiltonian derived by symmetry analysis, equation (3), when being expanded around the K -points.

Tight-binding model including π and σ -bands of silicene
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 Nearest-neighbor hopping matrix elements between given orbitals can then be expressed by parameters V 1 , V 2 and V 3 as functions of Slater bond parameters V ppπ , V ppσ , V spσ and curvature angle θ with (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: 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 0 the Stark element weighting the transition. Spin-orbit coupling is now included in the same way as in [12] by the term H SO = L · s where L is the angular momentum vector, s the spin operator and the coupling constant. The tightbinding model of π and σ bands allows us to estimate some of the parameters introduced in the previous section as a function of the parameters V 1−3 (that are in principle known) as well as the Stark element z 0 . To do so, we now apply a unitary transformation U to the Hamiltonian that corresponds to the following change of basis: with mixed orbitals, for example, ). Then, the transformed Hamiltonian can be analyzed perturbatively on the first (4 × 4) block corresponding to the basis {φ 1 ↑, φ 1 ↓, φ 4 ↑, φ 4 ↑}, 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: In the last step, the expressions were approximated under the assumption of low buckling, meaning θ → 90 • and V 3 → 0, where only the term of lowest order in V 3 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 3 , 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.

Numerical estimates
We now estimate the size of the spin-orbit terms coupling to an external electric field in silicene. The bond parameters V 1−3 , 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 0 = φ 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 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 .

Topological analysis
Triggered by the theoretical prediction [4][5][6] and experimental discovery [8] of the QSH effect, the study of topological effects in the physics of Bloch bands has been a major focus of condensed matter physics in recent years [1][2][3]. The QSH state is a bulk insulating state featuring metallic edge states that are protected by TRS. Due to Kramers theorem, these edge states appear in so called helical pairs. The bulk energy gap in the original proposal for the QSH effect in graphene by Kane and Mele [4] is due to an intrinsic SOI which preserves a residual U (1) spin symmetry. In the presence of this global spin quantization axis, the helical edge states are characterized by a perfect locking of spin and momentum. States with opposite spin move with opposite chirality along the edge. In the presence of (Rashba) SOI, no spin symmetry is present resulting in the absence of a global spin quantization axis. Rather, the spin quantization axis of the edge states can precess spatially. However, as TRS is present the Kramers pair of helical edge states is still protected from hybridizing, i.e. from a gap opening on the edge. This statement is true as long as the bulk gap is maintained implying that the system is adiabatically connected to the U (1) preserving case and hence is still in the QSH state. However, the role of Rashba SOI for the QSH effect in graphene is only detrimental. Upon increasing the strength of Rashba SOI, the system will go through a quantum phase transition destroying the QSH phase. Conversely, given the lattice symmetry of graphene, switching on Rashba spin-orbit coupling cannot generate the QSH phase starting from a semi-metallic or trivial insulator phase.
The slightly reduced symmetry of silicene stemming from its buckled structure entails several new couplings, at least two of which are of key relevance regarding QSH physics. Firstly, let us consider the staggered potential distinguishing sublattice A and B that has been introduced formally in [4] to open a trivial gap. While this term is symmetry forbidden in graphene, it has been shown [14] to be induced by a simple out of plane electric field in silicene as we confirmed in our symmetry analysis above. This provides a knob to experimentally switch off the QSHE. Secondly, and even more interestingly, it has very recently been conjectured [15] that the QSHE can be generated by virtue of SOI terms that are symmetry forbidden in graphene but allowed in silicene. The authors of [15] probe the sub-gap conductance as a fingerprint of the QSH state. While a non-vanishing conductance in a phase with a bulk gap is a promising signature of the QSH state, it is not in one to one correspondence with the Z 2 -invariant defining the QSHE [5].
For example, what is called the QSHE2-phase in [15] is a trivial insulator which only shows the characteristic conductance of 2 e 2 h expected for the QSH phase since one additional pair of edge states does not contribute due to a mini-gap at the ribbon sizes considered in [15]. In the thermodynamic limit, the conductance in this parameter regime would be 4 e 2 h signaling a Z 2 -trivial phase.
The main purpose of this section is to calculate the correct phase diagram of silicene in the presence of the (Rashba) SOI terms by rigorous calculation of the Z 2 invariant in the absence of any symmetries besides TRS. We note that our method can be applied to study the entire parameter space of the silicene band structure without any further complications. However, in this work, we would like to focus on a particularly interesting parameter regime where the QSH phase is generated by (Rashba) SOI. Previous work on the topology of the QSH state in silicene in the presence of (Rashba) SOI [14] has been focused on effective models valid close to the K -points. However, in the presence of (Rashba) SOI, the bulk gap can close away from the K -points, a phase transition which would be missed by such power series expansions. Our analysis does not suffer from these limitations since we directly calculate the Z 2 -invariant characterizing the QSH phase from its very definition [5,17]. We find that there is indeed a QSH phase in the absence of the original Kane-Mele term λ so . This QSH effect can be switched on by only tuning silicene specific s z symmetry breaking SOI terms and the staggered potential-all terms that are known to have only detrimental effect as to QSH physics in graphene. Our analysis, hence, settles in the affirmative the discussion on whether other SOI terms than the s z conserving Kane-Mele term can in principle generate a QSH state.

Manifestly gauge invariant calculation of the topological Z 2 -invariant
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 2 |. 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 Z 2 -invariant defining the QSH phase [5]. While the 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 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][29][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 where ∂ µ = ∂ ∂k µ 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 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 Z 2 -invariant for a rectangular lattice with lattice constants a x , a y , for completeness: Here, b x = π a x , b y = π a y are the boundaries of the BZ, θ (k x ,k y ) 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 describes the adiabatic evolution along the straight line from k (i) to k (f) with the path ordering operator 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, 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 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 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 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) n . 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 θ (k x ,k y ) 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.

From the honeycomb lattice of silicene to a rectangular super-cell
On the one hand, the procedure for the calculation of the 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 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 The Here, φ are atomic wavefunctions, in particular r |φ X j = φ( r − R X j ). R j denote lattice vectors in real space, connecting the reference points of unit cells, while R X j 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.

Results for the phase diagram of silicene
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 λ i E 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).
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 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 nontrivial 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,  The staggered potential λ e = 0.1 is fixed as well as λ so = 0. All couplings are measured in units of the hopping energy t. The white points close to the phase boundary denote critical regions where the bulk gap is so small that our numerical calculation did not converge properly.
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 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. 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 λ so > λ e /(3 √ 3), another nontrivial 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.

Conclusion
A symmetry analysis of the lattice of silicene close to the K points of the Brillouin zone was performed. With the help of the invariant expansion model, the π-band Hamiltonian was constructed by symmetry considerations only, including spin-orbit coupling and external electric fields perpendicular to the atomic plane. Thereby, we discovered an additional term that is allowed by symmetry and related to an interplay of spin-orbit coupling and an external electric field in a buckled honeycomb structure. We supplemented the symmetry analysis by a tight-binding model which allowed us to estimate the relevant parameters of the model. Subsequently, we carefully analyzed the topological properties of the band structure and proved that a topological phase transition can be generated by an external electric field. This analysis enabled us to plot a topological phase diagram of silicene employing Prodan's manifestly gauge-invariant method for the direct calculation of the Z 2 topological invariant defining the QSH phase. Since this method requires a rectangular lattice, we considered a rectangular silicene super-cell that contains four atoms. Interestingly, the tunable phase transition can happen in a destructive as well as a constructive way, i.e. a QSH phase cannot only be destroyed but also be generated by means of external (Rashba) spin-orbit coupling. x, y, z denote polar vector components and R x , R y , R z axial vector components. Here, axial and polar vector components appear on equal footing as the symmetry operations of D 3 do not provide mirror planes. Therefore, for simplicity, combinations with components of R are not listed again in higher orders but they are of course present.
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) 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.

A.2. Application to the π -band Hamiltonian of silicene near the K points
We now start to perform an invariant expansion of the group D 3 , 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 3 of the group D 3 . Since we wish to include spin in our symmetry analysis, we have to complement this IR by 4 [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 k and the electric field E z perpendicular to the plane. Note that since 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 H 45 = γ 1 (Ik x − ıσ z k y ) + ıγ 1 (Ik x − ıσ z k y ) + γ 2 (σ x k x + σ y k y ) + ıγ 2 (σ x k x + σ y k y ) 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 1 to a 16 : H K = a 1 I + a 2 σ z s z + a 3 σ z s 0 E z + a 4 σ 0 s z E z + a 5 (σ x s y − σ y s x )E z + a 6 (σ x s x + σ y s y ) +a 7 (σ x k x + σ y k y )s 0 + a 8 σ 0 (s x k x + s y k y ) + a 9 [σ x (s x k x − s y k y ) − σ y (s y k x + s x k y )] +a 10 E z [σ x (s x k y + s y k x ) + σ y (s x k x − s y k y )] + a 11 E z σ 0 (s x k y − s y k x ) +a 12 E z (σ x k y − σ y k x )s 0 + a 13 σ z (s x k y − s y k x ) + a 14 (σ x k y − σ y k x )s z +a 15 E z σ z (s x k x + s y k y ) + a 16 E z (σ x k x + σ y k y )s z .
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 4 , coupling directly the out-of-plane components of spin and electric field.

A.3. Consequences of time-reversal symmetry
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 3 + σ 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: where D(R y ) is the matrix representation of the operation R y . Polar ( x) and axial ( R) tensor components will then transform like R −1 y : (x, y, z) = (−x, y, z) and R −1 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 = σ 0 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 T = −σ 0 τ x s y C, where 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 where ξ = ±1 depending on whether the tensor component K changes sign under time reversal or not. Hence, both operators R y and 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]): (R y ) = T D(R y ) = ıs z C. This operator yields an additional symmetry constraint for all terms which forces the coefficients a 6 , a 8 , a 9 and a 12 in equation (A.7) to vanish. Our result is the low-energy Hamiltonian of silicene near the K -points (in first orders in k and E z ) presented in the basis (ψ A β, ψ A α, ψ B β, ψ B α) T . It is given by H K (K ) = a 1 I + a 2 τ z σ z s z + a 3 σ z s 0 E z + a 4 τ z σ 0 s z E z + a 5 (τ z σ x s y − σ y s x )E z + a 7 (τ z σ x k x + σ y k y )s 0 +a 10 E z [σ x (s x k y + s y k x ) + τ z σ y (s x k x − s y k y )] + a 11 E z σ 0 (s x k y − s y k x ) +a 13 σ z (s x k y − s y k x ) + a 16 E z (σ x k x + τ z σ y k y )s z . (A.9) The terms contributing to the tight-binding Hamiltonian in equation (23)