I tried using the python module Kwant for quantum transport calculation

Introduction

Investigating the properties of a substance by giving a stimulus (for example, an electric field) to the substance and observing the response (for example, an electric current) is the basis of condensed matter physics, but it is quite difficult to actually perform simulations and numerical calculations.

Therefore, we will introduce the open source python module "Kwant" that can easily analyze quantum transport phenomena.

What is Kwant?

Kwant is a tight binding model solver that is still being actively developed by researchers at the Grenoble Center for Nuclear Research and Delft University of Technology. It is also Paper. The original HP document is very rich, and you can easily reproduce and expand physically interesting calculation examples at hand just by clicking according to it. If it is a system that can be written with a tight binding model, it can be defined very easily including the electrode terminals, and it comprehensively covers the wave function, density of states, and transportation calculation using the Landauer-Buttiker formula and Kubo formula.

I had previously thought about using Kwant as part of my research, but I had to postpone the introduction because there were almost no articles or communities in Japanese and I didn't know anything about python at that time. It was. I don't even know how to pronounce it. Now that I have some free time and have become accustomed to python, I finally tried using it, so I would like to introduce how it looks.

• There is a great possibility that errors may be included in both the understanding of physics and the understanding of modules. Please see at your own risk. Those who can read English may want to refer to the head family.

Mesoscopic

Kwant mainly analyzes physical systems called mesoscopic systems. The point is that mesoscopic devices are small at the atomic level, and quantum effects such as confinement and interference dominate transport. It is especially interesting that the terminals affect the quantum state and transport.

As for the review of physics, if you read Supriyo Datta's "Electronic transport in mesoscopic systems" or "Quantum transport basics" series of "Basics of nanoscale physical properties", which is also recommended by the head family, I will not go into it here, but focus on implementation to place. First, try it. Then think.

Installation

It is assumed that the environment for python has been built. If you create a python environment using Anaconda

conda install -c conda-forge kwant

You can install it in one shot. The current latest version is v1.4. It seems that pip will have a hard time, so please refer to the Procedure of the head family.

Example 1: Quantization of conductance

First, as the simplest example, attach terminals to both ends of a two-dimensional rectangular device and calculate the conductivity between them. I refer to here and chapters 2.2.2 and 2.2.3 of the original document.

In Kwant, we first define (Build) the grid shape and system.

#Import base module
import kwant
import numpy as np
import scipy as sp
import math
import matplotlib.pyplot as plt
%matplotlib inline

#Parameters
W=10 #Device width
L=30 #Device length
t=1 #Recent hopping
a=1 #Lattice constant

#Lattice shape and system definition
lat = kwant.lattice.square(a)
syst = kwant.Builder()

This defines the grid shape lat and the system `` `syst. As it is, syst``` is an empty state defined only by Gawa, so we will fill in the contents. The content here is the Hamiltonian, which is a tight binding of the system.

#Onsite potential
syst[(lat(x, y) for x in range(L) for y in range(W))] = 0

#Recent hopping(y direction)
syst[((lat(x, y),lat(x,y-1)) for x in range(L) for y in range(1,W))] = -t
#(x direction)
syst[((lat(x, y),lat(x-1,y)) for x in range(1,L) for y in range(W))] = -t

This will store the tight binding Hamiltonian with a width of W and a length of L in syst. If you define hopping only in one direction, it also defines hopping in the opposite direction so that Hamiltonian becomes Hermitian. `for x in range (L)` is a comprehension of python.

Next, define the leads (electrode terminals). Assuming that the lead extends in one direction semi-infinitely, specify that direction. Here, we first define the terminal to be attached to the left side of the system, so specify the (-1,0) direction. Lead onsite potential and hopping can be defined in the same way as `syst```, but if you define only points with x-coordinates 0 and 1, TranslationalSymmetry ((-a, 0)) By ``, it repeats periodically to the left.

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

#Onsite potential
lead[(lat(0, j) for j in range(W))] = 0
#Recent hopping
lead[((lat(0, j), lat(0, j - 1) ) for j in range(1,W))] = -t
lead[((lat(1, j), lat(0, j    ) ) for j in range(0,W))] = -t

Now we have defined a lead that extends infinitely to the left, but it's still independent of syst, so we'll stick it together.

