[Kwant] Finite temperature, finite bias conductance, and units of energy

[John Eaves]

Hello,

My question is based on the reply to
https://www.mail-archive.com/kwant-***@kwant-project.org/msg00983.html.
Here it is written that one can extend the KWANT zero temperature
infinitesimal source drain bias conductance to finite bias, finite
temperature by numerical integration. The answer refers the user to the
Datta book to see how exactly this is done.



However, I have looked at the Datta book and I cannot figure out how to
implement it in KWANT. My problem has a few parts, so below I will build it
up step by step so that one can see if perhaps I make a mistake somewhere.
However, I can start with a very simple question:



In order to account for finite temp, you need the fermi dirac distribution
(or its derivative). Inside the distribution we have k_B*T; how do you
rescale k_B*T to agree with the units of the kinetic term? The kinetic term
is t = hbar^2/(2m a^2). If I set t = 1, what do I do to kB*T? I don’t
understand this as a (the lattice constant) enters the expression, but
changing a and keeping t fixed does not change any of my results. So
somehow I need to relate hbar^2/2m with kB*T as these are the fixed
constants? It is the fact that a enters into t but not into kb*T that makes
it difficult to understand I think.



Now, for my main question, lets start from the basics. If I am not
mistaken, in the default zero T zero bias case KWANT essentially calculates
G = 2e^2/h * T(E_fermi), where T(E_fermi) is the quantity calculated from

smatrix = kwant.smatrix(sys, energy)

T = smatrix.transmission(1, 0)



This makes sense in the context of linear response. Now if we want to go to
finite temp, finite bias we take the expression from Datta



I = 2e/h \int{dE T(E) (f_L(E) - f_R(E))}



where L,R is the left and right leads at potentials \mu_L, \mu_R abd f_L is
then the fermi-dirac distribution of the left lead.



Now this is the current; we don’t want the current, we want the
conductance. So let us set \mu_L = E_f – eV_sd/2 and \mu_R = E_f + eV_sd/2
(which I suppose one should be able to do without loss of generality, but
correct me if I am wrong). If we then want to find the conductance we
simply take the derivative with respect to V_sd giving us



G = 2e^2/h \int{dE T(E) d(f_L(E) - f_R(E))/dVsd}



One simple limit here is to take T = 0, in which case the fermi-dirac
distribution becomes a step function, its derivative a delta function, and
we (not forgetting to carry the minus sign and the œ) get that

