# 2D topological insulator

- The topological phase
- The trivial insulator phase

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.

Fig. 16 2d bulk bands for the 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) https://doi.org/10.1103/PhysRevB.83.235401

    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.

Fig. 17 Wannier centers for homo 2 in the topological phase.

Fig. 18 Wannier centers for homo 4 in the topological phase.

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.

Fig. 19 Surface states in the topological phase.

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.

Fig. 20 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

Fig. 21 Wannier centers for homo 2 in the trivial phase.

Fig. 22 Wannier centers for homo 4 in the trivial phase.

As a final confirmation let’s have a look at the Surface states in the trivial phase..

Fig. 23 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.