Dear Sir,
I am a PhD student of Hong Kong University of Science and Technology. I
want to use KWANT to caculate Hall resistance of a Hall bar structure.We
can get the conductance between 6 electrodes, but how to get hall
resistance? Can you give me some help? Thank you very much.
Best Regards,
Zhang Bing

Dear all,
I am using Kwant to study a system with both spins and electron and hole,
so the hopping matrix is 4X4. Now I want to know how to separate the spins,
electron and hole. I study this through the "2.6. Superconductors: orbital
degrees of freedom, conservation laws and symmetries"
smatrix = kwant.smatrix(syst, energy)
data.append(smatrix.submatrix((0, 0), (0, 0)).shape[0]
-smatrix.transmission((0, 0), (0, 0)) +smatrix.transmission((0, 1), (0, 0)))
I do not find more information about this for kwant.
For my system, if I have 3 leads, and the hopping and onsite energy is set
in this order:
(e↑,0,0,0)- spin up electron,first row of the 4X4matrix
(0,e↓,0,0)- spin down electron,second row of the 4X4matrix
(0,0, h↑ ,0)- spin up hole,third row of the 4X4matrix
(0,0, 0, h ↓ )- spin down hole,fourth row of the 4X4matrix
I want to obtain the transmissions from 0 →2.
smatrix = kwant.smatrix(syst, energy)
The transmission from h↑ to e↑ is: smatrix.transmission((2, 1), (0, 2))
The transmission from h ↓ to e↑ is: smatrix.transmission((2, 1), (0, 3))
Is my understanding correct?
Thanks very much in advance!
Hosein Khani

Hello Kwant community,
I write since I am trying to simulate a Kitaev chain with disorder,
here distributed uniformly in the range [1,-1],
and I am experiencing troubles to implement the disorder itself.
In the absence of disorder,
the function that defines the tight-binding (open) Hamiltonian can be written
e.g. as follows:
def kitaev(L ,t, mu, Delta) :
syst = kwant.Builder()
lat = kwant.lattice.chain(norbs=2)
syst[lat(0)] = - (mu/2)*pauli3
for i in range(1,L):
syst[lat(i)] = -(mu/2)*pauli3
syst[lat(i), lat(i - 1)] = -(t/2)*pauli3 + 1j*(Delta/2)*pauli2
return syst
where the Pauli matrices have been previously defined, t is the hopping,
mu is the chemical potential, and Delta is the gap.
Now, I want to add disorder, proportional to W and uniformly
distributed in [-1,1]. The tutorial online for Kwant states that the
function kwant.digest.uniform(input, salt='') creates a random offset
in the interval [0,1]. Therefore I thought that the Hamiltonian could
be changed e.g. as follows:
def kitaev(L ,t, mu, Delta, W) :
syst = kwant.Builder()
lat = kwant.lattice.chain(norbs=2)
syst[lat(0)] = -(1/2)*(mu + W*(2*onsite - 1))*pauli3
for i in range(1,L):
syst[lat(i)] = -(1/2)*(mu + W*(2*onsite - 1))*pauli3
syst[lat(i), lat(i - 1)] = -(t/2)*pauli3 + 1j*(Delta/2)*pauli2
return syst
after that the onsite function has been defined as
def onsite(site, phi, salt):
return uniform(repr(site), salt)
(I am still trying to follow the tutorial).
However, when trying to calculate eigenvalues and eigenvectors,
via the code
L = 50
W = 0
mub = 3
systb = kitaev(L ,1, mub, 1, W)
systb = systb.finalized()
hamatb = systb.hamiltonian_submatrix(sparse = False)
evalsb, evecsb = la.eigh(hamatb)
I encountered the problem
-----------------------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-10-12a574eefcd0> in <module>
3
4 mub = 3
----> 5 systb = kitaev(L ,1, mub, 1, W)
6 systb = systb.finalized()
7 hamatb = systb.hamiltonian_submatrix(sparse = False)
<ipython-input-9-661d5ddb539a> in kitaev(L, t, mu, Delta, W)
4 lat = kwant.lattice.chain(norbs=2)
5
----> 6 syst[lat(0)] = -(1/2)*(mu + W*(2*onsite - 1))*pauli3
7
8 for i in range(1,L):
TypeError: unsupported operand type(s) for *: 'int' and 'function'
—————————————————————————————------
that instead does no occur in the absence of disorder (then when the
first definition above for the Kitaev chain is assumed).
Probably I did not understand properly how to use the
function uniform(repr(site), salt).
Can you help me ?
Thank you very much in advance and best regards
L. L.

