superconductivity and current conservation
Dear all, in the tutorial https://kwant-project.org/doc/1/tutorial/superconductors one calculates the conductance from left lead to right lead, G_10, as G_00 = e^2/h (N-R_ee+R_he) assuming that, due to the current conservation, G_10 = G_00. My question is, what happens if we detach the superconducting lead L1. Intuitively, I would expect G_00 = 0. But using the formula above one obtains non-zero values, and making the finite superconducting region (the part of the scattering region) large enough, the subgap transport is exactly as if the lead L1 was attached. This means that the current conservation is violated, and something is wrong. Below is slightly modified code from the tutorial to generate a plot comparing G_00 for a system with and without lead L1. You can see, that if parameter L (the superconducting part of the scattering region) is large enough, the subgap transport is the same whether L1 is attached or not. Why is this happening, can someone point out the error? Best, Tibor import kwant import tinyarray import numpy as np %matplotlib inline from matplotlib import pyplot as plt tau_x = tinyarray.array([[0, 1], [1, 0]]) tau_y = tinyarray.array([[0, -1j], [1j, 0]]) tau_z = tinyarray.array([[1, 0], [0, -1]]) L=80 barrier=1.5 W=10 def make_system(a=1, W=W, L=L, barrier=barrier, barrierpos=(3, 4), mu=0.4, Delta=0.1, Deltapos=4, t=1.0, phs=True): lat = kwant.lattice.square(norbs=2) syst = kwant.Builder() syst[(lat(x, y) for x in range(0,Deltapos) for y in range(W))] = (4 * t - mu) * tau_z syst[(lat(x, y) for x in range(Deltapos, L) for y in range(W))] = (4 * t - mu) * tau_z + Delta * tau_x # The tunnel barrier syst[(lat(x, y) for x in range(barrierpos[0], barrierpos[1]) for y in range(W))] = (4 * t + barrier - mu) * tau_z # Hoppings syst[lat.neighbors()] = -t * tau_z sym_left = kwant.TranslationalSymmetry((-a, 0)) lead0 = kwant.Builder(sym_left, conservation_law=-tau_z, particle_hole=tau_y) lead0[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z lead0[lat.neighbors()] = -t * tau_z sym_right = kwant.TranslationalSymmetry((a, 0)) lead1 = kwant.Builder(sym_right) lead1[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z + Delta * tau_x lead1[lat.neighbors()] = -t * tau_z syst.attach_lead(lead0) syst.attach_lead(lead1) return syst syst = make_system() kwant.plot(syst) syst = syst.finalized() data = [] energies=[0.002 * i for i in range(0,100)] for energy in energies: smatrix = kwant.smatrix(syst, energy) # Conductance is N - R_ee + R_he data.append(smatrix.submatrix((0, 0), (0, 0)).shape[0] - smatrix.transmission((0, 0), (0, 0)) + smatrix.transmission((0, 1), (0, 0))) syst = make_system() del syst.leads[-1] # Check that the system looks as intended. kwant.plot(syst) # Finalize the system. syst = syst.finalized() data2 = [] for energy in energies: smatrix = kwant.smatrix(syst, energy) # Conductance is N - R_ee + R_he data2.append(smatrix.submatrix((0, 0), (0, 0)).shape[0] - smatrix.transmission((0, 0), (0, 0)) + smatrix.transmission((0, 1), (0, 0))) fig , (ax0) = plt.subplots(1, 1, sharey = True, figsize = [5, 5]) ax0.plot(energies, data, label = 'with lead 1') ax0.plot(energies, data2, label = 'without lead 1') ax0.set_xlabel("energy [t]") ax0.set_ylabel("conductance [e^2/h]") ax0.set_ylim(0,3.5) ax0.legend() fig.savefig('NS_conductance.png')
Hi Tibor, As soon as there is a superconducting region of any size, the charge can leave away as supercurrent, while the current carried by the quasiparticles isn't conserved anymore. Essentially all the superconductors in the system, both finite and infinite are an extra lead. I hope this answers your question, Anton On Thu, Aug 31, 2017 at 6:25 PM, Tibor Sekera <tibor.sekera@unibas.ch> wrote:
Dear all,
in the tutorial https://kwant-project.org/doc/1/tutorial/superconductors one calculates the conductance from left lead to right lead, G_10, as G_00 = e^2/h (N-R_ee+R_he) assuming that, due to the current conservation, G_10 = G_00.
My question is, what happens if we detach the superconducting lead L1. Intuitively, I would expect G_00 = 0. But using the formula above one obtains non-zero values, and making the finite superconducting region (the part of the scattering region) large enough, the subgap transport is exactly as if the lead L1 was attached. This means that the current conservation is violated, and something is wrong.
Below is slightly modified code from the tutorial to generate a plot comparing G_00 for a system with and without lead L1. You can see, that if parameter L (the superconducting part of the scattering region) is large enough, the subgap transport is the same whether L1 is attached or not.
Why is this happening, can someone point out the error?
Best, Tibor
import kwant import tinyarray import numpy as np %matplotlib inline from matplotlib import pyplot as plt
tau_x = tinyarray.array([[0, 1], [1, 0]]) tau_y = tinyarray.array([[0, -1j], [1j, 0]]) tau_z = tinyarray.array([[1, 0], [0, -1]])
L=80 barrier=1.5 W=10
def make_system(a=1, W=W, L=L, barrier=barrier, barrierpos=(3, 4), mu=0.4, Delta=0.1, Deltapos=4, t=1.0, phs=True): lat = kwant.lattice.square(norbs=2) syst = kwant.Builder()
syst[(lat(x, y) for x in range(0,Deltapos) for y in range(W))] = (4 * t - mu) * tau_z syst[(lat(x, y) for x in range(Deltapos, L) for y in range(W))] = (4 * t - mu) * tau_z + Delta * tau_x
# The tunnel barrier syst[(lat(x, y) for x in range(barrierpos[0], barrierpos[1]) for y in range(W))] = (4 * t + barrier - mu) * tau_z
# Hoppings syst[lat.neighbors()] = -t * tau_z sym_left = kwant.TranslationalSymmetry((-a, 0)) lead0 = kwant.Builder(sym_left, conservation_law=-tau_z, particle_hole=tau_y) lead0[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z lead0[lat.neighbors()] = -t * tau_z sym_right = kwant.TranslationalSymmetry((a, 0)) lead1 = kwant.Builder(sym_right) lead1[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z + Delta * tau_x lead1[lat.neighbors()] = -t * tau_z syst.attach_lead(lead0) syst.attach_lead(lead1) return syst
syst = make_system() kwant.plot(syst) syst = syst.finalized()
data = [] energies=[0.002 * i for i in range(0,100)] for energy in energies: smatrix = kwant.smatrix(syst, energy) # Conductance is N - R_ee + R_he data.append(smatrix.submatrix((0, 0), (0, 0)).shape[0] - smatrix.transmission((0, 0), (0, 0)) + smatrix.transmission((0, 1), (0, 0)))
syst = make_system() del syst.leads[-1] # Check that the system looks as intended. kwant.plot(syst) # Finalize the system. syst = syst.finalized() data2 = [] for energy in energies: smatrix = kwant.smatrix(syst, energy) # Conductance is N - R_ee + R_he data2.append(smatrix.submatrix((0, 0), (0, 0)).shape[0] - smatrix.transmission((0, 0), (0, 0)) + smatrix.transmission((0, 1), (0, 0)))
fig , (ax0) = plt.subplots(1, 1, sharey = True, figsize = [5, 5]) ax0.plot(energies, data, label = 'with lead 1') ax0.plot(energies, data2, label = 'without lead 1') ax0.set_xlabel("energy [t]") ax0.set_ylabel("conductance [e^2/h]") ax0.set_ylim(0,3.5) ax0.legend() fig.savefig('NS_conductance.png')
Hi Anton, I understand. Is this a feature of KWANT (and the way superconductivity is implemented) or of the formula for conductance? Then, for example, with KWANT we cannot calculate properly the conductance of a NSN system, where S is finite (not grounded) superconductor and N are leads. Because if we make S region very large we can get G_00 finite while G_10=0. In https://journals.aps.org/prb/abstract/10.1103/PhysRevB.53.16390 they claim that Eq.(31) is valid also for a 'floating' superconductor (not grounded finite piece), like in their Fig.2(b). In other words, how to treat properly with KWANT such a case, when we want to have finite piece of superconductor, that is not grounded? Best, Tibor ________________________________ From: anton.akhmerov@gmail.com [anton.akhmerov@gmail.com] on behalf of Anton Akhmerov [anton.akhmerov+kd@gmail.com] Sent: Thursday, August 31, 2017 6:58 PM To: Tibor Sekera Cc: kwant-discuss@kwant-project.org Subject: Re: [Kwant] superconductivity and current conservation Hi Tibor, As soon as there is a superconducting region of any size, the charge can leave away as supercurrent, while the current carried by the quasiparticles isn't conserved anymore. Essentially all the superconductors in the system, both finite and infinite are an extra lead. I hope this answers your question, Anton On Thu, Aug 31, 2017 at 6:25 PM, Tibor Sekera <tibor.sekera@unibas.ch<mailto:tibor.sekera@unibas.ch>> wrote: Dear all, in the tutorial https://kwant-project.org/doc/1/tutorial/superconductors one calculates the conductance from left lead to right lead, G_10, as G_00 = e^2/h (N-R_ee+R_he) assuming that, due to the current conservation, G_10 = G_00. My question is, what happens if we detach the superconducting lead L1. Intuitively, I would expect G_00 = 0. But using the formula above one obtains non-zero values, and making the finite superconducting region (the part of the scattering region) large enough, the subgap transport is exactly as if the lead L1 was attached. This means that the current conservation is violated, and something is wrong. Below is slightly modified code from the tutorial to generate a plot comparing G_00 for a system with and without lead L1. You can see, that if parameter L (the superconducting part of the scattering region) is large enough, the subgap transport is the same whether L1 is attached or not. Why is this happening, can someone point out the error? Best, Tibor import kwant import tinyarray import numpy as np %matplotlib inline from matplotlib import pyplot as plt tau_x = tinyarray.array([[0, 1], [1, 0]]) tau_y = tinyarray.array([[0, -1j], [1j, 0]]) tau_z = tinyarray.array([[1, 0], [0, -1]]) L=80 barrier=1.5 W=10 def make_system(a=1, W=W, L=L, barrier=barrier, barrierpos=(3, 4), mu=0.4, Delta=0.1, Deltapos=4, t=1.0, phs=True): lat = kwant.lattice.square(norbs=2) syst = kwant.Builder() syst[(lat(x, y) for x in range(0,Deltapos) for y in range(W))] = (4 * t - mu) * tau_z syst[(lat(x, y) for x in range(Deltapos, L) for y in range(W))] = (4 * t - mu) * tau_z + Delta * tau_x # The tunnel barrier syst[(lat(x, y) for x in range(barrierpos[0], barrierpos[1]) for y in range(W))] = (4 * t + barrier - mu) * tau_z # Hoppings syst[lat.neighbors()] = -t * tau_z sym_left = kwant.TranslationalSymmetry((-a, 0)) lead0 = kwant.Builder(sym_left, conservation_law=-tau_z, particle_hole=tau_y) lead0[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z lead0[lat.neighbors()] = -t * tau_z sym_right = kwant.TranslationalSymmetry((a, 0)) lead1 = kwant.Builder(sym_right) lead1[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z + Delta * tau_x lead1[lat.neighbors()] = -t * tau_z syst.attach_lead(lead0) syst.attach_lead(lead1) return syst syst = make_system() kwant.plot(syst) syst = syst.finalized() data = [] energies=[0.002 * i for i in range(0,100)] for energy in energies: smatrix = kwant.smatrix(syst, energy) # Conductance is N - R_ee + R_he data.append(smatrix.submatrix((0, 0), (0, 0)).shape[0] - smatrix.transmission((0, 0), (0, 0)) + smatrix.transmission((0, 1), (0, 0))) syst = make_system() del syst.leads[-1] # Check that the system looks as intended. kwant.plot(syst) # Finalize the system. syst = syst.finalized() data2 = [] for energy in energies: smatrix = kwant.smatrix(syst, energy) # Conductance is N - R_ee + R_he data2.append(smatrix.submatrix((0, 0), (0, 0)).shape[0] - smatrix.transmission((0, 0), (0, 0)) + smatrix.transmission((0, 1), (0, 0))) fig , (ax0) = plt.subplots(1, 1, sharey = True, figsize = [5, 5]) ax0.plot(energies, data, label = 'with lead 1') ax0.plot(energies, data2, label = 'without lead 1') ax0.set_xlabel("energy [t]") ax0.set_ylabel("conductance [e^2/h]") ax0.set_ylim(0,3.5) ax0.legend() fig.savefig('NS_conductance.png')
Hi Tibor, If you want to simulate a floating superconductor you only need to do an extra step. 1. Fix the voltages of all normal leads. 2. Assuming a certain voltage of a superconductor, compute current that is drained by it. 3. Numerically find a voltage when this current is zero by using any root-finding algorithm. Conceptually it isn't too hard, but it is somewhat more costly numerically. Also of course, if you have more than one superconductor and they are at different voltages, then this problem becomes much harder since there is no quasiparticle energy conservation anymore. Best, Anton On Thu, Aug 31, 2017 at 7:22 PM, Tibor Sekera <tibor.sekera@unibas.ch> wrote:
Hi Anton,
I understand. Is this a feature of KWANT (and the way superconductivity is implemented) or of the formula for conductance?
Then, for example, with KWANT we cannot calculate properly the conductance of a NSN system, where S is finite (not grounded) superconductor and N are leads. Because if we make S region very large we can get G_00 finite while G_10=0.
In https://journals.aps.org/prb/abstract/10.1103/PhysRevB.53.16390 they claim that Eq.(31) is valid also for a 'floating' superconductor (not grounded finite piece), like in their Fig.2(b).
In other words, how to treat properly with KWANT such a case, when we want to have finite piece of superconductor, that is not grounded?
Best, Tibor ------------------------------ *From:* anton.akhmerov@gmail.com [anton.akhmerov@gmail.com] on behalf of Anton Akhmerov [anton.akhmerov+kd@gmail.com] *Sent:* Thursday, August 31, 2017 6:58 PM *To:* Tibor Sekera *Cc:* kwant-discuss@kwant-project.org *Subject:* Re: [Kwant] superconductivity and current conservation
Hi Tibor,
As soon as there is a superconducting region of any size, the charge can leave away as supercurrent, while the current carried by the quasiparticles isn't conserved anymore. Essentially all the superconductors in the system, both finite and infinite are an extra lead.
I hope this answers your question, Anton
On Thu, Aug 31, 2017 at 6:25 PM, Tibor Sekera <tibor.sekera@unibas.ch> wrote:
Dear all,
in the tutorial https://kwant-project.org/doc/1/tutorial/superconductors one calculates the conductance from left lead to right lead, G_10, as G_00 = e^2/h (N-R_ee+R_he) assuming that, due to the current conservation, G_10 = G_00.
My question is, what happens if we detach the superconducting lead L1. Intuitively, I would expect G_00 = 0. But using the formula above one obtains non-zero values, and making the finite superconducting region (the part of the scattering region) large enough, the subgap transport is exactly as if the lead L1 was attached. This means that the current conservation is violated, and something is wrong.
Below is slightly modified code from the tutorial to generate a plot comparing G_00 for a system with and without lead L1. You can see, that if parameter L (the superconducting part of the scattering region) is large enough, the subgap transport is the same whether L1 is attached or not.
Why is this happening, can someone point out the error?
Best, Tibor
import kwant import tinyarray import numpy as np %matplotlib inline from matplotlib import pyplot as plt
tau_x = tinyarray.array([[0, 1], [1, 0]]) tau_y = tinyarray.array([[0, -1j], [1j, 0]]) tau_z = tinyarray.array([[1, 0], [0, -1]])
L=80 barrier=1.5 W=10
def make_system(a=1, W=W, L=L, barrier=barrier, barrierpos=(3, 4), mu=0.4, Delta=0.1, Deltapos=4, t=1.0, phs=True): lat = kwant.lattice.square(norbs=2) syst = kwant.Builder()
syst[(lat(x, y) for x in range(0,Deltapos) for y in range(W))] = (4 * t - mu) * tau_z syst[(lat(x, y) for x in range(Deltapos, L) for y in range(W))] = (4 * t - mu) * tau_z + Delta * tau_x
# The tunnel barrier syst[(lat(x, y) for x in range(barrierpos[0], barrierpos[1]) for y in range(W))] = (4 * t + barrier - mu) * tau_z
# Hoppings syst[lat.neighbors()] = -t * tau_z sym_left = kwant.TranslationalSymmetry((-a, 0)) lead0 = kwant.Builder(sym_left, conservation_law=-tau_z, particle_hole=tau_y) lead0[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z lead0[lat.neighbors()] = -t * tau_z sym_right = kwant.TranslationalSymmetry((a, 0)) lead1 = kwant.Builder(sym_right) lead1[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z + Delta * tau_x lead1[lat.neighbors()] = -t * tau_z syst.attach_lead(lead0) syst.attach_lead(lead1) return syst
syst = make_system() kwant.plot(syst) syst = syst.finalized()
data = [] energies=[0.002 * i for i in range(0,100)] for energy in energies: smatrix = kwant.smatrix(syst, energy) # Conductance is N - R_ee + R_he data.append(smatrix.submatrix((0, 0), (0, 0)).shape[0] - smatrix.transmission((0, 0), (0, 0)) + smatrix.transmission((0, 1), (0, 0)))
syst = make_system() del syst.leads[-1] # Check that the system looks as intended. kwant.plot(syst) # Finalize the system. syst = syst.finalized() data2 = [] for energy in energies: smatrix = kwant.smatrix(syst, energy) # Conductance is N - R_ee + R_he data2.append(smatrix.submatrix((0, 0), (0, 0)).shape[0] - smatrix.transmission((0, 0), (0, 0)) + smatrix.transmission((0, 1), (0, 0)))
fig , (ax0) = plt.subplots(1, 1, sharey = True, figsize = [5, 5]) ax0.plot(energies, data, label = 'with lead 1') ax0.plot(energies, data2, label = 'without lead 1') ax0.set_xlabel("energy [t]") ax0.set_ylabel("conductance [e^2/h]") ax0.set_ylim(0,3.5) ax0.legend() fig.savefig('NS_conductance.png')
participants (2)
-
Anton Akhmerov -
Tibor Sekera