I have been trying to evaluate current through a josephson junction using the green function formalism using the following formula: I=2 (e Kb T)/hbar(sum over matsubara freq.)(H_ij G_ij-H_ji G_ji) To calculate current I have mounted two virtual leads on the slabs where I want to evaluate the current using the following code:
def mount_vlead(sys, vlead_interface, norb): dim = len(vlead_interface)*norb zero_array = np.zeros((dim, dim), dtype=float) def selfenergy_func(energy, args=()): return zero_array
vlead = kwant.builder.SelfEnergyLead(selfenergy_func, vlead_interface,()) sys.leads.append(vlead)
where vlead_interace was supplied using the following code:
v1=kwant.Builder() v1[(lat(nx_1,i) for i in range(ny))]=(2*t-mu_n)*tau_z v1[lat.neighbors()]=-t*tau_z v2=kwant.Builder() v2[(lat(nx_2,i) for i in range(ny))]=(2*t-mu_n)*tau_z v2[lat.neighbors()]=-t*tau_z mount_vlead(syst,v2.sites(),2)
Here nx_1 and nx_2 are the two points where the points where I have mounted virtual leads. Then I have used the greens_function solver in kwant to evaluate green function for a given matsubara frequency. The issue that i am facing is that I am getting identical diagonal elements in green_function.submatrix(v1,v2) and green_function.submatrix(v2,v1), thus when taking trace it gives me zero as the hopping matrix is same.
Dear Mohit,
The Caroli et al formula gives the tunneling current in a system connected to electrodes (leads in our terminology). A lead is an infinite system which has the important characteristics that it exhibits a complex self energy with a negative imaginary part.
If you had constructed your leads normally, this imaginary part would have been present automatically in the conduction region (the bands). Since you put fictitious leads with a zero self energy, you will obtain no current.
Change your self energy to something like: -1j* Gamma/2 where Gamma is a norbs x norbs definite positive matrix. It should work like this.
I hope this helps, Adel
On Wed, Jun 10, 2020 at 4:22 AM mohit.gupta9607@gmail.com wrote:
I have been trying to evaluate current through a josephson junction using the green function formalism using the following formula: I=2 (e Kb T)/hbar(sum over matsubara freq.)(H_ij G_ij-H_ji G_ji) To calculate current I have mounted two virtual leads on the slabs where I want to evaluate the current using the following code:
def mount_vlead(sys, vlead_interface, norb): dim = len(vlead_interface)*norb zero_array = np.zeros((dim, dim), dtype=float) def selfenergy_func(energy, args=()): return zero_array
vlead = kwant.builder.SelfEnergyLead(selfenergy_func,vlead_interface,()) sys.leads.append(vlead)
where vlead_interace was supplied using the following code:
v1=kwant.Builder() v1[(lat(nx_1,i) for i in range(ny))]=(2*t-mu_n)*tau_z v1[lat.neighbors()]=-t*tau_z v2=kwant.Builder() v2[(lat(nx_2,i) for i in range(ny))]=(2*t-mu_n)*tau_z v2[lat.neighbors()]=-t*tau_z mount_vlead(syst,v2.sites(),2)
Here nx_1 and nx_2 are the two points where the points where I have mounted virtual leads. Then I have used the greens_function solver in kwant to evaluate green function for a given matsubara frequency. The issue that i am facing is that I am getting identical diagonal elements in green_function.submatrix(v1,v2) and green_function.submatrix(v2,v1), thus when taking trace it gives me zero as the hopping matrix is same.
Dear Abbout, Thanks for your reply. The idea was to evaluate green's function for two particular slices of the system I basically followed this discussion: https://rb.gy/v6bkpo Which says to add zero self energy to obtain green's function everywhere. So my next question is if I don't give zero self energy won't my Green's function depend on the choice of Gamma?
I think the other issue is that I am not 100% sure about how to read the output of the greens_function(syst,energy).
Dear Mohit,
The green's function as implemented in kwant returns the elements of this latter only between the sites connected to the leads. The idea behind mounting fictitious leads with zero self-energy on some specific sites inside the system is just to tell kwant that I want the Green's function for these sites too.
Now, in the answer I gave to you last time, I supposed that your system is finite and had no leads (other than the fictitious leads you mounted). That is why I said: the current will be zero if the fictitious self energy is zero. Is your system finite?
I also do not see from where did you get this formula.((H_ij G_ij-H_ji G_ji)) ? The usual formula for the current is Im[H_ij (G Gamma G^\dagger)]_ij, Gamma = 2Im[Self_energy]. So please check also your formula.
Most importantly, Kwant has functions which return local current which are very efficient. Why do you not use them?
Now, I come back to your question about the meaning of what the function of greensfunction returns:
If the mounted lead is number 3, and you want the greens function between sites site1, site2 than you do: mount_vlead(syst,[site1,site2],2) G=greens_function(syst.finalized(),energy=-1.8, in_leads=[3],out_leads=[3]).data
This will return the matrix G=([G11,G12],[G21,G22]) if the number of orbitals =1 otherwise it will double in case norbs=2.
G12 id the greens function element from site 1 to site 2.
I hope this helps. Adel
On Wed, Jun 17, 2020 at 3:02 AM mohit.gupta9607@gmail.com wrote:
Dear Abbout, Thanks for your reply. The idea was to evaluate green's function for two particular slices of the system I basically followed this discussion: https://rb.gy/v6bkpo Which says to add zero self energy to obtain green's function everywhere. So my next question is if I don't give zero self energy won't my Green's function depend on the choice of Gamma?
I think the other issue is that I am not 100% sure about how to read the output of the greens_function(syst,energy).
Dear Abbout, I apologize for not being clear in my question. My system is infinite (that is there are leads attached to the system) and I mounted virtual leads at points where I wanted kwant to evaluate green's function. Sorry about that.
The formula that I am using is from the supplementary material of the following paper: https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.120.047702 It is possible that I am misinterpreting the formula as given in the paper as am very new to the subject.
I want to evaluate critical current in the Josephson Junction by maximizing current-phase dependency, for that I need total current in the system as a function of phase difference phi, I am not sure how to do that with local current operator in kwant.
Thank you for clarifying the output of green's function I am now sure that I was reading the output wrong which must have been the primary reason for getting zero.
I again apologize if I am not clear about what I am saying.
Hi,
I am not sure if kwant can really compute supercurrent in JJ. I think kwant superconducting simulation does not compute pairing potentials self-consistently. As such, the current continuity condition is not guaranteed and the current is not accurate. I would suggest that you first plot current as a function of position (in 1D or 2D or whatever geometry you are considering) to see if the continuity condition is satisfied.
Best,
Amrit
On Wed, Jun 17, 2020 at 8:51 AM mohit.gupta9607@gmail.com wrote:
Dear Abbout, I apologize for not being clear in my question. My system is infinite (that is there are leads attached to the system) and I mounted virtual leads at points where I wanted kwant to evaluate green's function. Sorry about that.
The formula that I am using is from the supplementary material of the following paper: https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.120.047702 It is possible that I am misinterpreting the formula as given in the paper as am very new to the subject.
I want to evaluate critical current in the Josephson Junction by maximizing current-phase dependency, for that I need total current in the system as a function of phase difference phi, I am not sure how to do that with local current operator in kwant.
Thank you for clarifying the output of green's function I am now sure that I was reading the output wrong which must have been the primary reason for getting zero.
I again apologize if I am not clear about what I am saying.
Thanks. I will try to check that. I am basically trying to reproduce some results from the papers I have been reading which has the implementation done in kwant.
Dear Mohit,
I see the formula you are mentioning using the matsubara frequencies.
To use that expression, you need to use imaginary energies. The steps are;
Sites1=[lat(5,j) for j in range(W)] # the sites where to mount lead 2 Sites2=[lat(6,j) for j in range(W)] # the sites where to mount lead 3
mount_vlead(sys, Sites1, 4) # number of orbitals =4. Here we mount lead number 2 mount_vlead(sys, Sites2, 4) # number of orbitals =4. Here we mount lead number 3 syst=sys.finalized() G12=greens_function(syst, energy=-1.8*1j, in_leads=[3],out_leads=[2], check_hermiticity=False).data G21=greens_function(syst, ,energy=-1.8*1j, in_leads=[2],out_leads=[3], check_hermiticity=False).data H12=syst.hamiltonian_submatrix(to_sites=Sites1, from_sites=Sites2) H21=syst.hamiltonian_submatrix(to_sites=Sites2, from_sites=Sites1)
Notices that, we had to use "check_hermiticity=False" because for these complex energies, the system is not hermitian
Trace(np.dot(H21, G12)-np.dot(H12,G21) )) # check if the order is correct.
I hope this helps. Adel
On Wed, Jun 17, 2020 at 6:51 PM mohit.gupta9607@gmail.com wrote:
Dear Abbout, I apologize for not being clear in my question. My system is infinite (that is there are leads attached to the system) and I mounted virtual leads at points where I wanted kwant to evaluate green's function. Sorry about that.
The formula that I am using is from the supplementary material of the following paper: https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.120.047702 It is possible that I am misinterpreting the formula as given in the paper as am very new to the subject.
I want to evaluate critical current in the Josephson Junction by maximizing current-phase dependency, for that I need total current in the system as a function of phase difference phi, I am not sure how to do that with local current operator in kwant.
Thank you for clarifying the output of green's function I am now sure that I was reading the output wrong which must have been the primary reason for getting zero.
I again apologize if I am not clear about what I am saying.
Dear Abbout, Thank you for your helpful reply it clears up a lot. I tried to obtain H12 as you mentioned: H12=syst.hamiltonian_submatrix(to_sites=Sites1, from_sites=Sites2) But it gives me the error that "an integer is required". Although that is not a major issue I could try to construct H12 manually. Thanks a lot for your help.
Dear Mohit,
Sorry, I wrote it without testing. Indeed, we need to provide the indices of the sites and not the sites themselves. Please change Sites1, and Sites2 in the hamiltonian_submatrix as follows
Sites1=[syst.id_by_site[lat(5, j)] for j in range(W)]
I hope this helps. Adel
On Thu, Jun 18, 2020 at 12:52 AM mohit.gupta9607@gmail.com wrote:
Dear Abbout, Thank you for your helpful reply it clears up a lot. I tried to obtain H12 as you mentioned: H12=syst.hamiltonian_submatrix(to_sites=Sites1, from_sites=Sites2) But it gives me the error that "an integer is required". Although that is not a major issue I could try to construct H12 manually. Thanks a lot for your help.
Thanks. The code is working now. Though I am still getting zero current I will try to implement hamiltonian more carefully and see. Thanks for all your help.
-
Abbout Adel -
Amrit Poudel -
mohit.gupta9607＠gmail.com