syst.attach_lead(lead)

Only this.

Then attach the lead to the opposite right side as well. To flip a defined lead over and attach it to the other side

syst.attach_lead(lead.reversed())

Only this. it's simple.

By now, the rectangular sample and the two leads `` `leadshould be stored in syst```. Kwant is also excellent in system visualization, and you can see at a glance what kind of system is designed.

kwant.plot(syst)

The blue dots are the grid points of the sample, the red dots are the grid points of the leads, and the hopping is defined where there are lines. You can see that the leads are neatly attached to both ends of the sample. When you are satisfied, complete the design.

syst = syst.finalized()

This completes the system design including the terminals.

Next, let's calculate the conductivity between the terminals of the created system. According to the Landauer formula, the conductivity between terminals is $ e ^ 2 / h in terms of the quantum mechanical transmission probability between terminals (a plane wave is incident and a wave function is connected in the middle to obtain various values). Calculate the conductivity because it is equal to multiplying by $. The transmittance can be calculated by giving the finalized `syst``` and the energy of the incident wave to the kwant.smatrix () `class.

conductance = []
energies = np.linspace(-4,-3,50)
for energy in energies:
    #Scattering matrix(S matrix)Definition of
    smatrix = kwant.smatrix(syst, energy)
    #Calculates and stores the transparency probability from read 0 to 1.
    conductance.append(smatrix.transmission(1, 0))
#plot
plt.plot(energies, conductance)
plt.grid()
plt.xlabel("energy/t")
plt.ylabel("conductance($e^2/h$)")
plt.show()

it's simple. Leads are numbered 0, 1, 2 ... in the order of attach_lead, so 1 and 0 are specified this time. Looking at the results obtained, it can be seen that the conductivity increases stepwise as the energy of the incident wave is increased. This is an effect peculiar to the mesoscopic system called conductance quantization, and its origin is the quantum mechanical confinement effect in the y direction.

Looking back, I was able to perform the Hamiltonian definition of the sample, the Hamiltonian definition of the lead, the adhesion of the lead, and the calculation of the conductivity with about 20 lines of code. It's nice and easy to do a similar calculation from scratch, considering that it takes a lot of physical and coding effort to calculate the surface Green's function of the leads and define the lead-sample hopping.

Bonus (spin freedom)

The unit matrix of spin space is defined as s0 = np.array ([[1,0], [0,1]])` `` and written as "` t in this example. Can be expressed as spin degeneracy by replacing with t * s0, and the conductance is doubled. Similarly, if Pauli matrices are defined, spin-orbit interaction and Zeeman magnetic field can be given.

Example 2: Quantum Hall effect

As an example of dealing with non-uniform on-site potential and hopping in addition to Example 1, consider the quantum Hall effect of a disturbed system. By the way, you need to be careful when handling a large number of terminals. I refer to various documents of the head family (same below)

Assuming that the magnetic field is applied vertically in the two-dimensional plane where the device exists, attach 6 terminals according to the arrangement of the hole bars that are often seen in experiments. Import of kwant and numpy is the same as the above example, so it will be omitted. First is the definition of the shape of the device. Only the additional parameters are commented.

W=15
L_lead=6 #Hole lead width
M_lead=6 #Distance from the edge of the hole lead
L=30
t=1
a=1
V=0.05 #Random potential strength

lat = kwant.lattice.square(a)
syst = kwant.Builder()

Next, we define the Hamiltonian of the sample part. Onsite potential is given appropriately using normal random numbers. Again, we will consider the nearest-only hopping and take into account the Peierls phase due to the perpendicular magnetic field.

#Onsite random potential
for x in range(L):
    for y in range(W):
        syst[lat(x, y)] = np.random.randn()*V

#(xi, yi)From(xj, yj)Function that gives hopping to
def hopping(site_i, site_j, phi):
    xi, yi = site_i.pos
    xj, yj = site_j.pos
    return -t*np.exp(-0.5j * phi * (xi - xj) * (yi + yj))

#Recent hopping
syst[lat.neighbors()] = hopping

Unlike Example 1, I used `lat.neighbors ()` to specify the Hamiltonian of the reed. This is useful when defining the closest hopping of a site point that has already been defined, as you don't have to bother to specify the coordinates. By the way, if n = 2 is given as a parameter, the next proximity hopping can be defined in the same way.