Dear Rudolf,
indeed, manybody expectation values involving thermal averages
over the occupied states are only available in tkwant by now.
The "voltag_raise" tutorial you mentioned is a good reference, just make sure
that you stick with the new version (tkwant's API has slightly changed):
https://kwant-project.org/extensions/tkwant/tutorial/AC_josephson_without_s…
For spin currents, it should work similar to the kwant operator tutorial you mention, with
an operator kwant.operator.Current(syst, sigma_u)
and sigma_u = ((1,0),(0,0)) for the upp current,
respectively sigma_d = ((0,0),(0,1)) for the down current.
Best,
Thomas

Dear Kwant users,
Solving Schrodinger equantions describing quantum physics of noninteracting systems amounts to solving eigenproblems. As I know, the Kwant package has the superiority at calculating transport properties, and moreover, it can also solve eigenproblems of closed systems. To be precise, firstly, one can use "kwant.continuum.discretize" to discrete continous Hamiltonians into tight-binding models with a specified lattice constant, and secondly, diagonalize the lattice model to obtain the eigen-energies and eigen-wavefunctions. Actually, we can solve Schrodinger equantions, which are essentially second-order partial differential equations, by numerical methods such as the powerful finite-element-method with irregular mesh and specified boudary conditions. I know that many papers in codensed matter community studying finite-size models by the tight-binding method rather than the finite-element-method. Is this due to any restriction of the finite-element-method in solving Schrodinger equantion? Does anyone ever compare the performance of tight-binding-method and finite-element-method on solving eigenproblem of closed Hamiltonian. Is there any literature addressing this issue?
Regards,
Zhan

Dear kwant users
I was trying to implement the following Hamiltonian: H_j,j+1 = H_j+1,j = -t, H_j,j = -0.02*t*cos(k_y + 2pi (phi/W)*j)
for a square lattice of size W*L with attached leads in y-direction. Unfortunately it doesn't work. Could anyone give me a hint
how to do this the right way?
Thank you for your answers in advance
Best,
Simon
def make_system(a=1, L=60, W=30):
offdiagmat = np.zeros((W, W))
for i in range(W):
for j in range(W):
if i == j + 1 or j == i + 1:
offdiagmat[i][j] = 1
hamiltonian = "-1.0*offdiagmat - identity(4)*2*0.01*cos(k_y + (2*pi*phi/30)*x)"
hamiltonian = kwant.continuum.sympify(hamiltonian, locals=dict(offdiagmat=offdiagmat))
template = kwant.continuum.discretize(hamiltonian, grid=a)
def shape(site):
(x, y) = site.pos
return (0 <= x < W and 0 <= y < L)
def lead_shape(site):
(x, y) = site.pos
return (0 <= x < W)
syst = kwant.Builder()
syst.fill(template, shape, (0, 0))
lead = kwant.Builder(kwant.TranslationalSymmetry([0, -a]))
lead.fill(template, lead_shape, (0, 0))
syst.attach_lead(lead)
syst.attach_lead(lead.reversed())
syst = syst.finalized()
return syst
def analyze_system():
params = dict(phi=0.2)
syst = make_system()
kwant.plotter.bands(syst.leads[0], params=params, momenta=np.linspace(-pi, pi, 200), show=False)
plt.grid(True)
plt.xlim(-pi, pi)
plt.xlabel('momentum [1/A]')
plt.ylabel('energy [eV]')
plt.show()
def main():
analyze_system()
if __name__ == '__main__':
main()

Dear Sir,
I was trying to calculate the Integer quantum Hall effect in graphene with the help of Laughlinargument.ipynb. For square lattice the quantization is okay ( ne^2/h). But for graphene it is (2n+1) 2e^2/h. I am getting the quantization like 1, 3, 5, 7 e^2/h in graphene which seems there is missing of factor 2. Is the Hall conductivity I am getting is in 2e^2/h unit or should I multiply the data by 2 to get the quantization like 2, 6, 10 which is well-known IQHE in graphene.
It will be really helpful If you kindly clear my doubts regarding this.
Thanking You,
Priyanka Sinha

Dear kwant users
I'm trying to implement the interband-kubo formula for dc-conductivity.
The underlying system is a square lattice with attached leads in z-direction and perpendicular magnetic field.
My problem is that I couldn't manage to calculate the expectation value for the current operator for different bands.
Could anyone give me a hint what is the right way of doing it?
Best regards
Simon
def fermi(energy, Efermi):
f = np.real(1 / (exp((energy - Efermi) / temp) + 1))
return f
def make_system(a=1, t=1.0, w=15, d=60):
lat = kwant.lattice.square(a, norbs=1)
syst = kwant.Builder()
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 and yi != yj:
return (-0.01 * t) * exp(-0.5j * phi * (xi + xj))
syst[(lat(x, y) for x in range(w) for y in range(d))] = 4 * t
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
syst = make_system()
syst = syst.finalized()
def Interband_Kubo(syst):
w=15
d=60
Gs_inter = []
count = 0
phis = np.linspace(0, (2 * pi / 15) * 5, 100)
for phi in phis:
count += 1
bands = kwant.physics.Bands(syst.leads[1], params=dict(phi=phi))
momenta = np.linspace(-pi, pi, 300)
energydisp = [bands(k) for k in momenta]
elist = []
for m in range(0, len(energydisp), 1):
for n in range(0, len(energydisp[m]), 1):
energy = np.real(energydisp[m][n])
elist.append(energy)
selist = np.sort(elist)
EF = selist[int(len(elist) / 2)]
condsum = []
for l in range(0, len(energydisp), 1):
for s in range(0, len(energydisp[l]), 1):
count = 0
for q in range(0, len(energydisp[l]), 1):
count += 1
if s != q:
#print(count)
energy_s = np.real(energydisp[l][s])
energy_q = np.real(energydisp[l][q])
if EF -40*temp <= energy_s <= EF + 40*temp:
if EF - 40 * temp <= energy_q <= EF + 40* temp:
J_z = kwant.operator.Current(syst)
wf_s = kwant.solvers.default.wave_function(syst, energydisp[l][s], params=dict(phi=phi))
wf_q = kwant.solvers.default.wave_function(syst, energydisp[l][q], params=dict(phi=phi))
if len(wf_s(1)) == 1 and len(wf_q(1)) ==1:
psi_s = wf_s(1)[0]
psi_q = wf_q(1)[0]
print(psi_s)
J_qs = np.sum(J_z(psi_q, psi_s, params=dict(phi=phi)))
J_sq = np.sum(J_z(psi_s, psi_q, params=dict(phi=phi)))
else:
J_qs = 0
J_sq = 0
cond = 1j * ((fermi(energy_s, EF) - fermi(energy_q, EF))) * (
(J_sq * J_qs) / (energy_s - energy_q + (1 / 15 ** 2) * 1j))
condsum.append(cond)
conductance = (np.sum(condsum) / (15*60))
print('conductance', conductance)
Gs_inter.append(conductance)
return Gs_inter
fluxes = np.linspace(0, 5, 100)
plt.plot(fluxes, Interband_Kubo(syst), label = "inter")
plt.grid(True)
plt.title(" Conductance")
plt.legend()
plt.ylabel("conductance")
plt.xlabel(r'$\phi [2\pi / L]$')
plt.show()

Dear (t)kwant community,
applying an AC voltage in a time-dependent setting involves
tkwant.leads.add_voltage(sys,phase) as proposed in
https://gitlab.kwant-project.org/jbweston/tkwant-tutorial/blob/master/1_sys…,
where phase is the voltage phase, with which the hopping from a newly
created (set of) system-site(s) is transformed in a Peierls-substitution
like manner. The documentation says, that "phase" should be a callable
returning either a scalar, or a sequence of values, where the number of
elements in this sequence is corresponding to the number of orbitals on
a site.
For my purpose, I am working with norbs = 2, trying to simulate a system
of spinfull fermions on a lattice, of say, after setting the
quantization axis, spin-up and spin-down electrons. The phase I defined
like this by a minor change of the original tutorial:
def faraday_flux(time, V_dynamic):
omega = 0.5
t_upper = pi / omega
if time <= 0:
return [0,0]
if 0 < time < t_upper:
return [V_dynamic * (time - sin(omega * time) / omega)/2,V_dynamic *
(time - sin(omega * time) / omega)/2]
return [V_dynamic * (time - t_upper) + V_dynamic * t_upper/2,V_dynamic *
(time - t_upper) + V_dynamic * t_upper/2]
Then, after finishing my system and adding leads to it, I call
tkwant.leads.add_voltage(syst, 0, faraday_flux)
in order to apply this dynamic phase "faraday_flux" to lead 0.
*Is* this the proper way to define the same ac voltage phase for both
orbitals*?* If I call a single scalar only (instead of a list of
scalars), I get the same current-time characteristics. I would actually
expect calling a single scalar for a system with 2 orbitals per site
should not work in first place, and if it is handled cleverly, it should
definitively not yield the same result (since one part of the electron
hilbertspace has a vanishing probability of hopping onto the system and
thus a vanishing current, which in turn should be present in the total
current).
I observe the current with the following sequence:
1. generate system: make_syst(lat,params,...)
2. add ac voltage: tkwant.leads.add_voltage(syst,0,faraday_flux)
3. finalize system
4. define operator kwant.operator.Current(fsyst, spinc,
where=interface,sum=True)
5. define solver with tkwant.manybody.State(fsyst, tmax, occupation,
params=params)
6. propagate solver in time with solver.evolve(time) and evaluate the
current with the solver at a specific time solver.evaluate(operator)
I am using tkwant.manybody.State in order to incorporate a specific
chemical potential and temperature in the lead, which, to my knowledge,
is not possible with the kwant.wave_function (correct me if I am
wrong!). Another reason I would not use kwant.wave_function is because
it is just computing wave functions for specific energies, I am however
interested in quantities that are already integrated over the energy. As
a reference, I am using the voltage_raise tkwant tutorial, which can be
found here:
https://gitlab.kwant-project.org/jbweston/tkwant/blob/k-integration/tutoria…
.
One other thing which I have not been able to figure out yet: How to
access the individual spin components for a 2-orbital system of the
manybody-wavefunction? This is possible with kwant.wave_function as
shown in https://kwant-project.org/doc/dev/tutorial/operators , however
when it comes to the time-dependent equivalent, I am a bit clueless. My
aim is to observe a spin-resolved current, where I only address the
individual components of my observables originating from a single
contribution of spin-up (-down) electrons only. Subsequently I started
to ask myself, how to measure a spin-current? Is defining the operator
as above with setting spinc = sigma_z doing the job?
Thank you for your answers in advance!
Best,
Rudolf

Hi All,
Short question regarding the definition of the hamiltonian in an SNS
junction. I have the following hamiltonian:
### S
syst[(lat(x, 0) for x in range(-S, 0))] = (2 * t - mu) * tau_z + Delta *
tau_x
#### N
syst[(lat(x, 0) for x in range(0, N))] = (2 * t - mu + barrier) * tau_z
### S
syst[(lat(x, 0) for x in range(N, N+S))] = (2 * t - mu) * tau_z + Delta *
tau_x
with N=10 and S=20. However the wavefunction I get in return with the
boundstate algorithm (as this one, for example,
https://gitlab.kwant-project.org/kwant/boundstate) has twice the number of
sites in each region (N=20 and S =40).
Does anyone have any idea why it has twice the number of points and if I
can trust the solution in these points?
Here is an example that I've sent in the previous email.
[image: image.png]
Best,
Denise