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() On Mon, Apr 10, 2017 at 5:07 PM, John Eaves <johneaves02@gmail.com> wrote:

Hello,

My question is based on the reply to https://www.mail-archive.com/ kwant-discuss@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