2D topological insulator¶
Note
You need to use the newer xfbp/xfplo version,
which comes with pyfplo
in order for the cmd
scripts to work
properly.
This example tries to explain how to calculate the \(Z_2\) invariant of the
2d time reversal invariant topological insulator BHZ model (Bernevig,
Hughes and Zhang). The tutorial files are in
FPLO.../DOC/pyfplo/Examples/slabify/BHZmodel
where FPLO...
stands for your version’s FPLO directory, e.g. FPLO22.00-62
.
The method can be applied to any plane spanned by time reversal invariant points also in 3d lattices, e.g. in Weyl semi metals, where the \(Z_2\) invariant tells us how many pairs of surface states are to be expected to cross a given line – the projection of the plane onto the surface. (All assuming that the corresponding plane is gapped.)
We use a Hamiltonian file (+hamdata
) containing the model
data, which was created by bhz.py
. Usually +hamdata
is created by the Wannier function module of fplo. The model data
are taken from [Yu2011] and modified a bit (a=5, M_TI=-6). The model
defines a four orbital tight binding model. Depending on the parameter
choice the model is either a trivial insulator or a topological
insulator. The two cases are in subdiretories called I
and
TI
.
[Yu2011] | Rui Yu et.al. Phys. Rev. B 84, 075119 (2011). https://doi.org/10.1103/PhysRevB.84.075119 |
The topological phase¶
In python script for the 2d band structure and the \(Z_2\) calculation is
TI/2d/Z2.py
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 | #! /usr/bin/env python
from __future__ import print_function
import sys
# If your pyfplo is not found you could also
# explicitly specify the pyfplo version path:
#sys.path.insert(0,"/home/magru/FPLO/FPLO22.00-62/PYTHON/doc");
import numpy as np
import pyfplo.slabify as sla
print( '\npyfplo version=: {0}\nfrom: {1}\n'.format(sla.version,sla.__file__))
# protect against wrong version
#if fedit.version!='22.00': raise RuntimeError('pyfplo version is incorrect.')
# ===================================================================
#
# ===================================================================
def work():
hamdata='../+hamdata'
s=sla.Slabify()
s.object='3d'
s.printStructureSettings()
s.prepare(hamdata)
# prepare BandPlot
bp=sla.BandPlot()
bp.ndiv=100
bp.points=[
['$~G',[0,0,0]],
['X',[0.5,0,0]],
['M',[0.5,0.5,0]],
['$~G',[0,0,0]],
['Y',[0,0.5,0]],
]
bp.calculateBandPlotMesh(s.dirname)
s.calculateBandStructure(bp)
# wannier centers
s.calculateZ2Invariant([0,0,0],[0.5,0,0],[0,0.5,0],Nint=20,Nky=100,
homos=[2,4])
# ===================================================================
#
# ===================================================================
if __name__ == '__main__':
work()
|
Line 8 can be uncommented and edited to make python search for
pyfplo in a particular location (also see Setup).
Line 14 shows which pyfplo
version was loaded and from where.
If you uncomment line 16, the script is showing an error message if the
pyfplo
version does not match a particular version.
Please run the script (from TI/2d/README.rst
):
# run the script as in
./Z2.py | tee out
xfbp bands.cmd
xfbp Z2_homo2.cmd
xfbp Z2_homo4.cmd
# and have a look at the output and into slabifyres/
In line 27 we set the proper +hamdata
and in line 30 we set
the object to '3d'
, which in the case of a 2d lattice means 2d-bulk.
In lines 35-48 the band structure is calculated
(2d bulk bands for the topological phase.). Note that there is a small gap close to the
\(\Gamma\)-point where the orbital weight is inverted across the
gap which already indicates the possibility of a topological phase.
In line 52 the Wannier centers are calculated. Consult
calculateZ2Invariant
for more information. Note, that we have to
chose three time reversal invariant points (TRIM) in the Brillouin
zone in units of kscale
. At the end of the output you will see
the primitive reciprocal cell and the kscale factor (code lines
57-58). The TRIM points are always in the middle between two
reciprocal cell vectors. This should explain our choice.
The output contains
quite a lot of information. Please consider it all. For each homo
which was specified a file containing the Wannier centers
slabifyres/+Z2_homo...
and a file containing a reference line
+zgap_homo...
are created. The reference line usually follows
the largest gap between Wannier centers [Sol2011]. In our
modification a number of such lines for some of the largest WC-gaps is
created. Each yields a topological index. The majority result will
win. The percentage of such lines with the majority result is called
reliability in the output. The table of invariants for all homos in
the output looks like
[Sol2011] | A. A. Soluyanov, D. Vanderbilt Phys. Rev. B 83, 235401 (2011) |
Invariants:
homo Z2 (reliability) smallest gap [eV]
2 1 (100%) 0.20353
4 0 (100%) >10
If the reliability is 100% all WC-gap-following curves gave the same result. This does however not mean it is a save result. The exactness of this algorithm depends on the actual system. If the Berry curvature is strongly fluctuating a finer grid is needed. This often happens e.g. if the energy gaps are small or zero. In such cases the degeneracies at \(ky=0\) and \(ky=\pi\) can be broken. Also discontinuities can occur and the results can depend on the evenness/oddness of Nint. All these are warning signs that the \(Z_2\) invariant might be nonsensical.
Additionally, the smallest energy gap for each homo on the grid defined by the parameters Nint and Nky is printed. This helps you decide, if the time reversal invariant plane is really gaped. Also this information is not fool proof, since the grid is discrete.
If efhomo is given to calculateZ2Invariant
a * is printed after
this homo number, usually ment to be the Fermi level. If the
reliability is less then 99% a question mark is printed before the \(Z_2\)
invariant, as in in the following example:
Invariants:
homo Z2 (reliability) smallest gap [eV]
12 ? 1 (90%) 0.4
14 *? 0 (67%) 0.13
16 1 (100%) 0.56
Let’s have a look at the corresponding Wannier center files
(Wannier centers for homo 2 in the topological phase. and Wannier centers for homo 4 in the topological phase.) created by
loading Z2_homo...cmd
into xfbp.
The Soluyanov algorithm draws a curve through all wannier centers,
such that it always stays in the largest gap between Wannier centers
(which we modified to give more competing solutions). If the number of
centers crossed so far at a certain ky-value is even a blue symbol is
plotted if the number is odd a red symbol is plotted. If the last
symbol (\(ky=\pi\)) is red the \(Z_2\) invariant is odd. Of course the
WC-curves are discrete and hence we cannot a priori decide about their
connectedness, which leaves an uncertainty. Ultimately, the user has
to look at these pictures and decide for himself, how many WCs a
chosen reference line crosses. In xfbp a right mouse click close
to a curve tells you how many sets there are, which helps to verify
the number of crossed Wannier centers. Any reference curve (also a
straight horizontal line) is OK. The one shown in the plots is just
the automatically determined one. Since we modified the algorithm, the
automatic curve does not necessarily follow the largest WC-gaps. The
two figures tell us that homo 2 has a non trivial \(Z_2\) invariant in the
plane spanned by the TRIM points given to calculateZ2Invariant
.
Now, we make a calculation of the 1d surface states. Go into
TI/semi
and run the script (from TI/semi/README.rst
):
# run the script as in
./calcedc.py
xfbp edc.cmd
# and have a look at the output and into slabifyres/
The script looks like this:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 | #! /usr/bin/env python
from __future__ import print_function
import sys
# If your pyfplo is not found you could also
# explicitly specify the pyfplo version path:
#sys.path.insert(0,"/home/magru/FPLO/FPLO22.00-62/PYTHON/doc");
import numpy as np
import pyfplo.slabify as sla
print( '\npyfplo version=: {0}\nfrom: {1}\n'.format(sla.version,sla.__file__))
# protect against wrong version
#if fedit.version!='22.00': raise RuntimeError('pyfplo version is incorrect.')
# ===================================================================
#
# ===================================================================
def work():
hamdata='../+hamdata'
s=sla.Slabify()
s.object='semislab'
s.zaxis=[0,1,0]
s.numberoflayers=2
s.printStructureSettings()
s.prepare(hamdata)
bp=sla.BandPlot()
bp.ndiv=100
bp.points=[
['X',[0.5,0,0]],
['$~G',[0,0,0]],
['X',[0.5,0,0]],
]
bp.calculateBandPlotMesh(s.dirname)
ec=sla.EnergyContour(400,-6,6)
s.calculateEDC(bp,ec,penetrationdepth=-1,suffix='_1block')
# another zoom
s.dirname='zoom'
bp.points=[
['X $arrowleft',[0.15,0,0]],
['$~G',[0,0,0]],
['$arrowright X',[0.15,0,0]],
]
bp.calculateBandPlotMesh(s.dirname)
ec=sla.EnergyContour(200,-0.5,0.5)
s.calculateEDC(bp,ec,penetrationdepth=-1,suffix='_1block')
# ===================================================================
#
# ===================================================================
if __name__ == '__main__':
work()
|
Lines 30-32 setup a semislab
with [010] surface and two layers in
the primary layer. Since the bulk is 2d with a [001] surface this
results in a semi infinite ribbon. (Actually, in the moment the 2d
bulk is a 3d repeated monolayer slab.) After setup, the former [010]
direction is the new c-direction (surface of the semi infinite ribbon)
but still the cartesian y-axis (remember we use rotations internally
to keep the original orientation). The global x-axis is now the
(periodic) b-direction, while the global z-axis (as for the 2d bulk)
is the direction perpendicular to the ribbon. If we chose the high
symmetry points along the cartesian x-axis (lines 40-42, 52-54) we
sample k-points in the periodic b-direction.
In lines 37-46 we set up a bandplot path and an energy contour and
calculate the energy resolved Bloch spectral density. In line 50-59 a
zoomed in version is calculated and placed into a different directory
(line 50), because the +points
file is different for both
calculations.
The result looks like this Figure and clearly shows a Dirac cone at the edge of the surface.
The trivial insulator phase¶
In python script for the 2d band structure and the \(Z_2\) calculation is
I/2d/Z2.py
and is identical to the one for the topological
phase.
Please run the script (from I/2d/README.rst
):
# run the script as in
./Z2.py | tee out
xfbp bands.cmd
xfbp Z2_homo2.cmd
xfbp Z2_homo4.cmd
# and have a look at the output and into slabifyres/
Have a look at the 2d bulk bands for the trivial phase.
It looks pretty much like the 2d bulk bands for the topological phase. except that the
band inversion across the gap is missing. Consequently,
Wannier centers for homo 2 in the trivial phase. (created by xfbp Z2_homo2.cmd
) clearly
indicates a trivial phase as is also seen by the output:
Invariants:
homo Z2 (reliability) smallest gap [eV]
2 0 (100%) 0.46419
4 0 (100%) >10
As a final confirmation let’s have a look at the Surface states in the trivial phase..
Run the script (from TI/semi/README.rst
):
# run the script as in
./calcedc.py
xfbp edc.cmd
# and have a look at the output and into slabifyres/
The python script is the same as for the topological phase.