# 3D 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.

Here we discuss the calculation of the topological indices for 3d topological insulators. The algorithm is based on Wannier centers as in 2D topological insulator. The only difference is that we need to calculate four \(Z_2\) invariants for 4 planes in the 3d BZ. In principle there are 6 different planes forming the faces of a parallelepiped whose corners are the eight distinct time reversal invariant momenta (TRIM) of any chosen smallest primitive reciprocal cell. However, if the electronic spectrum is gapped above a certain band there are only 4 independent \(Z_2\) invariants, which are usually grouped into four topological indices \(\nu_0;(\nu_1\nu_2\nu_3)\). \(\nu_0=1\) indicates a strong topological insulator while any \(\nu_{1,2,3}=1\) indicate conditions on the number of surface state pairs for certain surface orientations. If \(\nu_0=0\) but some \(\nu_{1,2,3}\neq 0\) it is a weak topological insulator. If all indices are zero it is a trivial insulator. The mapping of \(Z_2\) invariants to indices is as follows: if any two parallel planes of the parallelepiped have differing invariants \(\nu_0=1\). We take the two planes spanned by the first two reciprocal lattice vectors (\(g_{1,2}\)) called \(z_0\) if the plane goes through the origin and \(z_1\) if it goes through \(g_{3}\). The three weak indices are identical to the \(Z_2\) invariants of the three planes, which do not go through the origin, which we call \(x_1\) (spanned by \(g_{2,3}\) through \(g_1\)), \(y_1\) (spanned by \(g_{3,1}\) through \(g_2\)) and \(z_1\) (spanned by \(g_{1,2}\) through \(g_3\)). The weak indices of course depend on the chosen TRIM points.

The algorithm selects the TRIM points and performs four \(Z_2\) calculations. It tries to determine the invariants via the automatic procedure explained in 2D topological insulator. The results including reliability and electronic gap estimates are printed to the output and convenience files `Z2_3dTI_homo...cmd` for use in xfbp are created.

This algorithm works for centro symmetric and non-centro symmetric compounds. A version of this algorithm is linked directly into fplo. The difference is that when used from fplo a all-orbital Wannier basis is created instead of the reduces Wannier models, which usually are the basis of `pyfplo.slabify`. This leads to differently locking Wannier center curves, since there are more orbitals (semi core and such). If surface state calculations are planned it is anyway necessary to create a Wannier model first. The resulting \(Z_2\)-calculations are faster than in fplo.

For illustration we chose a system, which is centro symmetric, since in this case we know the invariants from the parity algorithm, which allows comparison. This example also shows some fo the complications. We will use both the fplo version and the `pyfplo` version of the Wannier center algorithm to understand the pecularities.

The tutorial files are in `FPLO.../DOC/pyfplo/Examples/slabify/3dTI/Bi2Se3`.

To start with change into this directory and have a look into the topological insulator submenu in fedit. You will see that “Force wannier center” is switched on, to force the use of this algorithm despite the fact that this compound has an inversion center.

Now, have a look at `makewandef.py`. It illustrates how to create a `=.wandef` quickly. There is a particularity in the system, which is that the Bi6p and Se4p orbitals are separated from lower and higher bands through two gaps. This makes creating a Wannier model especially easy. We use the option lbands and ubands to define an energy window which excludes all other bands. Please execute

    python makewandef.py

You should have a =.wandef file now. Next run fplo. E.g.:

    fplo.... > out

to produce a bandstructure, `+wancoeff` and the output from the fplo Wannier center calculation. This will take some time, since the fplo version of the WC algorithm is not the fastest, due to the larger number of bands needed. Output related to the TI calculation can be found via:

    grep TI: out

First, the invariants from parities are printed and then the Wannier center algorithm is executed. Various files are created. Note, that the indices from parities can differ from the Wannier center results; first, because the WC algorithm is tricky or second, if the electronic gap for a certain homo is zero, in which case both algorithms makes no sense.

