Dear tkwant community,

Has anyone recently tried to simulate Josephson Junctions using tkwant?

I tried to follow the examples and tutorials, but several of them do not work in the last tkwant version. I am sending attached and you can find below the code that I am running but the current does not oscillate as a function of time and the IV curve is linear. Any suggestions on what may be wrong? Any help or comment is appreciated.

Here is the code:

import kwant
import tkwant
import tinyarray
import numpy as np
import matplotlib.pyplot as plt
import scipy.constants as ct
import time as comp_time
from mpi4py import MPI
comm=MPI.COMM_WORLD

start = comp_time.time()

   
tau_x = tinyarray.array([[0, 1], [1, 0]])
tau_y = tinyarray.array([[0, -1j], [1j, 0]])
tau_z = tinyarray.array([[1, 0], [0, -1]])
tau_0 = tinyarray.array([[1, 0], [0, 1]])
c_e = tinyarray.array([[1, 0], [0, 0]], complex)
c_h = tinyarray.array([[0, 0], [0, -1]], complex)


### parameters
Delta=0.02 #superconducting gap
barrier=0.9 #barrier in the N region
N=3 #size of the N region
S=10 #size of the S region
t=1
mu=t
I_time_avg=[]
V_values = np.linspace(-1.5*Delta, 1.5*Delta, 31)


def make_system(a=1):

    lat = kwant.lattice.square(norbs=2)
    syst = kwant.Builder()
   
    #### Define junction. ####
    ### S
    syst[(lat(x, 0) for x in range(-S, 0))] = (2 * t - mu) * tau_z + Delta * tau_x
    #### N
    syst[(lat(x, 0) for x in range(0, N))] = (2 * t - mu + barrier) * tau_z
    ### S
    syst[(lat(x, 0) for x in range(N, N+S))] = (2 * t - mu) * tau_z + Delta * tau_x


    # Hoppings
    syst[lat.neighbors()] = -t * tau_z
   
    ### add JJ phase
    syst[(lat(-1,0),lat(0,0))] =  -t*np.exp(1j*phase(time))*c_e -t*np.exp(-1j*phase(time))*c_h  ### JJ phase
   
    #### Define the leads. ####
    lead_left = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
    lead_left[lat(0, 0)] = (2 * t - mu) * tau_z + Delta * tau_x
    lead_left[lat.neighbors()] = -t * tau_z
    lead_right = lead_left.reversed()
    syst.attach_lead(lead_left)
    syst.attach_lead(lead_right)  
    return syst

def phase(time):
    ### int_{0}^{t}A(1-cos(pi*t'/tau))dt'
    def g(time,tau,A):
        return A*time-A*tau/np.pi*np.sin(np.pi*time/tau)
    if time<0:
        return 0
    elif time<tau:
        return g(time, tau=tau, A=V/2)
    else:
        return V*(time-tau)
       
for V in V_values:  
    tau=0.01/(np.abs(0.04*np.pi*Delta))
    tmax=4.01/(np.abs(0.04*np.pi*Delta))
     
    syst = make_system()
    lat = kwant.lattice.square(norbs=2)
   
    # Finalize the system.
    fsyst = syst.finalized()
   
    #bondaries condition
    boundaries=[tkwant.leads.MonomialAbsorbingBoundary(strength=8,degree=3,num_cells=1300)]
    #kwant.plot(fsyst)
    left_lead_interface = [(lat(1, 0), lat(0, 0))]


    current_operator = kwant.operator.Current(fsyst,where=left_lead_interface)

    # occupation -- for each lead
    occupation_l = tkwant.manybody.make_occupation(mu)
    occupation_r = tkwant.manybody.make_occupation(mu+V)
    occupation = [occupation_l, occupation_r]

    # Create the solver
    solver = tkwant.manybody.State(fsyst, tmax, occupation)
    I_time = []
   
    #get the current as a function of time
    real_time = np.linspace(0, tmax, 200)
    for time in real_time:
            solver.evolve(time)
            result = solver.evaluate(current_operator)
            I_time.append(result)
    data = np.column_stack((real_time, I_time))  
    I_time_avg.append(np.mean(I_time))        
   

end = comp_time.time()
print('Execution time is: '+str((end - start)/60)+' minutes')

# check current as a function of time
plt.figure()
plt.plot(real_time, I_time)

#plot IV curve
plt.figure()
plt.plot(V_values, np.asarray(I_time_avg)/(2*barrier*Delta), 'o')

data = np.column_stack((V_values, I_time_avg))
np.savetxt('IV_curve.txt', data)

Best regards,

Denise