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 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