You will see the following output from the Wannier centers:

    TI: Invariants:
    TI:   homo  invariants  E(k=0)                estim.gaps [eV] for
    TI:                                  z0        z1        x1        y1
    TI:     114  ? 0;(000)  -1.09567   0.05585  0.01917  0.01917  0.01917
    TI:     116  ? 1;(111)  -0.79835   0.01306  0.03873  0.03873  0.03873
    TI:     118 *  1;(000)  -0.22789   0.44699  0.43513  0.43513  0.43513
    TI:     120  ? 1;(111)   0.28647   0.15850  0.14845  0.14845  0.14845
    TI:     122  ? 0;(111)   1.42635   0.02117  0.17378  0.17378  0.17378
    TI:     124  ? 0;(000)   1.44752   0.21087  0.29692  0.29692  0.29692
    TI:
    TI:   homo  invariants  E(k=0)                  Z2 (reliability)
    TI:                                 z0        z1        x1        y1
    TI:     114  ? 0;(000)  -1.09567    0 ( 77%)  0 ( 77%)  0 ( 77%)  0 ( 77%)
    TI:     116  ? 1;(111)  -0.79835    0 ( 50%)  1 ( 73%)  1 ( 73%)  1 ( 73%)
    TI:     118 *  1;(000)  -0.22789    1 (100%)  0 (100%)  0 (100%)  0 (100%)
    TI:     120  ? 1;(111)   0.28647    0 ( 92%)  1 ( 65%)  1 ( 65%)  1 ( 65%)
    TI:     122  ? 0;(111)   1.42635    1 ( 69%)  1 (100%)  1 (100%)  1 (100%)
    TI:     124  ? 0;(000)   1.44752    0 ( 80%)  0 (100%)  0 (100%)  0 (100%)

The first table shows homos, indices the energy at the Gamma point and estimates of the electronic gap (taken from the values on the finite grid used for the WC algorithm), while the second shows the \(Z_2\) invariants and their estimated reliability for the four planes of the parallelepiped mentioned above. If the reliability is less than 100% a question mark is printed after the homo. The homo, which is likely the highest occupied band is marked by a * for orientation. This is of course only reliable if the system has a true gap.

In our case we see that there are rather small estimated gaps for lots of bands. Fortunately above the Fermi energy (homo 118) the gap is sizable. Not, surprisingly the reliabilities are all very high (actually 100%). We also find lots of question marks. Let’s have a look at the fplo and Wannier model bands (red/green).

Fig. 24 fplo and Wannier model bands (red/green)

We see the two gaps, separating the p-bands as discussed above and that nearly all bands are separated by gaps (although small). This means that topological invariants have a meaning. Note, that it is not needed that the gap is visible in the DOS. All one needs is that the bands of homo and homo+1 never cross (warped gap).

Since this table is not ment for straight forward consumption we have to inspect the Wannier centers by hand. Run:

    xfbp Z2_3dTI_homo118.cmd

and so forth. You will get several pictures. Let’s start with Wannier centers and reference line for homo 118. As already pointed out in the table his band has a nice electronic gap above itself, happens to be the last occupied band below the Fermi level and has high reliability of the automatically determined \(Z_2\) invariants (reference line winding number).

Fig. 25 Wannier centers and reference line for homo 118

It can be clearly seen that the \(x_1\), \(y_1\) and \(z_1\) planes have zero \(Z_2\) invariant, while plane \(z_0\) is non-trivial. To really verify this use the right mouse click in xfbp close to the higest Wannier center curve in the \(z_0\) panel close to \(\theta=1\) and \(k_y=0\). It will show that only one center (called Set118 in xfbp) comes down from it’s degenerate value of +-1. A clean straight reference line can be drawn for \(\theta=0.9\) (not shown), which only crosses this center and hence crosses an odd number of Wannier centers, which results in \(Z_2=1\). This and the fact that we can visually connect the Wannier center curves in a reasonable smooth way convinces us that the topological indices are 1;(000). (Remember \(z_0\neq z_1 \to \nu_0=1\) and \(\nu_{1}=Z_2(x_1), \ldots\).) In fact due to symmetry the tree planes \(x_1\), \(y_1\) and \(z_1\) give the same result in this compound.

Next we have a look at all the other homos, just for learning purposes.

Fig. 26 Wannier centers and reference line for homo 114

