![FIG. 9. (Color online) DFT gap with different exchangecorrelation functionals (LDA and PBE) and pseudopotentials [on-thefly ultrasoft (USP) and norm conserving (NC)15] for silicene subject to an electric field of 0.257 V/Å in a cell of length Lz = 13.35 Å with a 15 × 15 k-point grid and a plane-wave cutoff energy of 816 eV.](/figures/fig-9-color-online-dft-gap-with-different-351ai20z.webp)
PHYSICAL REVIEW B 85, 075423 (2012)
Electrically tunable band gap in silicene
N. D. Drummond, V. Z
´
olyomi, and V. I. Fal’ko
Department of Physics, Lancaster University, Lancaster LA1 4YB, United Kingdom
(Received 20 December 2011; revised manuscript received 10 February 2012; published 22 February 2012)
We report calculations of the electronic structure of silicene and the stability of its weakly buckled honeycomb
lattice in an external electric field oriented perpendicular to the monolayer of Si atoms. The electric field produces
a tunable band gap in the Dirac-type electronic spectrum, the gap being suppressed by a factor of about eight by
the high polarizability of the system. At low electric fields, the interplay between this tunable band gap, which is
specific to electrons on a honeycomb lattice, and the Kane-Mele spin-orbit coupling induces a transition from a
topological to a band insulator, whereas at much higher electric fields silicene becomes a semimetal.
DOI: 10.1103/PhysRevB.85.075423 PACS number(s): 73.22.Pr, 61.48.Gh, 63.22.Rc
I. INTRODUCTION
Two-dimensional (2D) carbon crystals are hosts for Dirac-
type electrons, whose unusual properties have been studied
extensively in graphene monolayers produced by mechanical
exfoliation from graphite.
1,2
A close relative of graphene,
a 2D honeycomb lattice of Si atoms called silicene,
3
does
not occur in nature, but nanoribbons of silicene have been
synthesized on metal surfaces.
4–6
Due to the similarity of
the lattice structures, the band structure of silicene resembles
that of graphene, featuring Dirac-type electron dispersion in
the vicinity of the corners of its hexagonal Brillouin zone
(BZ).
7
Moreover, silicene has been shown theoretically to be
metastable as a freestanding 2D crystal,
3
implying that it is
possible to transfer silicene onto an insulating substrate and
gate it electrically. In this work we predict the properties of
this 2D crystal.
The similarity between graphene and silicene arises from
the fact that C and Si belong to the same group in the periodic
table of elements. However, Si has a larger ionic radius,
which promotes sp
3
hybridization, whereas sp
2
hybridization
is energetically more favorable in C. As a result, in a 2D
layer of Si atoms, the bonding is formed by mixed sp
2
and
sp
3
hybridization. Hence silicene is slightly buckled, with
one of the two sublattices of the honeycomb lattice being
displaced vertically with respect to the other, as shown in Fig. 1.
Such buckling creates new possibilities for manipulating the
dispersion of electrons in silicene and opening an electrically
controlled sublattice-asymmetry band gap.
8
In this paper we
report density functional theory (DFT) calculations of the
band gap for Dirac-type electrons in silicene opened by
a perpendicular electric field using a combination of top and
bottom gates. We show that can reach tens of meV before
the 2D crystal transforms into a semimetal and then, at still
higher fields, loses structural stability. We also determine the
weak electric field at which electrons in silicene experience
a transition from a topological insulator regime
9,10
caused by
the Kane-Mele spin-orbit (SO) coupling
11
for electrons on a
honeycomb lattice into a conventional band insulator regime.
The rest of this paper is arranged as follows. In Sec. II
we report our results for the structural and electronic proper-
ties of freestanding silicene, and compare them with other
theoretical and experimental results in the literature. In
Sec. III we analyze the effects of a transverse electric field
on the structural and electronic properties of silicene, and
in Sec. IV we discuss the effects of SO coupling on the
electronic structure, arguing that a crossover from topological
insulating behavior to band insulating behavior must take place
as the transverse field increases in strength. In Sec. V we give
the technical details of our computational methodology and
demonstrate the convergence of our results with respect to
simulation parameters. Finally, we draw our conclusions in
Sec. VI.
II. STRUCTURAL AND ELECTRONIC PARAMETERS
OF FREESTANDING SILICENE
A. Comparison with theoretical and experimental results
in the literature
The lattice constant and the z (out-of-plane) coordinates
of the Si atoms lying on the 2D honeycomb lattice were
both fully relaxed using DFT (i) in the local density approx-
imation (LDA), (ii) with the Perdew-Burke-Ernzerhof (PBE)
exchange-correlation functional,
12
and (iii) with the screened
Heyd-Scuseria-Ernzerhof 06 (HSE06) hybrid functional.
13,14
Our DFT calculations were performed using the CASTEP
15,16
and VASP
17
plane-wave-basis codes, using ultrasoft pseudopo-
tentials and the projector-augmented-wave (PAW) method,
respectively. The z coordinates of the two Si atoms in the
unit cell (the A and B sublattices) differ by a finite distance
z. Our results for a free silicene monolayer are shown in
Table I. The metastable lattice that we find is the same as the
“low-buckled” structure found by Cahangirov et al.
3
The experimental results for the lattice parameter depend on
the choice of substrate on which the silicene is grown.
5,6
The
extent to which theoretical results obtained for freestanding
silicene are applicable to the silicene samples that have been
produced to date is therefore unclear.
B. Stability of freestanding silicene
The cohesive energy of bulk Si (including a correction for
the zero-point energy) has been calculated within DFT-LDA
as 5.34 eV.
19
Comparing this with our DFT-LDA cohesive
energy of silicene reported in Table I shows that bulk Si is
substantially (0.22 eV per atom) more stable than silicene,
implying that silicene would not grow naturally as a layered
bulk crystal like graphite. However, by calculating the DFT
phonon dispersion it has been verified both here and in
Ref. 3 that the structure is dynamically stable: No imaginary
075423-1
1098-0121/2012/85(7)/075423(7) ©2012 American Physical Society
N. D. DRUMMOND, V. Z
´
OLYOMI, AND V. I. FAL’KO PHYSICAL REVIEW B 85, 075423 (2012)
FIG. 1. (Color online) Atomic structure of silicene, together with
a sketch of the charge density for the highest occupied valence band
in the vicinity of the K point.
frequencies appear anywhere in the BZ. The results of such
an analysis are summarized in Fig. 2. This convinces us that,
as a metastable 2D crystal, silicene can be transferred onto
an insulating substrate, where its electronic properties can be
studied and manipulated as suggested below.
C. Electronic band structure
The calculated band structure of a “free” silicene layer is
shown in Fig. 3. As expected, it resembles the band structure
of graphene; in particular it shows the linear Dirac-type
dispersion of electrons near the K points, where we find the
Fermi level in undoped silicene. The Fermi velocity v of
electrons in silicene is lower than that in graphene (see Table I).
Although the lattice parameters and sublattice buckling found
in the different DFT calculations are in good agreement, our
results for the Fermi velocity are very much smaller than the
Fermi velocity reported in Ref. 3.
III. APPLICATION OF A TRANSVERSE ELECTRIC FIELD
A. Breaking the sublattice symmetry
To exploit the weak buckling of silicene, we consider its be-
havior in an external electric field E
z
applied in the z direction,
TABLE I. Silicene structural and electronic parameters: lattice
constant a, sublattice buckling z (the difference between the
z coordinates of the A and B sublattices), cohesive energy E
c
,and
Fermi velocity v. The calculated cohesive energy of silicene includes
the DFT-PBE zero-point energy, which we found to be 0.10 eV
per atom. The theoretical results are for freestanding silicene; the
experimental results are for silicene nanoribbons on Ag substrates.
Method a (
˚
A) z (
˚
A) E
c
(eV) v (10
5
ms
−1
)
PBE (CASTEP)3.86 0.45 4.69 5.27
PBE (
VA S P )3.87 0.45 4.57 5.31
PBE
8
3.87 0.46
LDA (
CASTEP)3.82 0.44 5.12 5.34
LDA (
VA S P )3.83 0.44 5.00 5.38
LDA
3
3.83 0.44 5.06 ≈10
LDA
18
3.86 0.44
HSE06 (
VA S P )3.85 0.36 4.70 6.75
Exp. [on Ag(110)]
5
3.88
Exp. [on Ag(111)]
6
3.30.2
Γ
KM
Γ
k
0
100
200
300
400
500
600
ω (cm
-1
)
E
z
= 0
E
z
= 0.51 VÅ
-1
FIG. 2. (Color online) DFT-PBE phonon dispersion curves for
silicene in zero external field and at E
z
= 0.51 V
˚
A
−1
. In both
cases the calculations were performed using the method of finite
displacements, with the atomic displacements being 0.0423
˚
A, in a
supercell consisting of 3 × 3 primitive cells with a 20 × 20 k-point
grid in the primitive cell.
as shown in Fig. 1. The main effect of such an electric field
is to break the symmetry between the A and B sublattices of
silicene’s honeycomb structure and hence to open a gap in
the band structure at the hexagonal BZ points K and K
.Inthe
framework of a simple nearest-neighbor tight-binding model,
this manifests itself in the form of an energy correction to the
on-site energies that is positive for sublattice A and negative for
B. This difference in on-site energies = E
A
− E
B
leads to a
spectrum with a gap for electrons in the vicinity of the corners
of the BZ: E
±
=±
(/2)
2
+|vp|
2
, where p is the electron
“valley” momentum relative to the BZ corner. Opening a gap
in graphene by these means would be impossible because the
A and B sublattices lie in the same plane.
B. First-order perturbation theory
Ana
¨
ıve estimate of the electric-field-induced gap in
silicene can be made using first-order perturbation theory by
diagonalizing a 2 × 2 Hamiltonian matrix at p → 0,
δH(E
z
) = eE
z
ψ
−
K
|z|ψ
−
K
ψ
−
K
|z|ψ
+
K
ψ
+
K
|z|ψ
−
K
ψ
+
K
|z|ψ
+
K
. (1)
Here, ψ
±
K
are the degenerate lowest unoccupied and highest
occupied Kohn-Sham orbitals at the K point at E
z
= 0, and
z = 0 corresponds to the midplane of the buckled lattice. This
suggests a band gap which opens linearly with the electric field
at a rate d/dE
z
= 0.554 and 0.573 e
˚
A for the wave functions
ψ
K
found using the LDA and PBE functionals, respectively.
C. Self-consistent DFT calculations in the presence of the field
The estimate given in Sec. III B is in fact only an upper
limit for the rate at which the band gap opens, since it neglects
screening by the polarization of the A and B sublattices. In
order to obtain an accurate value of the rate at which a band
gap can be opened with an electric field, we have performed
fully self-consistent calculations of the DFT band structure
in the presence of an electric field. A typical result of such
a calculation is shown in Fig. 3(b). At small electric fields,
075423-2
ELECTRICALLY TUNABLE BAND GAP IN SILICENE PHYSICAL REVIEW B 85, 075423 (2012)
Γ
KM
Γ
k
-14
-12
-10
-8
-6
-4
-2
0
ε
(k) (eV)
K
-3.6
-3.4
-3.2
-3.0
ε
(k) (eV)
E
z
= 0
(a)
E
F
Γ
KM
Γ
k
-14
-12
-10
-8
-6
-4
-2
0
ε
(k) (eV)
1.05
1.1
K
k (Å
-1
)
-3.5
-3.4
-3.3
ε
(k) (eV)
E
z
= 0.26 VÅ
-1
(b)
E
F
Γ
KM
Γ
k
-14
-12
-10
-8
-6
-4
-2
0
ε
(
k) (eV)
No reopt.
Reopt.
1.05
1.1
1.15
K
k (Å
-1
)
-4.1
-4.0
-3.9
ε
(
k
) (eV)
E
z
= 0.51 VÅ
-1
(c)
E
F
FIG. 3. (Color online) DFT-PBE band structures for silicene in a
cell of length L
z
= 26.5
˚
A with a plane-wave cutoff energy of 816 eV
and a 53 × 53 k-point grid: (a) in zero external electric field, (b) with
E
z
= 0.26 V
˚
A
−1
, and (c) with E
z
= 0.51 V
˚
A
−1
(shown both with
and without the relaxation of the atomic coordinates in the electric
field). The zero of the external potential is in the center of the silicene
layer. The dashed line shows the Fermi energy in each case and the
insets show the spectrum near the Fermi level in the vicinity of the
K point.
relaxing the structure in the presence of the field does not have
a significant effect on the band gap, but the screening of the
electric potential by the sublattice polarization of the electron
states makes a substantial difference. The DFT-calculated gaps
are gathered in Fig. 4. The variation of the band gap at
K with electric field E
z
is almost perfectly linear for fields
up to E
z
≈ 1V
˚
A
−1
. The results for the rate d/dE
z
at
which a gap is opened are shown in the table inset in Fig. 4.
The eightfold difference between the self-consistent and the
0 0.1 0.2 0.3 0.4
0.5
Electric field E
z
(VÅ
-1
)
0
0.01
0.02
0.03
0.04
0.05
Band gap Δ (eV)
DFT-PBE
DFT-LDA
Method dΔ / dE
z
(eÅ)
0.0742
0.0693
L
z
=13.35 Å
L
z
=18.52 Å
L
z
=26.46 Å
L
z
=∞
L
z
=∞; LDA
Unscreened
FIG. 4. (Color online) DFT gap against applied electric field E
z
for silicene with a plane-wave cutoff energy of 816 eV and a 53 ×
53 k-point grid. Unless otherwise stated, the PBE functional was used.
The box length in the z direction was varied from L
z
= 13.35
˚
Ato
26.46
˚
A. The results have been extrapolated to the limit L
z
→∞of
infinite box length (solid lines) as described in Sec. VD4. Unscreened
band gaps calculated using perturbation theory are also shown. The
inset table shows the calculated rate at which the band gap opens.
unscreened values of d/dE
z
indicates that the system
exhibits a strong sublattice polarizability.
Our value for the rate at which the band gap opens within
DFT-PBE is 0.0742 e
˚
A. This is substantially lower than the
result obtained by Ni et al.,
8
which is 0.157 e
˚
A. Part of the
reason for the discrepancy is that we extrapolated our results
to infinite box length, whereas Ni et al. used a fixed amount
of vacuum between the periodic images of the layers. Another
possible reason for the difference is that we used a plane-wave
basis set, whereas Ni et al. used a localized basis set. An
incomplete localized basis set would tend to undermine the
extent to which the electrons can adjust to screen the electric
field.
D. Stability of the silicene lattice in an electric field
The narrow-gap silicene band structure shown in Fig. 3
persists over a broad range of electric fields E
z
. However,
for electric fields of more than E
z
≈ 0.5V
˚
A
−1
, the band gap
starts to close due to an overlap of the conduction band at and
the valence band at K, and silicene becomes a semimetal, as
shown in Fig. 3(c). According to our calculations, the buckled
honeycomb crystal is still metastable at this electric field, as
can be seen in Fig. 2. The main effects of the electric field on the
phonon dispersion curve are (i) to lift some degeneracies at K
and M and (ii) to soften one of the acoustic branches, but with-
out making the frequency imaginary. Under much higher elec-
tric fields, the honeycomb structure of silicene becomes unsta-
ble. We found that E
z
2.6V
˚
A
−1
causes the lattice parameter
to increase without bound when the structure is relaxed.
IV. SO COUPLING IN SILICENE
A. SO-induced gap
We have also performed a study of the effects of SO
coupling (which is more pronounced in Si than in C) on the
075423-3
N. D. DRUMMOND, V. Z
´
OLYOMI, AND V. I. FAL’KO PHYSICAL REVIEW B 85, 075423 (2012)
Γ
KM
Γ
k
-3
-2
-1
0
1
2
3
4
ε
(k) - E
F
(eV)
PBE
PBE, SO
LDA
LDA, SO
K
-5
0
5
ε
(k) - E
F
(meV)
FIG. 5. (Color online) DFT-PBE and DFT-LDA band structures
with and without SO coupling taken into account. The inset shows
the bands around the K point, revealing a small band gap induced by
SO coupling. The width of the bottom panel corresponds to 1/200 of
the K line.
band structure. The SO coupling term is explicitly included in
the Hamiltonian in the DFT calculations. The results obtained
with the LDA and PBE functionals are shown in Fig. 5.Both
functionals predict an SO gap of the order of a few meV at
the K point, while the rest of the band structure barely differs
from the nonrelativistic case. Our calculated LDA and PBE
SO gaps are 1.4 meV and 1.5 meV, respectively, in agreement
with the recent literature.
20
B. Crossover from topological to band insulating behavior
In the theory of Dirac electrons on the honeycomb lattice,
the SO gap is accounted for by the Kane-Mele term describing,
e.g., intrinsic SO coupling in graphene.
11
The Kane-Mele
SO coupling and the electric-field induced A-B sublattice
asymmetry for electrons in the vicinity of the BZ corners
K
±
= (±4π/(3a),0) in silicene can be incorporated in the
Hamiltonian
H
K
±
= vp · σ +
SO
s
z
σ
z
+
1
2
ξ
z
σ
z
, (2)
where ξ =±1 distinguishes between the two valleys, K
+
and
K
−
, in silicene’s spectrum. Here, the Pauli matrices σ
x
, σ
y
,
and σ
z
act in the space of the electrons’ amplitudes on orbitals
attributed to the A and B sublattices, (ψ
A
,ψ
B
) for the valley
at K
+
and (ψ
B
, −ψ
A
) for the valley at K
−
.InEq.(2), s
z
is the
electron spin operator normal to the silicene plane, and
SO
and
z
are the DFT-calculated SO-coupling and electric-field
induced gaps.
The Hamiltonian of Eq. (2) generically describes the
transition between the 2D topological and band-gap insulators.
Its spectrum,
E
↑±
=±
1
4
(
SO
+ ξ
z
)
2
+ v
2
p
2
,
(3)
E
↓±
=±
1
4
(
SO
− ξ
z
)
2
+ v
2
p
2
,
includes two gapped branches, one with a larger gap |
SO
+
z
| and another with a smaller gap |
SO
−
z
|. At a critical
external electric field E
c
z
≈ 20 mV
˚
A
−1
,
SO
=
z
, and the
smaller gap closes, marking a transition from a topological
insulator
9–11
at
SO
>
z
to a simple band insulator at
SO
<
z
. The difference between these two states of silicene is that
the topological insulator state supports a gapless spectrum of
edge states for the electrons, in contrast to a simple insulator,
where the existence of gapless edge states is not protected by
topology. However, one may expect something reminiscent of
the topological properties of Dirac electrons to show up even
in the band insulator state of silicene: An interface between
two differently gated regions, with electric fields E
z
and −E
z
(where E
z
E
c
z
), should support a one-dimensional gapless
band with an almost linear dispersion of electrons.
21
V. COMPUTATIONAL DETAILS
A. Cohesive energy
All our plane-wave DFT total energies were corrected
for finite-basis error
22
and it was verified that the residual
dependence of the total energy on the plane-wave cutoff energy
is negligible. We used ultrasoft pseudopotentials throughout,
except where otherwise stated. The silicene system was made
artificially periodic in the z direction (normal to the silicene
layer) in our calculations. The atomic structure was obtained
by relaxing the lattice parameter and atom positions within
DFT, subject to the symmetry constraints and at fixed box
length L
z
in the z direction. The cohesive energy was then
evaluated using this optimized structure.
The energy of an isolated Si atom (needed when evaluating
the cohesive energy) was obtained in a cubic box of side-length
L subject to periodic boundary conditions. We extrapolated the
energy of the isolated atom to the limit of infinite box size by
fitting
E(L) = E(∞) + cL
−8
(4)
to the DFT energies E(L) obtained in a range of box sizes,
where E(∞) and c were parameters determined by fitting.
Equation (4) gave a very good fit to our data.
We have also calculated the DFT zero-point correction
to the energy of silicene. This is expected to be largely
independent of the exchange-correlation functional used.
Indeed, our calculations show that the zero-point correction
is 0.103 eV within the LDA and 0.101 eV with the PBE
functional.
12
We used the PBE result in our final calculations
of the cohesive energy reported in Table I.
B. Evaluation of the Fermi velocity
To evaluate the Fermi velocity shown in Table I we
evaluated the DFT band structure using a 53 × 53 k-point
grid and a plane-wave cutoff energy of 816 eV in a cell of
length L
z
= 26.46
˚
A. We then fitted Eq. (17) of Ref. 23 to
the highest occupied and lowest unoccupied bands within a
circular region around the K point; the Fermi velocity is one
of the fitting parameters. The radius of the circular region was
6% of the length of the reciprocal lattice vectors; we verified
that the Fermi velocity was converged with respect to this
radius.
075423-4
ELECTRICALLY TUNABLE BAND GAP IN SILICENE PHYSICAL REVIEW B 85, 075423 (2012)
C. Geometry optimization and phonon
dispersion curves
The phonon dispersion curves shown in Sec. II B were
calculated using the method of finite displacements, with atom
displacements of 0.042
˚
A, in a supercell consisting of 3 × 3
primitive cells with a 20 × 20 k-point grid in the primitive cell.
In the results with the external electric field, the box length was
L
z
= 19.05
˚
A and the plane-wave cutoff energy was 435 eV. In
the results without the field, the box length was L
z
= 13.35
˚
A
and the plane-wave cutoff was 816 eV. This choice was
made because the error due to a finite box length L
z
is
potentially much larger in the presence of a transverse electric
field.
The geometry optimization and band-structure calculations
at zero external field were performed with both the
CASTEP
15,16
and VASP
17
codes, to verify that the results are in good
agreement. This check was necessary because it was only
possible to perform the electric-field calculations with
CASTEP,
while for the SO calculations we had to use
VASP. In principle
the only difference between the calculations performed using
the two codes arises from the Si pseudopotentials used.
The PAW method
24
was used in the VASP calculations,
whereas ultrasoft pseudopotentials were used in the
CASTEP
calculations. As can be seen in Table I, the geometries
predicted by the two codes agree well. We have also verified
that the band structures are in good agreement. Finally, in
Fig. 6 we show that the phonon dispersions obtained with the
two codes are virtually identical when the same parameters are
used.
Figure 6 also demonstrates that our phonon dispersion
curves are converged with respect to supercell size.
Γ
KM
Γ
k
0
100
200
300
400
500
600
ω (cm
-1
)
CASTEP (LDA, 3×3 cell)
CASTEP (PBE, 3×3 cell)
VASP (LDA, 3×3 cell)
VASP (LDA, 7×7 cell)
FIG. 6. (Color online) Phonon dispersion curves for silicene
obtained with
CASTEP and VA S P using different exchange-correlation
functionals and supercell sizes. The results for a 3 × 3 supercell were
obtained with a box length L
z
= 13.35
˚
A, a 20 × 20 k-point grid in
the primitive cell, and a plane-wave cutoff energy of 816 eV. The
matrix of force constants was evaluated using the method of finite
displacements, with the displacements being 0.042
˚
A. The results for
a7× 7 supercell were obtained with a box length L
z
= 15.0
˚
A, a
12 × 12 k-point grid in the supercell, and a plane-wave cutoff energy
of 500 eV. The matrix of force constants was evaluated using the
method of finite displacements, with displacements of 0.09
˚
A.
400
600
800 1000 1200
PW cutoff ener
gy
(
eV
)
18.568
18.570
18.572
18.574
18.576
18.578
18.580
Band gap (meV)
FIG. 7. (Color online) DFT-PBE gap against plane-wave (PW)
cutoff energy for silicene subject to an electric field of 0.257 V/
˚
Ain
a cell of length L
z
= 13.35
˚
Awitha15× 15 k-point grid including K.
D. Band gap in the presence of an external electric field
1. Plane-wave cutoff energy
The convergence of the calculated band gap with respect to
the plane-wave cutoff energy for a particular applied field is
shown in Fig. 7. The gap converges extremely rapidly.
2. k-point sampling
The convergence of the calculated band gap at BZ point
K with respect to the k-point grid used in the self-consistent
field calculations is shown in Fig. 8. The finite-sampling error
falls off as the reciprocal of the total number of k points. The
prefactor of the finite-sampling error is vastly greater when K
or K
is included in the grid of k points for the self-consistent
field calculations.
3. Choice of pseudopotential
The dependence of the calculated gap on the exchange-
correlation functional and pseudopotential is shown in Fig. 9.
The difference between the results obtained with different
pseudopotentials is much smaller than the gap, but is not
0
0.005
0.01
1 / N
k
2
0
5
10
15
20
25
Band gap (meV)
K not included
K included
FIG. 8. (Color online) DFT-PBE gap against the reciprocal of the
number N
2
k
of k points for silicene subject to an electric field of
0.257 V/
˚
A in a cell of length L
z
= 13.35
˚
A with a plane-wave cutoff
energy of 816 eV.
075423-5