In [1]:
import sys

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

init_notebook()
%output size = 150
import scipy
from matplotlib import cm

bhz_parameters = {
"topo": {"A": 0.5, "B": 1.00, "D": 0.0, "M": 1.0, "del_z": 0.0},
"triv": {"A": 0.5, "B": 1.00, "D": 0.0, "M": -1.0, "del_z": 0.0},
"topo2": {"A": 0.5, "B": 1.00, "D": 0.3, "M": 1.0, "del_z": 0.0},
"slowed": {"A": 0.05, "B": 0.08, "D": 0.15, "M": -0.3, "del_z": 0.5},
}

# Onsite and hoppings for bhz model

def onsite(site, p):
return (p.M - 4 * p.B) * pauli.s0sz - 4 * p.D * pauli.s0s0 + p.del_z * pauli.sysy

def hopx(site1, site2, p):
return p.B * pauli.s0sz + p.D * pauli.s0s0 + 1j * p.A * pauli.szsx

def hopy(site1, site2, p):
return p.B * pauli.s0sz + p.D * pauli.s0s0 - 1j * p.A * pauli.s0sy

def bhz(w=20):
""" Make ribbon system with bhz model.

slowed parameters are used on the edge for finite size system.
"""
lat = kwant.lattice.square()

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

slowed_par = SimpleNamespace(Bz=0, **bhz_parameters["slowed"])
if w is None:
syst = kwant.Builder(kwant.TranslationalSymmetry(*lat.prim_vecs))
syst[lat.shape(lambda pos: True, (0, 0))] = onsite
syst[kwant.HoppingKind((1, 0), lat)] = hopping_x
syst[kwant.HoppingKind((0, 1), lat)] = hopy
else:
syst = kwant.Builder(kwant.TranslationalSymmetry((1, 0)))
syst[(lat(0, i) for i in range(w))] = onsite
syst[kwant.HoppingKind((1, 0), lat)] = hopping_x
syst[kwant.HoppingKind((0, 1), lat)] = hopy

syst[lat(0, -1)] = lambda site, p: onsite(site, slowed_par)
syst[lat(1, -1), lat(0, -1)] = lambda site1, site2, p: hopping_x(
site1, site2, slowed_par
)
syst[lat(0, 0), lat(0, -1)] = hopy

syst[lat(0, w)] = lambda site, p: onsite(site, slowed_par)
syst[lat(1, w), lat(0, w)] = lambda site1, site2, p: hopping_x(
site1, site2, slowed_par
)
syst[lat(0, w), lat(0, w - 1)] = hopy

return syst

def bhz_cylinder(w=3):
""" Make cylinder system with bhz model. """

def ribbon_shape(pos):
(x, y) = pos
return 0 <= y < w

lat = kwant.lattice.square(norbs=4)
sym = kwant.TranslationalSymmetry((1, 0))
syst = kwant.Builder(sym)

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

def hopy_phase(site1, site2, p):
x1, y1 = site1.pos
x2, y2 = site2.pos
return hopy(site1, site2, p) * np.exp(1j * p.ky)

syst[lat.shape(ribbon_shape, (0, 0))] = onsite
syst[kwant.HoppingKind((1, 0), lat)] = hopping_x
syst[kwant.HoppingKind((0, 1), lat)] = hopy
syst[kwant.HoppingKind((0, -w + 1), lat)] = hopy_phase

syst[lat(0, 0)] = onsite
syst[lat(0, w - 1)] = onsite

return syst

def ribbon_shape(pos):
(x, y) = pos
return 0 <= y < 2

lat = kwant.lattice.square(norbs=4)
sym = kwant.TranslationalSymmetry((1, 0))
syst = kwant.Builder(sym, time_reversal=1j * pauli.sys0)

syst[lat.shape(ribbon_shape, (0, 0))] = 1.8 * t * pauli.s0sz
syst[kwant.HoppingKind((1, 0), lat)] = -t * pauli.s0sz
return syst

def make_scatter_sys():
def shape(pos):
x, y = pos
return (x == 0) * (0 <= y < 3)

lat = kwant.lattice.square(norbs=4)
syst = kwant.Builder()
syst[lat.shape(shape, (0, 0))] = onsite
syst[kwant.HoppingKind((0, 1), lat)] = hopy

syst = syst.finalized()
return syst

def scattering_det_pfaff(syst, p):
def pfaffian(syst, p, ky):
p.ky = ky
smat = kwant.smatrix(syst, energy=0.0, params=dict(p=p)).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)

pfaff = [pfaffian(syst, p, 0), pfaffian(syst, p, np.pi)]

ks = np.linspace(0.0, np.pi, 50)
det = [np.linalg.det(kwant.smatrix(syst, energy=0.0, params=dict(p=p)).data) for p.ky in ks]
det = np.array(det)

phase = np.angle(pfaff[0]) + 0.5 * np.cumsum(np.angle(det[1:] / det[:-1]))
kdims = ["$k_y$", "phase"]
plot = holoviews.Path((ks[1:], phase), kdims=kdims).opts(style={"color": "b"})
plot *= holoviews.Points(([0, np.pi], np.angle(pfaff)), kdims=kdims).opts(
style={"color": "g"}
)
xlims, ylims = slice(-0.2, np.pi + 0.2), slice(-np.pi - 0.2, np.pi + 0.2)
pi_ticks = [(-np.pi, r"$-\pi$"), (0, "$0$"), (np.pi, r"$\pi$")]
ticks = {"xticks": [(0, "0"), (np.pi, "$\pi$")], "yticks": pi_ticks}
return plot.relabel("Winding", depth=1)[xlims, ylims].opts(plot=ticks)

def title(p):
title = r"$A={:.2}$, $B={:.2}$, $D={:.2}$, $M={:.2}$"
return title.format(p.A, p.B, p.D, p.M)

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:59.467032.