The Wannier centers and reference line for homo 114 are quite messy. However, you should convince yourselfs that \(z_1\) has a region (\(\theta\approx0.5\)) without any centers crossing \(\to\) trivial, while for the \(z_0\)-plane an even number of curves cross any reference line. Hence we get 0;(000).

Fig. 27 Wannier centers and reference line for homo 116

The Wannier centers and reference line for homo 116 are also messy and reveal another possible complication. First, note that both 114 and 116 yield curves with fast varying (steep slopes) curves. These steep curves generate only very few points on their steep section. Depending on the density of the ky mesh there might even be zero points on such sections, which would make it look as if certain curves end for some ky value to reemerge for a later ky-value. Any algorithm including our own judgment can fail in such cases. All we can do is to take a denser ky-mesh. In this case the automatic reference line algorithm is correct as can be visually verified. There is an even/odd number of curves crossed by a suitably chosen reference curve for the \(z_0\)/\(z_1\) planes (\(x_1,y_1=z_1\)). Hence we get 1;(111). Note, however, that the electronic gaps above homos 114 and 116 are rather small, they could well be zero. If this homo were important we would try to use much finer meshes to make absolutely sure that the gap is finite.

Homo 122 also has a small electronic gap. If we assume that it is finite we can proceed analysing the invariants. The Wannier centers and reference line for homo 120 and the Wannier centers and reference line for homo 122 both show that the indices printed in the output table are wrong (as can be verified in this case due to the parity algorithm).

Fig. 28 Wannier centers and reference line for homo 120

For homo 120 the \(z_0\) plane is clearly trivial since we can draw a non crossing reference line at \(\theta=\pm1\) (the automatic blue reference line is just that). For the other three planes the algorithm made a mistake around ky=0.33, where a single red dot appears. This mistake is due to a steep section at one curve. One can find hand drawn reference curves, which cross an even number of times and hence these planes are trival. The result should be 0;(000). Please also note that the reliability for homo 120 is high for the \(z_0\) plane but only 65% for the other three, which fits to the findings explained above. This is of course only an indicator.

Fig. 29 Wannier centers and reference line for homo 122

For homo 122 the situation is even trickier. According to the reliabilities the \(z_0\) panel might have a problem. Indeed, the \(z_1\) plane (and \(x_1\), \(y_1\)) is clearly topological, but the \(z_0\) plane was determined as non trivial although it is not. Again due to small electronic gaps in the \(z_0\) plane we get very steep curves close to ky=0. In fact this curve is so steep that we see only one point in the steep section. In xfbp use the right mouse click at the dot in the upper left corner. It will show that there are two curves. Given the rest of the dots it is clear that one dot is connected with the declining curve and the other must fall very steeply to connect to a curve close to \(\theta=0.28\). Hence the reference line actually crosses two curves and not one. Consequently, \(z_0\neq z_1 \to\) 1;(111) in contrast with the result printed in the output table.

Homo 124 is again a simple case and will not be discussed here.

In summary we saw that the fplo based version of the procedure is rather slow, which is bad if a high grid density is needed to obtain reliable results. Let’s continue with the `pyfplo` version then. We already set up the Wannier basis production in our first fplo run. Now, we can proceed by re-running:

    fplo...

Note that we do not use an output file (we don’t want to overright the TI results). This will produce among others the file `+hamdata`, which we will use for the `pyfplo` based version of the TI procedure. If the fplo run is over we confirm the Wannier basis by checking:

    xfbp wband.cmd

which should show green bands (finite cut-off tight-binding Wannier bands) on top of red bands (exact Wannier bands) on top of black bands (fplo bands). Please, change into `slabify/Z2` and execute:

    python -u Z2.py | tee out

