In [1]:
import sys

sys.path.append("../../code")
from init_mooc_nb import *

init_notebook()


def bhz(X=None, Y=None, Z=None, system_type="infinite"):
    """A cuboid region of BZZ material with two leads attached.

    parameters for leads and scattering region can be defined separately
    """
    # Onsite and hoppings matrices used for building BZZ model
    def onsite(site, p):
        return (p.C + 2 * p.D1 + 4 * p.D2) * pauli.s0s0 + (
            p.M + 2 * p.B1 + 4 * p.B2
        ) * pauli.s0sz

    def hopx(site1, site2, p):
        return -p.D2 * pauli.s0s0 - p.B2 * pauli.s0sz + p.A2 * 0.5j * pauli.sxsx

    def hopy(site1, site2, p):
        return -p.D2 * pauli.s0s0 - p.B2 * pauli.s0sz + p.A2 * 0.5j * pauli.sysx

    def hopz(site1, site2, p):
        return -p.D1 * pauli.s0s0 - p.B1 * pauli.s0sz + p.A1 * 0.5j * pauli.szsx

    def hopx_phase(site1, site2, p):
        x1, y1, z1 = site1.pos
        x2, y2, z2 = site2.pos
        return hopx(site1, site2, p) * np.exp(-0.5j * p.Bz * (x1 - x2) * (y1 + y2))

    def shape_slab(pos):
        (x, y, z) = pos
        return 0 <= z < Z

    def shape_lead(pos):
        (x, y, z) = pos
        return (0 <= y < Y) and (0 <= z < Z)

    def shape_cube(pos):
        (x, y, z) = pos
        return (0 <= x < X) and (0 <= y < Y) and (0 <= z < Z)

    lat = kwant.lattice.general(np.identity(3), norbs=4)

    if system_type == "slab":
        syst = kwant.Builder(kwant.TranslationalSymmetry([1, 0, 0], [0, 1, 0]))
        syst[lat.shape(shape_slab, (0, 0, 0))] = onsite
    if system_type == "lead":
        syst = kwant.Builder(kwant.TranslationalSymmetry((1, 0, 0)))
        syst[lat.shape(shape_lead, (0, 0, 0))] = onsite
    elif system_type == "cuboid":
        syst = kwant.Builder()
        syst[lat.shape(shape_cube, (0, 0, 0))] = onsite
    elif system_type == "infinite":
        syst = kwant.Builder(kwant.TranslationalSymmetry(*lat.prim_vecs))
        syst[lat.shape(lambda pos: True, (0, 0, 0))] = onsite

    syst[kwant.HoppingKind((1, 0, 0), lat)] = hopx_phase
    syst[kwant.HoppingKind((0, 1, 0), lat)] = hopy
    syst[kwant.HoppingKind((0, 0, 1), lat)] = hopz
    return syst


def title(p):
    return r"$M={:.3}$".format(p.M)


def make_lead():
    lat = kwant.lattice.general(np.identity(3), norbs=4)
    syst = kwant.Builder(
        kwant.TranslationalSymmetry((-1, 0, 0)), time_reversal=1j * pauli.sys0
    )
    syst[lat(0, 0, 0)] = 1.5 * pauli.s0sz
    syst[kwant.HoppingKind((-1, 0, 0), lat)] = -1 * pauli.s0sz
    return syst


def make_scatter_sys():
    syst = kwant.wraparound.wraparound(bhz(Z=1, system_type="slab"))
    syst.attach_lead(make_lead())
    syst.attach_lead(kwant.wraparound.wraparound(bhz(), keep=0))
    syst = syst.finalized()
    return syst


def scattering_det_pfaff(syst, p):
    def pfaffian(syst, p, k_x, k_y):
        smat = kwant.smatrix(
            syst, energy=0.0, params=dict(p=p, k_x=k_x, k_y=k_y, k_z=0)
        ).data
        # since we get relatively large numerical errors we project the matrix on
        # the space of antisymmetric matrices
        smat = 0.5 * (smat - smat.T)
        return pf.pfaffian(smat)

    xdim, ydim = "$k_y$", "phase"

    def plot_k_x(syst, p, k_x, label, col):
        pfaff = [pfaffian(syst, p, k_x, 0), pfaffian(syst, p, k_x, np.pi)]
        ks = np.linspace(0.0, np.pi, 50)
        det = [
            np.linalg.det(
                kwant.smatrix(
                    syst, energy=0.0, params=dict(p=p, k_x=k_x, k_y=k_y, k_z=0)
                ).data
            )
            for k_y in ks
        ]
        det = np.array(det)
        phase = np.angle(pfaff[0]) + 0.5 * np.cumsum(np.angle(det[1:] / det[:-1]))
        kdims = [xdim, ydim]
        plot = holoviews.Path((ks[1:], phase), kdims=kdims, label=label).opts(
            style={"color": col}
        )
        plot *= holoviews.Points(([0, np.pi], np.angle(pfaff)), kdims=kdims).opts(
            style={"color": col}
        )
        return plot

    plot = plot_k_x(syst, p, 0, r"$k_x=0$", "g") * plot_k_x(
        syst, p, np.pi, r"$k_x=\pi$", "b"
    )
    xlims, ylims = (-0.2, np.pi + 0.2), (-np.pi - 0.2, np.pi + 0.2)
    pi_ticks = [(-np.pi, r"$-\pi$"), (0, "$0$"), (np.pi, r"$\pi$")]
    style_overlay = {
        "xticks": [(0, "0"), (np.pi, "$\pi$")],
        "yticks": pi_ticks,
        "show_legend": True,
        "legend_position": "top",
    }
    style_path = {"show_legend": True}
    return plot.redim.range(**{xdim: xlims, ydim: ylims}).opts(
        plot={"Overlay": style_overlay, "Path": style_path}
    )
Populated the namespace with:
np, matplotlib, kwant, holoviews, init_notebook, SimpleNamespace, pprint_matrix, scientific_number, pretty_fmt_complex, plt, pf, display_html
from code/edx_components:
MoocVideo, MoocDiscussion, MoocCheckboxesAssessment, MoocMultipleChoiceAssessment, MoocSelfAssessment
from code/functions:
spectrum, hamiltonian_array, h_k, pauli
Using kwant 1.4.2 and holoviews 1.13.2
Executed on 2021-03-06 at 14:55:47.534125.