G = e^2/h (T(E_f+V_sd/2) + T(E_f-V_sd/2)



So if I am not mistaken, one can calculate the finite-bias zero-temperature
conductance without numerical integration; you simply have to evaluate
smatrix.transmission(1,0) at two different energies.This seems correct as
it recovers the standard KWANT expression for zero bias.



Now for finite temperature and zero bias things are different of course.
One can taylor expand (f_L(E) - f_R(E)) as df_L/dV_sd *eV_sd = -df_L/dE
*eV_sd so that



G = e^2/h \int{dE T(E) -df_L(E) /dE}


Now here things get a bit tricky: the interval is over all energies. Of
course the fermi-dirac derivative is strongly peaked around E_fermi, but
numerically this will still be annoying, will it not? If I were to
numerically calculate the above current I’d need to evaluate
smatrix.transmission(1, 0) for a huge range of energies, right?

Could you tell me something about how one would go about this in a smart
way? Do you use the trapezoid method? How do you keep it from taking ages?
What I would come up with myself is something like



data = []

integ = []

for energy in energies:

for erange in np.linspace(0.5*energy,1.5*energy,num=100):

smatrix = kwant.smatrix(sys, erange,args=[params])

integ.append(smatrix.transmission(1,0)*-dfd(erange,energy,T))


data.append(1.5*energy/100*(integ[0]+integ[100-1]+2*sum(integ[1:100-2])))

integ = []



where in the above dfd(energy,mu,temp) is the derivative of the fermi-dirac
distribution. But this seems horribly inefficient, and I don't see how this
would recover the KWANT limit (as we get some delta-like peak around E_f),
but perhaps that is too much to ask.



Finally we return to the full expression for finite bias and finite T:



G = 2e^2/h \int{dE T(E) d(f_L(E) - f_R(E))/dVsd}



I suppose evaluating this is equally difficult/easy as finite T with zero
bias; all that changes is that one now has to calculate two derivatives of
the fermi-dirac distribution which will form some sort of window around
E_fermi. Here my question this repeats itself; how does one compute such
an integral efficiently?



Kind regards,

John

[Abbout Adel]

Dear John,



At T=0, you do not need the derivative over v: in fact, F_L(E)-F_R(E)=
V_sd
(because at T=0, F(E) is 0 or 1 depending on E>V or E<V)
So you will have:
I = 2e^2/h \int{dE T(E) } V_sd

and therefore, your conductance is G=2e^2/h \int{dE T(E) } . The
integration is done from -V_sd/2 to V_sd/2 (So, your result is beyond the
linear response theory).

Case T different from zero (V_sd=0):
To include temperature, you need to have an idea about the order of
energies in your system (for example E_F, the broadening Gamma of
resonances, KT ...).
To compare with experiment, we need realistic values. The starting point is
the hopping 't'.
t=hbar^2/(2m* a^2) where m* is the effective mass. For example m*=0.067 m_e
in 2DEG. If we take a=2nm, we find t=147 mev.

We know that: 1mev=11.6 Kelvin. With this relation, we can get the
corresponding temperature for the energies we consider.

Temperature has a small effect if the transmission is not varying in the
range KT, and has a very noticeable effect on resonances. So, if you have a
resonance of width Gamma, you need an interval of integration of length
~5KT (or more ) and a step in the energies ~min(KT,Gamma) /10 (or
smaller). So in total you need around 50 KT/min(KT,Gamma) energies.

I put below a small example: beta=1/KT

regards,
Adel



# In[3]:

import kwant
from numpy import cosh,trapz
from matplotlib import pyplot


def make_system(a=1, t=1.0,tc=0.1, Pot=-1.8):

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

sys[(lat(x, 0) for x in range(-1,2) )] = 0
sys[lat.neighbors()] = -t

sys[lat(0,0)]=Pot
sys[(lat(0,0),lat(+1,0))]=-tc
sys[(lat(0,0),lat(-1,0))]=-tc


lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
lead[lat(0, 0)] = 0
lead[lat.neighbors()] = -t

sys.attach_lead(lead)
sys.attach_lead(lead.reversed())

return sys


def main():

sys = make_system()
kwant.plot(sys)

# derivative of Fermi-Dirac
def dF(energy,Ef,beta):
return 1/4.* beta/(cosh(beta*(energy-Ef)/2.))**2

sys = sys.finalized()


#KT=1/beta so you need an energy step much smaller than KT.if
beta=2500, KT=0.0004 ! so your step needs to be
# around 0.0004/15 or even smaller, that is why I put 1/(15*beta).
#more over, you need a total range of the energies = few KT here it is
300/(15*beta)=20*KT which is OK

def get_trans(Ef,beta):
energies=[ Ef+i/(15*beta) for i in range (-300,300)]
trans=[]
for energy in energies:
smatrix = kwant.smatrix(sys, energy)
T=smatrix.transmission(0, 1)
trans.append(T)
return trans,energies

#Fermi_E=[ Ef+i/(15*beta) for i in range (-150,150)] #Here I chose
the energies of the Fermi Dirac in the same ensemple as the one of the
transmission
#this way I do not have to recalculate the transmissions again.
# I stop at 150 step in order to let more energies around the Fermi
energies

def Use_temp(beta,Ef=-1.818): # Ef=-1.818 is the energy at which I
have a resonance
trans,energies=get_trans(Ef,beta)
couple_T_E=list(zip(trans,energies)) # here, I have a list
of (T,E)

Conductance=[]
Fermi_E=[ Ef+i/(15.*beta) for i in range (-150,150)] #Here
I chose the energies of the Fermi Dirac in the same ensemple as the one of
the transmission
#this way I do not have to recalculate the transmissions again.
# I stop at 150 in order to let more energies around the Fermi
energies

for Ef in Fermi_E:
Tr=[T*dF(energy,Ef=Ef,beta=beta) for T,energy in
couple_T_E]
g=trapz(Tr, energies,dx=1./(15*beta)) # this is the
conductance at one energy.
Conductance.append(g)
return Conductance, Fermi_E,trans,energies

fig = pyplot.figure()
ax = pyplot.subplot(111)
for beta in [2000.,1500.,500.,300.,200.]:
print(beta)
Conductance, Fermi_E,trans,energies= Use_temp(beta)
pyplot.plot(Fermi_E,Conductance,label= "beta="+str(beta),linewidth=3)
#pyplot.legend()

pyplot.plot(energies, trans,label= "beta="+r'$\infty$',linewidth=3)
ax.legend()
pyplot.xlabel("energy [t]")
pyplot.ylabel("conductance [e^2/h]")
#pyplot.savefig('fig.png')
pyplot.show()





if __name__ == '__main__':
main()

Post by John Eaves
Hello,
My question is based on the reply to https://www.mail-archive.com/
one can extend the KWANT zero temperature infinitesimal source drain bias
conductance to finite bias, finite temperature by numerical integration.
The answer refers the user to the Datta book to see how exactly this is
done.
However, I have looked at the Datta book and I cannot figure out how to
implement it in KWANT. My problem has a few parts, so below I will build it
up step by step so that one can see if perhaps I make a mistake somewhere.
In order to account for finite temp, you need the fermi dirac distribution
(or its derivative). Inside the distribution we have k_B*T; how do you
rescale k_B*T to agree with the units of the kinetic term? The kinetic term
is t = hbar^2/(2m a^2). If I set t = 1, what do I do to kB*T? I don’t
understand this as a (the lattice constant) enters the expression, but
changing a and keeping t fixed does not change any of my results. So
somehow I need to relate hbar^2/2m with kB*T as these are the fixed
constants? It is the fact that a enters into t but not into kb*T that makes
it difficult to understand I think.
Now, for my main question, lets start from the basics. If I am not
mistaken, in the default zero T zero bias case KWANT essentially calculates
G = 2e^2/h * T(E_fermi), where T(E_fermi) is the quantity calculated from
smatrix = kwant.smatrix(sys, energy)
T = smatrix.transmission(1, 0)
This makes sense in the context of linear response. Now if we want to go
to finite temp, finite bias we take the expression from Datta
I = 2e/h \int{dE T(E) (f_L(E) - f_R(E))}
where L,R is the left and right leads at potentials \mu_L, \mu_R abd f_L
is then the fermi-dirac distribution of the left lead.
Now this is the current; we don’t want the current, we want the
conductance. So let us set \mu_L = E_f – eV_sd/2 and \mu_R = E_f + eV_sd/2
(which I suppose one should be able to do without loss of generality, but
correct me if I am wrong). If we then want to find the conductance we
simply take the derivative with respect to V_sd giving us
G = 2e^2/h \int{dE T(E) d(f_L(E) - f_R(E))/dVsd}
One simple limit here is to take T = 0, in which case the fermi-dirac
distribution becomes a step function, its derivative a delta function, and
we (not forgetting to carry the minus sign and the œ) get that
G = e^2/h (T(E_f+V_sd/2) + T(E_f-V_sd/2)
So if I am not mistaken, one can calculate the finite-bias
zero-temperature conductance without numerical integration; you simply have
to evaluate smatrix.transmission(1,0) at two different energies.This seems
correct as it recovers the standard KWANT expression for zero bias.
Now for finite temperature and zero bias things are different of course.
One can taylor expand (f_L(E) - f_R(E)) as df_L/dV_sd *eV_sd = -df_L/dE
*eV_sd so that
G = e^2/h \int{dE T(E) -df_L(E) /dE}
Now here things get a bit tricky: the interval is over all energies. Of
course the fermi-dirac derivative is strongly peaked around E_fermi, but
numerically this will still be annoying, will it not? If I were to
numerically calculate the above current I’d need to evaluate
smatrix.transmission(1, 0) for a huge range of energies, right?
Could you tell me something about how one would go about this in a smart
way? Do you use the trapezoid method? How do you keep it from taking ages?
What I would come up with myself is something like
data = []
integ = []
smatrix = kwant.smatrix(sys, erange,args=[params])
integ.append(smatrix.transmission(1,0)*-dfd(erange,energy,T))
data.append(1.5*energy/100*(integ[0]+integ[100-1]+2*
sum(integ[1:100-2])))
integ = []
where in the above dfd(energy,mu,temp) is the derivative of the
fermi-dirac distribution. But this seems horribly inefficient, and I don't
see how this would recover the KWANT limit (as we get some delta-like peak
around E_f), but perhaps that is too much to ask.
G = 2e^2/h \int{dE T(E) d(f_L(E) - f_R(E))/dVsd}
I suppose evaluating this is equally difficult/easy as finite T with zero
bias; all that changes is that one now has to calculate two derivatives of
the fermi-dirac distribution which will form some sort of window around
E_fermi. Here my question this repeats itself; how does one compute such
an integral efficiently?
Kind regards,
John