import sys
sys.path.append("../../code")
from init_mooc_nb import *
init_notebook()
import scipy.sparse.linalg as sl
my = 0.5 * (pauli.sys0 + pauli.sysz)
s0s0sz = np.kron(pauli.s0s0, pauli.sz)
s0szsz = np.kron(pauli.s0sz, pauli.sz)
mys0 = np.kron(my, pauli.s0)
s0s0sx = np.kron(pauli.s0s0, pauli.sx)
s0s0sy = np.kron(pauli.s0s0, pauli.sy)
szsxsz = np.kron(pauli.szsx, pauli.sz)
s0sysz = np.kron(pauli.s0sy, pauli.sz)
sxsxsz = np.kron(pauli.sxsx, pauli.sz)
sysxsz = np.kron(pauli.sysx, pauli.sz)
def make_qshe_sc(l=40, w=10, lead=False):
def shape(pos):
(x, y) = pos
return (1.0 * y ** 2 / l ** 2 + 1.0 * x ** 2 / w ** 2) <= 0.25
def onsite(site, p):
(x, y) = site.pos
return (
(p.M - 4 * p.B) * s0szsz
- 4 * p.D * s0s0sz
+ p.gaps(x, y)[1] * mys0
+ p.gaps(x, y)[0] * s0s0sx
)
def hopx(site1, site2, p):
return p.B * s0szsz + p.D * s0s0sz + 0.5j * p.A * szsxsz
def hopy(site1, site2, p):
return p.B * s0szsz + p.D * s0s0sz - 0.5j * p.A * s0sysz
lat = kwant.lattice.square()
syst = kwant.Builder()
syst[lat.shape(shape, (0, 0))] = onsite
syst[kwant.HoppingKind((1, 0), lat)] = hopx
syst[kwant.HoppingKind((0, 1), lat)] = hopy
if lead:
sym = kwant.TranslationalSymmetry((0, 1))
lead = kwant.Builder(sym)
lead[lat(0, 0)] = 1.5 * p.B * s0szsz
lead[kwant.HoppingKind((0, 1), lat)] = -p.B * s0szsz
syst.attach_lead(lead)
syst.attach_lead(lead.reversed())
return syst.finalized()
def make_qshe_sc_ribbon(w=3):
def ribbon_shape(pos):
(x, y) = pos
return 0 <= x < w
def onsite(site, p):
(x, y) = site.pos
return (
(p.M - 4 * p.B) * s0szsz
- 4 * p.D * s0s0sz
+ p.gaps(x, y)[1] * mys0
+ p.gaps(x, y)[0] * s0s0sx
)
def hopx(site1, site2, p):
return p.B * s0szsz + p.D * s0s0sz + 0.5j * p.A * szsxsz
def hopy(site1, site2, p):
return p.B * s0szsz + p.D * s0s0sz - 0.5j * p.A * s0sysz
lat = kwant.lattice.square()
sym = kwant.TranslationalSymmetry((0, 1))
syst = kwant.Builder(sym)
syst[lat.shape(ribbon_shape, (0, 0))] = onsite
syst[kwant.HoppingKind((1, 0), lat)] = hopx
syst[kwant.HoppingKind((0, 1), lat)] = hopy
return syst.finalized()
def make_2d_pwave(w, l):
def shape(pos):
(x, y) = pos
return (1.0 * y ** 2 / l ** 2 + 1.0 * x ** 2 / w ** 2) <= 0.25
def hopx(site1, site2, p):
(x1, y1) = site1.pos
(x2, y2) = site2.pos
phi = p.phase(0.5 * (x1 + x2), 0.5 * (y1 + y2))
return -p.t * pauli.sz + 1j * p.delta * (
np.cos(phi) * pauli.sx + np.sin(phi) * pauli.sy
)
def hopy(site1, site2, p):
(x1, y1) = site1.pos
(x2, y2) = site2.pos
phi = p.phase(0.5 * (x1 + x2), 0.5 * (y1 + y2))
return -p.t * pauli.sz - 1j * p.delta * (
np.cos(np.pi / 2 + phi) * pauli.sx + np.sin(np.pi / 2 + phi) * pauli.sy
)
def onsite(site1, p):
return (4 * p.t - p.mu) * pauli.sz
lat = kwant.lattice.square()
syst = kwant.Builder()
syst[lat.shape(shape, (w / 2 - 1, 0))] = onsite
syst[kwant.HoppingKind((1, 0), lat)] = hopx
syst[kwant.HoppingKind((0, 1), lat)] = hopy
return syst.finalized()
def bhz_slab(l, w, h):
lat = kwant.lattice.general(np.identity(3))
syst = kwant.Builder()
def shape(pos):
(x, y, z) = pos
return (0 <= z < h) and (1.0 * y ** 2 / l ** 2 + 1.0 * x ** 2 / w ** 2) <= 0.25
def onsite(site, p):
(x, y, z) = site.pos
phi = p.phase(x, y)
return (
(p.C + 2 * p.D1 + 4 * p.D2) * s0s0sz
+ (p.M + 2 * p.B1 + 4 * p.B2) * s0szsz
+ p.delta * (np.cos(phi) * s0s0sx + np.sin(phi) * s0s0sy)
)
def hopx(site1, site2, p):
return -p.D2 * s0s0sz - p.B2 * s0szsz + p.A2 * 0.5j * sxsxsz
def hopy(site1, site2, p):
return -p.D2 * s0s0sz - p.B2 * s0szsz + p.A2 * 0.5j * sysxsz
def hopz(site1, site2, p):
return -p.D1 * s0s0sz - p.B1 * s0szsz + p.A1 * 0.5j * szsxsz
syst[lat.shape(shape, (0, 0, 0))] = onsite
syst[kwant.HoppingKind((1, 0, 0), lat)] = hopx
syst[kwant.HoppingKind((0, 1, 0), lat)] = hopy
syst[kwant.HoppingKind((0, 0, 1), lat)] = hopz
return syst.finalized()
def calc_energies(syst, p, num_orbitals, num_states):
ham = syst.hamiltonian_submatrix(params=dict(p=p), sparse=True).tocsc()
energies, states = sl.eigsh(ham, sigma=0, k=num_states)
densities = (
np.linalg.norm(states.reshape(-1, num_orbitals, num_states), axis=1) ** 2
)
return energies, states, densities