Dear kwant users

I'm trying to calculate magnetic field dependent conductance of a system. For this I used the scattering matrix and adding temperature dependence

by multiplying the transmission value by the derivative of the fermi-dirac function. Because of the results I got, I'm not sure if this is the right way of doing it.

Can someone please help me?

Thank you very much

Best regards

Simon

Code:

fluxes = np.linspace(0, 5, 100)

def make_system(a=1, t=1.5, w=15, d=30):

    lat = kwant.lattice.square(a)

    syst = kwant.Builder()

    def onsite(site, U0, salt):

        return U0 * (uniform(repr(site), repr(salt)) - 0.5) +  4*t

    syst[(lat(x, y) for x in range(w) for y in range(d))] = onsite

    def hopping(site_i, site_j, phi):

        xi, yi = site_i.pos

        xj, yj = site_j.pos

        if yi == yj:

            return -t

        if xi == xj:

            return -0.001 * t * exp(-0.5j * phi * (xi + xj) * (yi - yj))

    syst[lat.neighbors()] = hopping

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

    lead[(lat(j, 0) for j in range(w))] = 4 * t

    lead[lat.neighbors()] = hopping

    syst.attach_lead(lead)

    syst.attach_lead(lead.reversed())

    return syst

def plot_conductance(syst):

    temp = 0.01

    w=15

    d=30

    U0=0

    salt = 0.2

    phis = np.linspace(0, (2 * pi / w) * 5, 100)

    count = 0

    values = []

    for phi in phis:

        count += 1

        ham_mat = syst.hamiltonian_submatrix(params=dict(phi=phi, U0=U0, salt=salt), sparse=True)

        evals = np.real(np.sort(np.linalg.eigvals(ham_mat.todense())))

        EF = evals[int(len(evals) / 2)]

        liste = []

        for k in range(0,len(evals)):

            energy = evals[k]

            if EF-temp <= energy <= EF+temp:

                df = np.real(exp((energy - EF) / temp) / (temp * (exp((energy - EF) / temp) + 1) ** 2))

                Smatrix = kwant.smatrix(syst, energy, params=dict(phi=phi, U0=U0, salt=salt))

                liste.append(Smatrix.transmission(0,1)*df)

        val = sum(liste)/(len(liste))

        values.append(val/(w*d))

        print(count)

    plt.plot(fluxes,values, label = w)

    plt.grid(True)

    plt.title("Conductance eqhop temp=0.1")

    plt.legend()

    plt.ylabel("G_[G_0]")

    plt.xlabel(r'$\phi [2\pi / L]$')

    plt.show()

def main():

    syst = make_system()

    syst = syst.finalized()

    plot_conductance(syst)

if __name__ == "__main__":

    main()