Conductivity in 2D system with Zeeman and Rashba spin-orbit interaction
Dear Kwant community,
I would like to calculate the DOS and conductivity in 2D system with Zeeman and Rashba spin-orbit interaction in a square lattice. Regarding the chapter 2.3.1 in KWANT tutorial I defined the function make_system(). Then, according to the chapter 2.9 in the tutorial, I used the kwant.kpm.conductivity() method to calculate the longitudinal and transverse conductivity in the system. But I see that my results are wrong and I do not understand some things:
1. question: in chapter 2.3.1 we did not need to define the norbs parameter but I suppose that KWANT "knows" that we have here two orbitals per site (considering the spin) because we use the Pauli matrices. But, to obtain some results using function kwant.kpm.conductivity(), I see that the argument norbs=2 in the system deifinition is required (without it I get the error: 'NoneType' object cannot be interpreted as an integer). Thus, my quesion is: how does KWANT interpret my system if I add or not this argument (nobrs=2)? I think this could be the problem causing my wrong results.
2. question: To use the method kwant.kpm.conductivity() we don't need to attach the leads. Maybe it is obvious and trivial quesion, but I do not understand how we can calculate the conductivity (eg. induced by the external electric field) in a closed system? And what about the boundary conditions?
Below I present the whole code.
#-------------------------------------------------------------------------------------- #-------------------------------------------------------------------------------------- import kwant from matplotlib import pyplot import tinyarray
# we define the identity matrix and Pauli matrices sigma_0 = tinyarray.array([[1,0], [0,1]]) sigma_x = tinyarray.array([[0,1], [1,0]]) sigma_y = tinyarray.array([[0,-1j], [1j,0]]) sigma_z = tinyarray.array([[1,0], [0,-1]])
def make_system(t=1.8, alpha=0.1, e_z=0.0018, W=30, L=30, a=1): lat = kwant.lattice.square(norbs=2) # norbs - the number of orbitals per site syst = kwant.Builder()
syst[(lat(x,y) for x in range(L) for y in range(W))] = 4*t*sigma_0 + e_z*sigma_z # Zeeman energy adds to the on-site term syst[kwant.builder.HoppingKind((a,0), lat, lat)] = -t*sigma_0 - 1j*alpha*sigma_y/(2*a) # hoppings in x-direction; (a,0) means hopping form (i,j) to (i+1,j) syst[kwant.builder.HoppingKind((0,a), lat, lat)] = -t*sigma_0 + 1j*alpha*sigma_x/(2*a) # hoppings in y-direction; (0,a) means hopping form (i,j) to (i,j+1)
return syst
def plot_dos_and_curves(dos, labels_to_data): pyplot.figure() pyplot.fill_between(dos[0], dos[1], label="DoS [a.u.]", alpha=0.5, color='gray') for label, (x, y) in labels_to_data: pyplot.plot(x, y, label=label, linewidth=2) pyplot.legend(loc=2, framealpha=0.5) pyplot.xlabel("energy [t]") pyplot.ylabel("$\sigma [e^2/h]$") pyplot.show()
def main(): syst = make_system() kwant.plot(syst); # check if the system looks as intended syst = syst.finalized()
spectrum = kwant.kpm.SpectralDensity(syst, rng=0, num_vectors=30) # object that represents the density of states for the system energies, densities = spectrum()
# component 'xx' cond_xx = kwant.kpm.conductivity(syst, alpha='x', beta='x', rng=0) # component 'xy' cond_xy = kwant.kpm.conductivity(syst, alpha='x', beta='y', rng=0)
energies = cond_xx.energies cond_array_xx = np.array([cond_xx(e, temperature=0.01) for e in energies]) cond_array_xy = np.array([cond_xy(e, temperature=0.01) for e in energies])
plot_dos_and_curves( (spectrum.energies, spectrum.densities * 8), # why multiplied by 8? [ (r'Longitudinal conductivity $\sigma_{xx} / 4$', (energies, cond_array_xx.real / 4)), (r'Hall conductivity $\sigma_{xy}$', (energies, cond_array_xy.real))], )
# standard Python construct to execute 'main' if the program is used as a script if __name__ == '__main__': main() #-------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------
Regards, Anna Krzyzewska
Dear Anna,
1. In the newer versions of Kwant you must define the `norbs`, it may not be up to date everywhere in the documentation. Certainly the `kpm` module makes use of this parameter of your system.
2. The Kubo conductivity is calculated as a trace per unit cell, you can see the details in the references cited in the documentation of the `kpm.conductivity` function. For a finite system with open boudary conditions, the results will make sense if you compute this local quantity in the bulk of your system far away from the edges. To do this, you can make use of the `vector_factory` argument in the `kpm` module, for example with the `kwant.kom.LocalVectors` class.
The finite system does not need to have leads as we don't compute currents or scattering states, we compute the linear response (in some direction `x_alpha`) of the system to an electric field in some other direction (`x_beta`). These directions enter in the conductivity function with the `alpha`, and `beta` arguments.
For systems with periodic boundary conditions you need to define periodic velocity operatos and pass them to `kpm.conductivity`. This is not implemented in Kwant and can be tricky to define.
In both cases, you need to normalize the results by the area that each of your vectors cover.
By the way, I tried running your sistem and does not show a gap. Therefore, the conductivities are not quantized to a topological invariant.
Here a snippet of how you could define the local vectors in your case
``` edge_width = min(L, W) / 4
def center(s): return (edge_width <= s.pos[0] < L - edge_width) and (edge_width <= s.pos[1] < W - edge_width)
center_vectors = np.array(list(kwant.kpm.RandomVectors(syst, where=center))) # here kwant.kpm.LocalVectors can also be used, but will produce much more vectors num_vectors = len(center_vectors) norm = np.linalg.norm(center_vectors[0]) ** 2
num_vectors=num_vectors, vector_factory=center_vectors ```
I hope this helps, regards, Pablo
Dear Anna,
I modified your Hamiltonian and exaggerated the value of the parameters in order to have a larger gap.
If you want to use the values that you had the first time, your gap will be much smaller. In that case, you need more energy resolution in the KPM expansion, and that will slow down things a lot (just beware). Likely you will need a larger sample when the gap is small.
``` import kwant from matplotlib import pyplot import tinyarray import numpy as np # we define the identity matrix and Pauli matrices
sigma_0 = tinyarray.array([[1,0], [0,1]]) sigma_x = tinyarray.array([[0,1], [1,0]]) sigma_y = tinyarray.array([[0,-1j], [1j,0]]) sigma_z = tinyarray.array([[1,0], [0,-1]])
def make_system(t=1.8, alpha=0.1, e_z=0.0018, W=30, L=30, a=1): lat = kwant.lattice.square(norbs=2) # norbs - the number of orbitals per site syst = kwant.Builder()
syst[(lat(x,y) for x in range(L) for y in range(W))] = 4*t*sigma_0 + e_z*sigma_z # Zeeman energy adds to the on-site term syst[kwant.builder.HoppingKind((a,0), lat, lat)] = -t*sigma_z+ 1j*alpha*sigma_y/(2*a) # # hoppings in x-direction; (a,0) means hopping form (i,j) to (i+1,j) syst[kwant.builder.HoppingKind((0,a), lat, lat)] = -t*sigma_z + 1j*alpha*sigma_x/(2*a) # # hoppings in y-direction; (0,a) means hopping form (i,j) to (i,j+1)
return syst
def plot_dos_and_curves(dos, labels_to_data): pyplot.figure() pyplot.fill_between(dos[0], dos[1], label="DoS [a.u.]", alpha=0.5, color='gray') for label, (x, y) in labels_to_data: pyplot.plot(x, y, label=label, linewidth=2) pyplot.legend(loc=2, framealpha=0.5) pyplot.ylim(-1, 3) pyplot.xlabel("energy [t]") pyplot.ylabel("$\sigma [e^2/h]$") pyplot.show()
def main(random_vecs=True): num_moments = 400
W=30 L=30 syst = make_system(W=W, L=L, alpha=1, e_z=0.5) syst = syst.finalized()
fam = syst.sites[0].family # assuming all sites from the same (sub)lattice area_per_orb = np.cross(*fam.prim_vecs) / fam.norbs
if random_vecs: edge_width = min(L, W) / 4
def center(s): return (edge_width <= s.pos[0] < L - edge_width) and (edge_width <= s.pos[1] < W - edge_width)
center_vectors = kwant.kpm.RandomVectors(syst, where=center) one_vector = next(center_vectors) num_vectors = 10 # len(center_vectors) else: # use local vectors edge_width = min(L, W) / 4
def center(s): return (edge_width <= s.pos[0] < L - edge_width) and (edge_width <= s.pos[1] < W - edge_width) center_vectors = np.array(list(kwant.kpm.LocalVectors(syst, where=center))) one_vector = center_vectors[0] num_vectors = len(center_vectors)
norm = area_per_orb * np.linalg.norm(one_vector) ** 2
# check the area that the vectors cover kwant.plot(syst, site_color=np.abs(one_vector[::2]) + np.abs(one_vector[1::2]), cmap='Reds');
spectrum = kwant.kpm.SpectralDensity(syst, rng=0,num_moments=num_moments, num_vectors=num_vectors, vector_factory=center_vectors) # object # that represents the density of states for the system energies, densities = spectrum()
# component 'xx' cond_xx = kwant.kpm.conductivity(syst, alpha='x', beta='x', num_moments=num_moments, rng=0, num_vectors=num_vectors, vector_factory=center_vectors) # component 'xy' cond_xy = kwant.kpm.conductivity(syst, alpha='x', beta='y', num_moments=num_moments, rng=0, num_vectors=num_vectors,vector_factory=center_vectors)
energies = cond_xx.energies cond_array_xx = np.array([cond_xx(e, temperature=0.01) for e in energies]) / norm cond_array_xy = np.array([cond_xy(e, temperature=0.01) for e in energies]) / norm
scale_to_plot = 1 / (edge_width) plot_dos_and_curves( (spectrum.energies, spectrum.densities / norm), [(r'Longitudinal conductivity $\sigma_{xx}$ (scaled)', (energies, cond_array_xx.real * scale_to_plot)), (r'Hall conductivity $\sigma_{xy}$', (energies, cond_array_xy.real))], ) ```
Dear Pablo,
Thank You very much for Your really helpful replies!
I looked into Your advices and solution, I understood the code, and now some new issues appear.
Thanks to You, I can see now a beautiful QHE in the energy gap, but I am worry about the large noises in transverse conductivity beyond the energy gap (I suppose in this region it should be zero). The second thing is that the energy gap is not in the zero-energy, but around 4t. Thus, I want to check how the band structure looks like for this system with given parametres. According to the KWANT documentation, I found a method to calculate the band structure: kwant.physics.Bands() but it works only for infinite system.
Therefore, I've created two systems: 1. with leads - to calculate a band structure 2. closed (the previous one) - to calculate the components of conductvity tensor with local vectors defined away from the edges (as You adviced).
But here I encounter a problem with the first one: despite I know that the system do has attached leads (so it is infinite), I get an error that the method kwant.plotter.bands() is 'Expecting an instance of InfiniteSystem.'
The second question is (assuming that this solution would works): are the results obtained for this two systems can be related with each other? Means, can I assume that the band structure plotted for the first system will corresponds to the second system, where I obtained the result for conductivites?
I think also about a different approach: to create a system with four leads and to calculate (longitudinal and transverse) conductivities as a transmissions throught a proper leads using the scattering matrix. Maybe that solution could be more convenient for this problem.
Best regards, Anna
Below I attach the whole code.
#------------------------------------------------------------------ import kwant from matplotlib import pyplot import tinyarray import numpy as np
# we define the identity matrix and Pauli matrices sigma_0 = tinyarray.array([[1,0], [0,1]]) sigma_x = tinyarray.array([[0,1], [1,0]]) sigma_y = tinyarray.array([[0,-1j], [1j,0]]) sigma_z = tinyarray.array([[1,0], [0,-1]])
def make_system(t=1, alpha=1, e_z=0.5, W=30, L=30, a=1): # this parameters can be overwritten in the 'main()' section lat = kwant.lattice.square(norbs=2) # norbs - the number of orbitals per site syst = kwant.Builder()
# Zeeman energy adds to the on-site term syst[(lat(x,y) for x in range(L) for y in range(W))] = 4*t*sigma_0 + e_z*sigma_z
# hoppings in x-direction; (a,0) means hopping form (i,j) to (i+1,j) syst[kwant.builder.HoppingKind((a,0), lat, lat)] = -t*sigma_z+ 1j*alpha*sigma_y/(2*a)
# hoppings in y-direction; (0,a) means hopping form (i,j) to (i,j+1) syst[kwant.builder.HoppingKind((0,a), lat, lat)] = -t*sigma_z + 1j*alpha*sigma_x/(2*a)
return syst
def make_system_leads(t=1, alpha=1, e_z=0.5, W=30, L=30, a=1): # this parameters can be overwritten in the 'main()' section lat = kwant.lattice.square(norbs=2) # norbs - the number of orbitals per site syst = kwant.Builder()
# Zeeman energy adds to the on-site term syst[(lat(x,y) for x in range(L) for y in range(W))] = 4*t*sigma_0 + e_z*sigma_z
# hoppings in x-direction; (a,0) means hopping form (i,j) to (i+1,j) syst[kwant.builder.HoppingKind((a,0), lat, lat)] = -t*sigma_z+ 1j*alpha*sigma_y/(2*a)
# hoppings in y-direction; (0,a) means hopping form (i,j) to (i,j+1) syst[kwant.builder.HoppingKind((0,a), lat, lat)] = -t*sigma_z + 1j*alpha*sigma_x/(2*a)
# create the leads (left and right) lead = kwant.Builder(kwant.TranslationalSymmetry((-a,0))) lead[(lat(0,j) for j in range(W))] = 4*t*sigma_0 + e_z*sigma_z lead[kwant.builder.HoppingKind((1,0), lat, lat)] = -t*sigma_0 - 1j*alpha*sigma_y/(2*a) # hoppings in x-direction lead[kwant.builder.HoppingKind((0,1), lat, lat)] = -t*sigma_0 + 1j*alpha*sigma_x/(2*a) # hoppings in y-direction
# attach the leads (left and right) syst.attach_lead(lead) syst.attach_lead(lead.reversed())
# create the leads (top and bottom) #lead2 = kwant.Builder(kwant.TranslationalSymmetry((0,-a))) #lead2[(lat(i,0) for i in range(L))] = 4*t*sigma_0 + e_z*sigma_z #lead2[kwant.builder.HoppingKind((1,0), lat, lat)] = -t*sigma_0 - 1j*alpha*sigma_y/(2*a) # hoppings in x-direction #lead2[kwant.builder.HoppingKind((0,1), lat, lat)] = -t*sigma_0 + 1j*alpha*sigma_x/(2*a) # hoppings in y-direction
# attach the leads (top and bottom) #syst.attach_lead(lead2) #syst.attach_lead(lead2.reversed())
return syst
def plot_dos_and_curves(dos, labels_to_data): pyplot.figure() pyplot.fill_between(dos[0], dos[1], label="DoS [a.u.]", alpha=0.5, color='gray') for label, (x, y) in labels_to_data: pyplot.plot(x, y, label=label, linewidth=2) pyplot.legend(loc=2, framealpha=0.5) pyplot.ylim(-1, 3) #! pyplot.xlabel("energy [t]") pyplot.ylabel("$\sigma [e^2/h]$") pyplot.show()
def main(random_vecs=True): num_moments = 400 # used in the method kwant.kpm.SpectralDensity # 'num_moments' is the number of moments, order of the KPM expansion. Mutually exclusive with 'energy_resolution' # default num_moments = 100
# dimensions and parameters for the both systems W = 30 L = 30 alpha = 1 e_z = 0.5
# calculate the band structure (infinite system: lead on the left and on the right hand side) syst_infty = make_system_leads(W=W, L=L, alpha=alpha, e_z=e_z) kwant.plot(syst_infty); # check that the system looks as intended (with the left and right leads attached) syst_infty = syst_infty.finalized() kwant.plotter.bands(syst_infty) # HERE I GET AN ERROR: 'Expecting an instance of InfiniteSystem.'
# the closed system to calculate conductivities using the method: kwant.kpm.conductivity() syst = make_system(W=W, L=L, alpha=alpha, e_z=e_z) syst = syst.finalized()
fam = syst.sites[0].family # assuming all sites from the same (sub)lattice area_per_orb = np.cross(*fam.prim_vecs) / fam.norbs # area that each of our vectors covers
#---------------------------------------------------------- if random_vecs: edge_width = min(L, W) / 4
def center(s): return (edge_width <= s.pos[0] < L - edge_width) and (edge_width <= s.pos[1] < W - edge_width)
center_vectors = kwant.kpm.RandomVectors(syst, where=center) # here kwant.kpm.LocalVectors can also be used, but will produce much more vectors one_vector = next(center_vectors) num_vectors = 10 # len(center_vectors)
else: # use local vectors edge_width = min(L, W) / 4
def center(s): return (edge_width <= s.pos[0] < L - edge_width) and (edge_width <= s.pos[1] < W - edge_width)
center_vectors = np.array(list(kwant.kpm.LocalVectors(syst, where=center))) one_vector = center_vectors[0] num_vectors = len(center_vectors) # len(x) - length of the object 'x' #----------------------------------------------------------
# we need to normalize the results by the area that each of our vectors covers norm = area_per_orb * np.linalg.norm(one_vector) ** 2
# check the area that the vectors cover kwant.plot(syst, site_color=np.abs(one_vector[::2]) + np.abs(one_vector[1::2]), cmap='Reds');
# object that represents the density of states for the system spectrum = kwant.kpm.SpectralDensity(syst, rng=0, num_moments=num_moments, num_vectors=num_vectors, vector_factory=center_vectors) # vector_factory - optional; it should contain (or yield) vectors of the size of the system # (here we use only the vector in the center of the system, thus we need to define this argument; # otherwise the KWNAT will randomize vectors for whole sample area)
energies, densities = spectrum()
# component 'xx' cond_xx = kwant.kpm.conductivity(syst, alpha='x', beta='x', num_moments=num_moments, rng=0, num_vectors=num_vectors, vector_factory=center_vectors) # component 'xy' cond_xy = kwant.kpm.conductivity(syst, alpha='x', beta='y', num_moments=num_moments, rng=0, num_vectors=num_vectors,vector_factory=center_vectors)
energies = cond_xx.energies cond_array_xx = np.array([cond_xx(e, temperature=0.01) for e in energies]) / norm cond_array_xy = np.array([cond_xy(e, temperature=0.01) for e in energies]) / norm
scale_to_plot = 1 / (edge_width) plot_dos_and_curves( (spectrum.energies, spectrum.densities / norm), [(r'Longitudinal conductivity $\sigma_{xx}$ (scaled)', (energies, cond_array_xx.real * scale_to_plot)), (r'Hall conductivity $\sigma_{xy}$', (energies, cond_array_xy.real))], )
# standard Python construct to execute 'main' if the program is used as a script if __name__ == '__main__': main()
Dear Anna,
I'm glad it works now.
but I am worry about the large noises in transverse conductivity beyond the energy gap
In an infinite 2d plane, the xy conductivity should decay to zero outside the gap. We can approximate this behavior with a finite system with the KPM expansion. Here, the number of moments (that is, the energy resolution) will broaden the delta peaks of your finite system.
In the example above there are two distict sources of noise outside the gap: a) due to the energy resolution and b) due to random vectors.
a) The `num_moments = 400` makes that we have enough energy resolution to resolve single states in your finite system, therefore you see the oscillation outside the gap, for example, using `LocalVectors`. In this case, you can lower the energy resolution, so that the system "appears" infinite at that broadenin. You can try with `num_moments=100`. A note on the `LocalVectors`: you just need to trace over one unit cell if the system has no disorder. So setting ``` if random_vecs: edge_width = min(L, W) / 4 ... else: # use local vectors edge_width = min(L, W) / 2 - 1 # only 1 unit cell ``` is enough to get standard infinite 2d plane result, approximated by the KPM expansion (it will use only 2 vectors, and will be very fast).
b) even if you have a lower energy resolution (`num_moments=100`) and you use random vectors over some central area of your system, you might see noise outside the gap. This is due to the large off-diagonal elements in the Kubo-Bastin operator. There are a few tricks to get around it, but still, random vectors introduce noise by definition. To ease the noise from random vectors you can either use more random vectors, or use a larger system / central area (computation scales linearly with the number of vectors *and* with the number of system sites). However, if you don't have a reason to average over a central area (such as disorder in your system), I advice to use `LocalVectors`.
According to the KWANT documentation, I found a method to calculate the band structure: kwant.physics.Bands() but it works only for infinite system.
Actually, a kwant system with leads attached can be used for both the bands and the kpm computation. The `kwant.kpm` module will only take the `FiniteSystem` part, and work just the same. To compute the bands you have to point which lead you want to compute: ``` kwant.physics.Bands(syst_infty.leads[0]) ``` Note that `syst_infty` contains the finite system and the two leads, but it is an instance of `FiniteSystem`, while `syst_infty.leads[0]` is an instance of `InfiniteSystem`, which is what `kwant.physics.Bands` expects.
Dear Pablo,
thank You for Your answer - it clarifies a lot, but at the same time a few questions in my head appear.
(At the beginning I'll mention that I'm more familiar with analytical calculation thus maybe my questions are trivial for the people familiar with numerical methodes, but I just want to understand some issues and figure out what is going in the results which I've obtained here.)
1. Now I have noticed that in Your previous reply You had changed the hoppings in the system in this way: <code>-t*sigma_0</code> -> <code>-t*sigma_z</code> and without this change the energy gap does not appear. I don't understand this change, because for me it creates a different system: this part originates from the kinetic term of the Hamiltonian and this term is defined with the identity matrix (sigma_0). Could You explain me why we can introduce that change?
2. Regarding the noises in xy-conductivity outside the gap in the finite system. You wrote about the tricks how to reduce this noise by changing the parametres of computation in KPM method - I checked, the noise is in fact smaller. But I want to find out how the xy-component will behave in the infinite system using the scattering-matrix method (will I get the results expected for the infinite system?). Thus, I've added the top and bottom leads and with the <code>smatrix.transmission()</code> method I've calulated the xx- and xy-conductance (based on 2.3.1 tutorial). Unfortunately, the results are unclear for me, because: a) adding the top and bottom leads to the system cancels the quantization in xx-conductance - I do not understand why adding leads perpendicular to the current flow does have the influence on the longitudinal conductance, b) the energy dependence of the both, the xx- and xy-conductance, increasing with oscillations - I don't understand this behaviour either, c) if I replace the hopping terms as You recommended, means <code>-t*sigma_0</code> -> <code>-t*sigma_z</code>, the results are unexpected either (0 in conductance till some energy, then random positive values).
3. The third thing, where I have a problem, is about the band structure. Taking a pen and a piece of paper: if I defined the system describing by the Hamiltonian which is the matrix 2x2, I will get two eigenvalues corresponding to two branches of the dispersion relation. Here, by calculating the band structure of the lead (describing by 2x2 matrix written out as on-site and hopping parameters), I get a lot of bands (which number increases with the size of the lead/system). I suppose that in this notation every unit cell/site produces a vector, and thus eigenvalue, and in that way we obtain so plentiful band structure, but it's not clear for me.
Below I attach the code for the 2. point.
Best regards, Anna
----------------------------------------------------- <code> import kwant from matplotlib import pyplot import tinyarray
# we define the identity matrix and Pauli matrices sigma_0 = tinyarray.array([[1,0], [0,1]]) sigma_x = tinyarray.array([[0,1], [1,0]]) sigma_y = tinyarray.array([[0,-1j], [1j,0]]) sigma_z = tinyarray.array([[1,0], [0,-1]])
def make_system(t=1.0, alpha=1, e_z=0.5, W=30, L=30, a=1): # for simplicity we set 'a' equal to 1 lat = kwant.lattice.square() syst = kwant.Builder()
# Zeeman energy adds to the on-site term syst[(lat(x,y) for x in range(L) for y in range(W))] = 4*t*sigma_0 + e_z*sigma_z
# hoppings in x-direction; (a,0) means hopping form (i,j) to (i+1,j) syst[kwant.builder.HoppingKind((a,0), lat, lat)] = -t*sigma_0 - 1j*alpha*sigma_y/(2*a) #syst[kwant.builder.HoppingKind((a,0), lat, lat)] = -t*sigma_z+ 1j*alpha*sigma_y/(2*a) #change from the Pablo's code
# hoppings in y-direction; (0,a) means hopping form (i,j) to (i,j+1) syst[kwant.builder.HoppingKind((0,a), lat, lat)] = -t*sigma_0 + 1j*alpha*sigma_x/(2*a) #syst[kwant.builder.HoppingKind((0,a), lat, lat)] = -t*sigma_z + 1j*alpha*sigma_x/(2*a) #change from the Pablo's code
lead = kwant.Builder(kwant.TranslationalSymmetry((-a,0))) lead[(lat(0,j) for j in range(W))] = 4*t*sigma_0 + e_z*sigma_z lead[kwant.builder.HoppingKind((1,0), lat, lat)] = -t*sigma_0 - 1j*alpha*sigma_y/(2*a) # hoppings in x-direction lead[kwant.builder.HoppingKind((0,1), lat, lat)] = -t*sigma_0 + 1j*alpha*sigma_x/(2*a) # hoppings in y-direction
# attach the left and right leads syst.attach_lead(lead) # left syst.attach_lead(lead.reversed()) # right
# create the leads (top and bottom) lead2 = kwant.Builder(kwant.TranslationalSymmetry((0,a))) lead2[(lat(i,0) for i in range(L))] = 4*t*sigma_0 + e_z*sigma_z lead2[kwant.builder.HoppingKind((1,0), lat, lat)] = -t*sigma_0 - 1j*alpha*sigma_y/(2*a) # hoppings in x-direction lead2[kwant.builder.HoppingKind((0,1), lat, lat)] = -t*sigma_0 + 1j*alpha*sigma_x/(2*a) # hoppings in y-direction
# attach the leads (top and bottom) syst.attach_lead(lead2) # top - ADDING THIS ELECTRODE changes the xx-component of conductivity tensor syst.attach_lead(lead2.reversed()) # bottom
# return the finalized system return syst
def plot_conductance(syst, energies, firstLead, secondLead): data = [] for energy in energies: smatrix = kwant.smatrix(syst, energy) data.append(smatrix.transmission(secondLead,firstLead)) # form lead 'firstLead' to 'secondLead'
pyplot.figure() pyplot.plot(energies, data) pyplot.xlabel("energy [t]") pyplot.ylabel(r"conductance [e$^2$/h]") pyplot.title("form " + str(firstLead) + " lead to " + str(secondLead) + " lead") pyplot.show()
def main(): syst = make_system()
kwant.plot(syst); # check that the system looks as intended
syst = syst.finalized()
# calculate the band structure (of the lead) # the results are the same for each (left,right,top,bottom) lead kwant.plotter.bands(syst.leads[0])
print('leads:\t 0 - left \n\t 1 - right \n\t 2 - top \n\t 3 - bottom')
# longitudinal conductivity (right-left = top-bottom) plot_conductance(syst, energies=[0.01*i-0.3 for i in range(100)], firstLead=0, secondLead=1) # from right to left #plot_conductance(syst, energies=[0.01*i-0.3 for i in range(100)], firstLead=2, secondLead=3) # from bottom to top
# transverse conductivity (reversible values, means top-left = left-top) plot_conductance(syst, energies=[0.01*i-0.3 for i in range(100)], firstLead=0, secondLead=2) # form top to left #plot_conductance(syst, energies=[0.01*i-0.3 for i in range(100)], firstLead=2, secondLead=0) # form left to top #plot_conductance(syst, energies=[0.01*i-0.3 for i in range(100)], firstLead=1, secondLead=2) # form top to right
# standard Python construct to execute 'main' if the program is used as a script if __name__ == '__main__': main() </code>
Dear Anna,
Thanks for the interest. I just want to point out that this mailing list should be to Kwant related questions, not purely physics questions. Altough, of course physics questions will arise in this context.
For purely physics questions you may try https://physics.stackexchange.com/, for example.
I answer below your questions.
1 - My apologies! I forgot to mention the change! I just saw it and changed it, because I assumed that you wanted to simulate a Chern insulator. Now I have read more carefully that you want Zeeman and Rashba spin-orbit interaction, following the Kwant tutorial.
2 - This system (Zeeman + Rashba) is not fully gapped, there is a spin-orbit gap that means an avoided crossing between the different spin bands, but there is no spectral gap. Therefore, the Hall conductivity (sigma_xy) is not quantized.
The 2 or 4 terminal conductance may differ from the sigma_xx, sigma_xy Kubo conductivtiy, depending on the transport setup, geometry, etc. The two terminal conductance (with perfectly metallic leads) should give you the sum of steps of value 1 for each band that has modes at the specific energy. This should be reproduced by the Kubo sigma_xx conductivity (or sigma_yy, it may differ). To simulate the Hall conductivity you need a 4 terminal setup (see Ref [1]), and this should be reproduced by the Kubo sigma_xy conductivtiy.
3 - I imagine that with pen and paper you derive the Bloch Hamiltonian of the system that depends on k_x, k_y. This means that you have the bandstructure for a 2D plane, and this bands from a set of surfaces in 3D, if you plot the energy along the z axis. In a system with 1D translational symmetry (a Kwant lead) you have many modes, depending on the with of the lead, transverse to the propagation direction. The bands of the lead are (more or less) like perpendicular cuts of the 3D dispersion relation of the bulk 2D system, but projected along the k-direction of the lead. See [2] for graphene nanoribbons, or the Kwant tutorial 2.4.
I hope this helps, and points you to the right sources to fill-in the gaps.
Best, Pablo
[1] Four-Terminal Phase-Coherent Conductance. M. Büttiker Phys. Rev. Lett. 57, 1761 https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.57.1761 [2] https://en.wikipedia.org/wiki/Graphene_nanoribbon
Dear Pablo,
I'm sorry, I don't know well the TB method so I couldn't distinguish weather my questions concern strictly physics or the way KWANT works.
Thank You again (also for the references) and I appreciate Your patience in answering all my questions. I'm sure it helps me to better understand what is going on in the KWANT world.
Best regards, Anna
participants (2)
-
anna.krzyzewska＠amu.edu.pl -
pablo.perez.piskunow＠gmail.com