%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Calculates Tunnel current by

% integrating over tip Fermi sphere

% Uses files C1.m and C2.m

% Must have Octave version 2.1.34

%

% Ian Appelbaum

% Harvard University

% October 2002

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;

m=5.6854e-16;

hbar=6.582e-16;

hb2o2m=3.8100e-16;

Ef=5; kf=sqrt(Ef/hb2o2m);

global V;

global dgap=10e-8;

ii=1;

for V=(1e-15:0.1:5)

VV(ii)=V;

II(ii)=quad('C1',sqrt(kf^2-V/hb2o2m),kf)+quad('C2',0,sqrt(kf^2-V/hb2o2m));

ii=ii+1;

end

gset nokey;

plot(VV,II,'k');

hold on;

plot(-VV,-II,'k');

xlabel('Voltage [V]')

ylabel('Current Density [nA/nm^2]')

hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% First integrand for tunnel current calculation

% valid for T=0 only because infinitesimal element is constant K_perp (not E).

% Ian Appelbaum

% Harvard University

% October 2002

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function el=C1(ktip)

m=5.11e5/(2.998e10)^2;

Ef_tip=5;

Wf_tip=5;

hbar=6.582e-16;

hb2o2m=hbar*hbar/(2.0*m);

global V;

global dgap;

a=-V/dgap;

b=Wf_tip+Ef_tip;

Etipz=hb2o2m*ktip*ktip;

kf=sqrt(Ef_tip/hb2o2m);

q=1.6e-19;

gamma=2*sqrt(2*m)/(3*a*hbar)*((a*dgap+b-Etipz)^1.5-(b-Etipz)^1.5);

T=exp(-2*gamma);

el=1e-5*q*hbar/m*ktip*T*pi*(kf^2-ktip^2); %1e-5 for conversion to nA/nm^2

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Second integrand for tunnel current calculation

% valid for T=0 only because infinitesimal element is constant K_perp (not E).

% Ian Appelbaum

% Harvard University

% October 2002

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function el=C2(ktip)

m=5.11e5/(2.998e10)^2;

Ef_tip=5;

Wf_tip=5;

hbar=6.582e-16;

hb2o2m=hbar*hbar/(2.0*m);

global dgap;

global V;

a=-V/dgap;

b=Wf_tip+Ef_tip;

Etipz=hb2o2m*ktip*ktip;

q=1.6e-19;

gamma=2*sqrt(2*m)/(3*a*hbar)*((a*dgap+b-Etipz)^1.5-(b-Etipz)^1.5);

T=exp(-2*gamma);

el=1e-5*q*hbar/m*ktip*T*pi*(V/hb2o2m);