The function hopping that defines hopping contains the site information of the move source and destination in the first and second arguments, and the subsequent arguments are interpreted as parameters given from the outside. Here, the parameter `phi```, which indicates the magnitude of the magnetic field, has no value at this point, and can be given a value from the outside after finalizing syst```. For example, if you want to see the parameter dependency of some response such as conductivity, you don't have to redefine `` syst``` for each parameter value and finalize it (see below).

Then define the 6 leads in order. Note that the four terminals other than the two terminals attached to both ends are `kwant.TranslationalSymmetry ((0, a))` because they extend in the y direction. The code for each ``` lat.neighbors ()` `` written in a straightforward manner is also commented out.

lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
lead[(lat(0, j) for j in range(W))] = 0
# lead[((lat(0, j), lat(0, j - 1) ) for j in range(1,W))] = hopping
# lead[((lat(1, j), lat(0, j    ) ) for j in range(0,W))] = hopping
lead[lat.neighbors()] = hopping

lead_hall24 = kwant.Builder(kwant.TranslationalSymmetry((0, a)))
lead_hall24[(lat(j, 0) for j in range(M_lead,M_lead+L_lead))] = 0
# lead_hall24[((lat(j,0), lat(j-1,0) ) for j in range(M_lead+1,M_lead+L_lead))] = hopping
# lead_hall24[((lat(j,1), lat(j  ,0) ) for j in range(M_lead,M_lead+L_lead))] = hopping
lead_hall24[lat.neighbors()] = hopping

lead_hall35 = kwant.Builder(kwant.TranslationalSymmetry((0, a)))
lead_hall35[(lat(j, 0) for j in range(L-M_lead-L_lead,L-M_lead))] = 0
# lead_hall35[((lat(j,0), lat(j-1,0) ) for j in range(L-M_lead-L_lead+1,L-M_lead))] = hopping
# lead_hall35[((lat(j,1), lat(j  ,0) ) for j in range(L-M_lead-L_lead,L-M_lead))] = hopping
lead_hall35[lat.neighbors()] = hopping

#Installation of all 6 leads. Pay attention to the order.
syst.attach_lead(lead)
syst.attach_lead(lead.reversed())
syst.attach_lead(lead_hall24)
syst.attach_lead(lead_hall35)
syst.attach_lead(lead_hall24.reversed())
syst.attach_lead(lead_hall35.reversed())

So far, the whole device has been defined. Visualize and finalize

kwant.plot(syst)
syst = syst.finalized()

I hand-numbered the leads.

From here, we will calculate values such as Hall conductivity and Hall resistance. The basic function is `smatrix.conductance_matrix ()`. It can be used to calculate $ G $ of the proportional relation $ I = GV $ between the current $ I $ and the voltage $ V $ from the S matrix. Since $ I and V $ are defined for each of the 6 terminals, they are vectors of $ 6 $ dimension respectively, so $ G $ is a $ 6 \ times 6 $ matrix. When the state where the current is passed only between terminals 0 and 1 is expressed as $ I_0 = -I_1 = I, \ \ I_2 = I_3 = I_4 = I_5 = 0 $, the vertical resistance and the hall resistance are $ R_ {long}, respectively. You can write = \ frac {V_4-V_5} {I}, \ R_ {Hall} = \ frac {V_3-V_5} {I} $, so it seems that you should first find the inverse of the matrix $ G $.

However, the matrix $ G $ is known to be downgraded due to the conservation law of current and the indefiniteness of the origin of voltage. So here we assume $ V_5 = 0 $ and ignore the rightmost column of $ G $. Then $ G $ becomes a $ 6 \ times 5 $ matrix, and the inverse matrix problem becomes overdetermined. So I'll ignore the bottom row and treat $ G $ as a $ 5 \ times 5 $ matrix.

