I am trying to plot the wavefunction for a 2D NISIN Hamiltonian however I am unable to get a result, this is my code:

 

 

import kwant

import tinyarray as ta

import numpy as np

import matplotlib.pyplot as plt

from tqdm import tqdm

import winsound

import os

from datetime import datetime

import scipy.sparse.linalg as sla

 

 

sig_0 = np.identity(4)

s_0 = np.identity(2)

s_z = np.array([[1., 0.], [0., -1.]])

s_x = np.array([[0., 1.], [1., 0.]])

s_y = np.array([[0., -1j], [1j, 0.]])

 

tau_z = ta.array(np.kron(s_z, s_0))

tau_x = ta.array(np.kron(s_x, s_0))

tau_y = ta.array(np.kron(s_y, s_0))

sig_z = ta.array(np.kron(s_0, s_z))

tau_zsig_x = ta.array(np.kron(s_z, s_x))

tau_zsig_y = ta.array(np.kron(s_z, s_y))

 

# Lorentzian definition

 

 

def lorentzianxy(x, y, ind,y0):

        gam = 3

        L = params["L"]

        #if 0<ind<L: 

        return gam ** 2 /( gam ** 2 + (x - ind) ** 2 + (y - y0) ** 2)

       # else :

        #    return 0

 

# Define potential

       

def potential(site, pot,ind,y0):

        (x, y) = site.pos

        return pot * lorentzianxy(x,y,ind,y0)

 

 

# Make system

 

def makeNISIN2D(params):

 

# List all the params

    W=params["W"]

    L=params["L"]

    pot = params["pot"]

    ind = params["ind"]

    mu=params["mu"]

    B=params["B"]

    Delta=params["Delta"]

    alpha=params["alpha"]

    t=params["t"]

    barrier = params["barrier"]

    Smu = params["Smu"]

   

   

 

    def lorentzianxy(x, y, ind,y0):

        gam = 3

        L = params["L"]

        #if 0<ind<L: 

        return gam ** 2 /( gam ** 2 + (x - ind) ** 2 + (y - y0) ** 2)

       # else :

        #    return 0

 

# Define potential

       

    def potential(site, pot,ind,y0):

        (x, y) = site.pos

        return pot * lorentzianxy(x,y,ind,y0)

#        if y < lorentzian(x,ind)*4:

#               

#                return (tru_pot)

#        else:

#                return 0

#           

    def Straightpotential(site, pot,ind):

        (x, y) = site.pos

        if x == ind:

            return pot

        else:

            return 0

       

    def onsiteSc(site, pot,ind,y0):

        #(x,y)=site.pos

        #L=params["L"]

       # if x>L/2:

       return (4 * t - Smu) * tau_z + B * sig_z + Delta * tau_x - potential(site, pot,ind,y0)*tau_z

       # else:

       #return (4 * t ) * tau_z + B * sig_z + Delta * tau_x - potential(site, pot,ind,y0)*tau_z

       

    def onsiteNormal(site, pot,ind,y0):

        return (4 * t - mu) * tau_z #- potential(site, pot,ind,y0)*tau_z

    def onsiteBarrier(site, pot,ind,y0):

        return (4 * t - mu + barrier) * tau_z #- potential(site, pot,ind,y0)*tau_z

  

    def hop_x(site, pot, ind,y0):

       

        return -t * tau_z + 0.5j * alpha * tau_zsig_y

   

    def hop_y(site, pot, ind,y0):

   

           return -t * tau_z - .5j * alpha * tau_zsig_x

   

 

   

    # On each site, electron and hole orbitals.

    lat = kwant.lattice.square(norbs=4)

    syst = kwant.Builder()

   

    # S

    syst[(lat(i, j) for i in range(1,L-1) for j in range(W))] = onsiteSc

   

    barrierLen = 1;

    # I's

    syst[(lat(i, j) for i in range(barrierLen) for j in range(W))] = onsiteBarrier

    syst[(lat(i, j) for i in range(L-barrierLen, L)for j in range(W))] = onsiteBarrier

    syst[kwant.builder.HoppingKind((1, 0), lat, lat)]= hop_x

    syst[kwant.builder.HoppingKind((0, 1), lat, lat)]= hop_y

   

    # N's

    lead = kwant.Builder(kwant.TranslationalSymmetry((-1,0)),

                         conservation_law=-tau_z#, particle_hole=tau_y

                         )

    lead[(lat(0, j) for j in range(W))] = onsiteNormal

    lead[kwant.builder.HoppingKind((1, 0), lat, lat)]= hop_x

    lead[kwant.builder.HoppingKind((0, 1), lat, lat)]= hop_y

   

    syst.attach_lead(lead)

    syst.attach_lead(lead.reversed())

 

    return syst

 

 

def sorted_eigs(ev):

    evals, evecs = ev

    evals, evecs = map(np.array, zip(*sorted(zip(evals, evecs.transpose()))))

    return evals, evecs.transpose()

 

params = dict(mu=0.3, Delta=.1, alpha=.8, t=1.0, barrier = 2.0, pot = 0.0,W = 5, L = 30, ind = 0, B = 0.3, Smu=0.00,y0=0)

 

sys = makeNISIN2D(params)

sys = sys.finalized()

 

    # Calculate the wave functions in the system.

ham_mat = sys.hamiltonian_submatrix(sparse=True, params = params)

 

evals, evecs = sorted_eigs(sla.eigsh(ham_mat.tocsc(), k=20, sigma=0))

 

# Plot the probability density of the 10th eigenmode.

kwant.plotter.map(sys, np.abs(evecs[:, 9])**2,

                  colorbar=False, oversampling=1)

 

 

 

 

And the error I get is this: “ValueError: The number of sites doesn't match the number of provided values.”

 

When I print the values they indeed are more than what they should be. I’ve tried a couple of things but I’m quite at a loss of what to do. I want to plot the wavefunction so I can see the Majoranas at the edges of the superconductor so that I can see how they respond to a local perturbation in the potential. Any tips would be very helpful.

 

 

 

Kind regards,

 

Hanifa