BHZ model

This example shows the calculation of the \(\mathbb{Z}_2\) invariant for the periodic BHZ model.

import logging

import matplotlib.pyplot as plt
import numpy as np
import z2pack

logging.getLogger("z2pack").setLevel(logging.WARNING)

pauli_y = np.array([[0, 1], [1, 0]], dtype=complex)
pauli_x = np.array([[0, -1j], [1j, 0]])
pauli_z = np.diag([1, -1])
pauli_vector = [pauli_x, pauli_y, pauli_z]


def bhz(A, B, C, D, M):
    def h(k):
        kx, ky = 2 * np.pi * np.array(k)
        d = [
            A * np.sin(kx),
            -A * np.sin(ky),
            -2 * B * (2 - (M / (2 * B)) - np.cos(kx) - np.cos(ky)),
        ]
        H = sum(di * pi for di, pi in zip(d, pauli_vector))
        epsilon = C - 2 * D * (2 - np.cos(kx) - np.cos(ky))
        H += epsilon
        return H

    def Hamiltonian(k):
        return np.vstack(
            [
                np.hstack([h(k), np.zeros((2, 2))]),
                np.hstack([np.zeros((2, 2)), h(-np.array(k)).conjugate()]),
            ]
        )

    return Hamiltonian


def run(**kwargs):
    system = z2pack.hm.System(bhz(**kwargs), dim=2)
    res = z2pack.surface.run(system=system, surface=lambda s, t: [s / 2, t])
    z2pack.plot.wcc(res)
    plt.savefig("wcc_plot.pdf", bbox_inches="tight")
    print("Z2 invariant:", z2pack.invariant.z2(res))


if __name__ == "__main__":
    run(A=0.5, B=1.0, C=0.0, D=0.0, M=1.0)