The above-mentioned parameter phi can be given as a dictionary when defining the S matrix (in this example, it was not necessary to assign `` `phi```, so it is necessary to go out. I didn't have it ...)

#A function to obtain the vertical resistance, Hall resistance, conductivity, and Hall conductivity from the S matrix.
from numpy.linalg import LinAlgError
def calc_conductances2(smatrix):
    #Calculate the conductance matrix G
    gmat = smatrix.conductance_matrix()
    #Exception handling when the matrix is accidentally singular
    try:
        Rmat = np.linalg.inv(gmat[:5,:5])
    except LinAlgError:
        R_long = np.nan
        R_Hall = np.nan
        sigma_xx = 0
        sigma_xy = 0
    else:
        #V_5=I put it as 0, so R_long=V_4/I
        R_long = Rmat[4,0]-Rmat[4,1]
        R_Hall = Rmat[3,0]-Rmat[3,1]
        sigma_xx = R_long/(R_long**2+R_Hall**2)
        sigma_xy =-R_Hall/(R_long**2+R_Hall**2)
    return R_long , R_Hall , sigma_xx ,sigma_xy 

phi=2*np.pi*1/11 #The value of the parameter phi.
R_Halls =[]
R_longs =[]
sigma_xxs =[]
sigma_xys =[]
doss=[]
energies = np.linspace(-4,0,50)
for energy in energies:
    #Calculation of S matrix. Give parameters as dictionary type.
    smatrix = kwant.smatrix(syst, energy, params=dict(phi=phi) )
    R_long, R_Hall, sigma_xx, sigma_xy = calc_conductances(smatrix)
    R_Halls.append(R_Hall)
    R_longs.append(R_long)
    sigma_xxs.append(sigma_xx)
    sigma_xys.append(sigma_xy)    

    #Also calculate the density of states.
    dos = kwant.ldos(syst, energy, params=dict(phi=phi) ).mean()
    doss.append(dos)

If you plot conductivity and Hall conductivity as a function of energy

plt.plot(energies, doss, label="dos")
plt.plot(energies, sigma_xxs, label="sigma_xx")
plt.plot(energies, sigma_xys, label="sigma_xy")
plt.grid()
plt.ylim(-1,5)
plt.xlim(-4,-0.8)
plt.xlabel("energy/t")
plt.legend()
plt.show()

It can be confirmed that the Hall conductivity is quantized to an integral multiple of $ e ^ 2 / h $, and the density of states and the longitudinal conductivity are similar to the Landau level.

Now you can make the on-site potential and hopping non-uniform, give parameters from the outside, and add a lot of terminals.

Bonus (Kubo formula)

You can calculate the Hall conductivity more easily using the Kubo formula without attaching a terminal with the following code. I borrowed only the `` `hopping``` function from the top, but it can be calculated in about 10 lines. Wow!

W=100
L=W
t=1
a=1

lat = kwant.lattice.square(a)
syst = kwant.Builder()

syst[(lat(x, y) for x in range(L) for y in range(W))] = 0
syst[lat.neighbors()] = hopping
syst = syst.finalized()

sigma_xy = kwant.kpm.conductivity(syst, alpha='x', beta='y',params=dict(phi=2*np.pi/11))
conductivities = [sigma_xy(mu=mu, temperature=0)/L/W for mu in np.linspace(-4,0,50)]
plt.plot(np.linspace(-4,0,50),conductivities)
plt.grid()

However, if the default settings are used, the quantization will be loose and the calculation accuracy will be difficult. It seems that the kernel polynomial method (KPM) is used for the calculation here, and the accuracy can be improved by finely adjusting the calculation parameters. You can increase the device size, but it will take a lot of time. In this article, due to the author's knowledge and motivation limits, I will leave the details to the original document ...

Bonus (gauge)

In v1.4, a function was added to calculate the Peierls phase of the gauge field and the hopping caused by it when a magnetic field is applied. See the documentation for details.

Example 3: Anomalous quantum Hall effect

As a final example, let's look at the anomalous quantum Hall effect in the Haldane model (special case). When the honeycomb lattice-shaped spinless fermions perform pure imaginary order proximity hopping, the Hall conductivity is quantized without an external magnetic field.

So far, we have only considered rectangular devices with square grids arranged in a circle, but here we will create a device with honeycomb grids arranged in a circle. First, the grid shape is changed and defined as `` `lat```. it's simple.

#Parameters
a=1
t=1.0
tt=0.5 #Next proximity hopping
L_lead=10

#Define a honeycomb grid. norbs is a parameter that represents the number of orbitals per site (other than the final current distribution, it works without it). A on the secondary grid,Name it b.
lat = kwant.lattice.honeycomb(a, norbs=1, name=['a', 'b'])
#Here together
syst = kwant.Builder()

Next, let's define the shape of the device. To define the shape, use a function that returns a bool value for the coordinates. By specifying one point inside the defined shape, the simply connected area can be filled with the site.

#Define shape
def circle(pos):
    x, y = pos
    return x**2 + y**2 < 100
#Specify one point inside the shape. here(0,0).. Onsite potential is 0
syst[lat.shape(circle, (0, 0))] = 0

Hopping can be defined as well. Since the sub-lattices a and b are defined, the closest contact (= next proximity of the honeycomb lattice) can be easily specified only on it. Here, the next proximity hopping is a pure imaginary number.

syst[lat.neighbors()] = -t
syst[lat.a.neighbors()] = -1j*tt
syst[lat.b.neighbors()] = 1j*tt

The definition of lead is similar

lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
def lead_shape(pos):
        x, y = pos
        return (-L_lead/2 < y < L_lead/2)
lead[lat.shape(lead_shape,(0,0))] = 0

lead[lat.neighbors()] = -t
lead[lat.a.neighbors()] = -1j*tt
lead[lat.b.neighbors()] = 1j*tt
syst.attach_lead(lead)
syst.attach_lead(lead.reversed())

Let's take a look at the shape of the device.

kwant.plot(syst)
syst = syst.finalized()

Certainly I was able to attach leads to both ends of the circular honeycomb grid device. If you look closely, you can also see the black line that represents the next proximity hopping. The reed has a slightly strange shape as a result of being able to satisfy the semi-infinite translational symmetry to the left and right, but the shape of this side is likely to change in the armchair direction and the zigzag direction. (I haven't tried)

You can also add multiple terminals to calculate the quantization of the Hall conductivity, but since it is a big deal, let's visualize the quantum Hall state by another method. First, let's calculate the transmittance and density of states between the two terminals. The code is almost the same as in Example 1.

trans=[]
doss=[]
energies = np.linspace(-4,4,100,endpoint=False)
for energy in energies:
    smatrix = kwant.smatrix(syst, energy )
    trans.append(smatrix.transmission(0,1)) 
    dos = kwant.ldos(syst,energy ).mean()
    doss.append(dos)

plt.plot(energies, trans, label="conductance")
#It is roughly multiplied by 10 for good appearance.
plt.plot(energies, [d*10 for d in doss], label="dos")
plt.grid()
plt.legend()
plt.xlabel("energy/t")
plt.show()

Looking at this, we can see that the conductance is quantized when the incident energy is near 0. In other words, it means that there is only one conduction channel without scattering.

Next, looking at the density of states, it becomes smaller in this energy region. I used the spatial averaging of the ``` kwant.ldos ()` `` function to calculate the density of states, but let's plot the distribution (local density of states) in real space before averaging. ..

local_dos = kwant.ldos(syst,energy=0.1 )
kwant.plotter.map(syst,local_dos)

Then you can see that the density of states is finite only at the edges of the device. That is, the current can only flow through the edges of the device, and the interior is insulating. This represents the edge conduction due to the edge state peculiar to the quantum Hall state.

You can also actually plot the current distribution.

psi = kwant.wave_function(syst,energy=0.1)(0)
J = kwant.operator.Current(syst)
current=sum(J(p) for p in psi)
kwant.plotter.current(syst,current,min_linewidth=.01,density=0.3)

It is a vector field representation of how the electrons incident from lead 0 are flowing. Certainly, it has reached lead 1 on the right end while being localized on the edge of the device.

Now you know that you can design odd-shaped devices and grid-shaped devices.

Summary

I looked at the quantum conduction calculation by Kwant focusing on the two-dimensional system. There are many functions that I haven't understood the behavior yet, and many functions that I haven't mastered, so it seems that I can do more things by reading the document. In addition to the documentation, you can find a mailing list for discussion and a sample jupyter notebook from the Kwant homepage, so please play with it.

In addition to the ones illustrated here, there are other interesting things that can be done.

• Calculation of spin current and heat flow
• Calculation of band structure of periodic system
• Superconductor, Josephson junction?
• Exchange coupling with static magnetism by local Zeeman magnetic field
• 3D system, other grid system there is. If you have the time and motivation, I will check it.