The flag `-u` un-buffers the output, such that the progress is updated in real time. The file’s content is

     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
    70
    71
    72
    73
    #! /usr/bin/env python
    from __future__ import print_function
    import sys,copy
    import numpy as np
    import pyfplo.slabify as sla
    import pyfplo.common as com
    import pyfplo.fploio as fploio

    # ===================================================================
    #
    # ===================================================================
    def work():

        print( '\nversion: ',sla.version,'\n')

        np.set_printoptions(precision=5,suppress=True,linewidth=120)

        hamdata='../../+hamdata'

        s=sla.Slabify()
        s.object='3d'
        s.printStructureSettings()
        s.prepare(hamdata)

        p=fploio.INParser()
        p.parseFile('../../=.in')
        d=p()('special_sympoints')
        l=[]
        for i in range(d.size()):
            l.append([ d[i]('label').S,d[i]('kpoint').listD])

        bp=com.BandPlot()
        bp.points=l
        bp.ndiv=100
        bp.calculateBandPlotMesh(s.dirname)

        #help(sla)
        s.calculateBandStructure(bp)

        # 14 -> 0;(000)
        # 16 -> 1;(111)
        # 18 -> 1;(000)
        # 20 -> 0;(000)
        # 22 -> 1;(111)
        # 24 -> 0;(000)

        s.calculate3dTIInvariants(200,200,homos=range(14,25,2),efhomo=18)

        return

    # ===================================================================
    #
    # ===================================================================

    if __name__ == '__main__':

        work()

In lines 34-47 we calculate the band structure of the Wannier model along the same high symmetry points as were used in the fplo calculation. With the help of:

    xfbp slabifyres/+band_sf

and a right mouse click we find out that the Fermi level is above homo 18 (just a coincidence that 118-18=100). This is used in line 56 to define a set of homos and efhomo. The grid is chosen finer for both the “integration” and the ky direction. See `calculate3dTIInvariants`.

The output table:

    Invariants:
      homo  invariants  E(k=0)                estim.gaps [eV] for
                                     z0        z1        x1        y1
         14  ? 0;(000)  -1.09645   0.02338  0.01563  0.01563  0.01563
         16  ? 1;(111)  -0.79820   0.00531  0.02098  0.02098  0.02098 <-small
         18 *  1;(000)  -0.22063   0.42925  0.43611  0.43611  0.43611
         20    0;(000)   0.27854   0.14357  0.06937  0.06937  0.06937
         22    1;(111)   1.42038   0.00911  0.17176  0.17176  0.17176 <-small
         24    0;(000)   1.42948   0.20549  0.24304  0.24304  0.24304

      homo  invariants  E(k=0)                  Z2 (reliability)
                                    z0        z1        x1        y1
         14  ? 0;(000)  -1.09645    0 (100%)  0 ( 93%)  0 ( 93%)  0 ( 93%)
         16  ? 1;(111)  -0.79820    0 (100%)  1 ( 86%)  1 ( 93%)  1 ( 93%)
         18 *  1;(000)  -0.22063    1 (100%)  0 (100%)  0 (100%)  0 (100%)
         20    0;(000)   0.27854    0 (100%)  0 (100%)  0 (100%)  0 (100%)
         22    1;(111)   1.42038    0 (100%)  1 (100%)  1 (100%)  1 (100%)
         24    0;(000)   1.42948    0 (100%)  0 (100%)  0 (100%)  0 (100%)

looks a bit different from the previous one. First of all the gap values are different since there are always tiny shifts in a reduced Wannier model compared to a full band structure and second we used a finer grid. Furthermore, the reliabilities increased considerably, again mostly due to the finer grid. Finally, all topological indices are correct.

We can now load all the `Z2_3dTI_homo...cmd` files into xfbp and analyse the validity by hand. We will look at three Wannier center graphs. The Wannier centers and reference line for homo 16 show a quick change in the \(z_0\)-panel around ky=0.16, which is most likely due to the small electronic gap. Furthermore, the steep sections in the \(z_1\)-panel are now represented by more points, which makes the automatic reference line algorithm more reliable.

Fig. 30 Wannier centers and reference line for homo 16

The Wannier centers and reference line for homo 20 has the correct reference line due to the increased mesh size.

Fig. 31 Wannier centers and reference line for homo 20

The Wannier centers and reference line for homo 22 now clearly shows that two WC curves meet close to the upper left corner of the \(z_0\)-panel, which makes this plane trivial. Consequently, the invariants are correctly determined as 1;(111).

Fig. 32 Wannier centers and reference line for homo 22

In summary, the use of a reduce Wannier function model leads to a much faster TI procedure, which allows to use finer grids, which in turn increase the accuracy of the results. Additionally, one could calculate surface spectra using the same model.