diff options
Diffstat (limited to 'figs/fowler-nordheim.fig')
| -rw-r--r-- | figs/fowler-nordheim.fig/FN_base.jl | 170 | ||||
| -rw-r--r-- | figs/fowler-nordheim.fig/Makefile | 29 | ||||
| -rw-r--r-- | figs/fowler-nordheim.fig/asymptotic.gnuplot | 59 | ||||
| -rw-r--r-- | figs/fowler-nordheim.fig/asymptotic.jl | 46 |
4 files changed, 304 insertions, 0 deletions
diff --git a/figs/fowler-nordheim.fig/FN_base.jl b/figs/fowler-nordheim.fig/FN_base.jl new file mode 100644 index 0000000..af2a1ee --- /dev/null +++ b/figs/fowler-nordheim.fig/FN_base.jl @@ -0,0 +1,170 @@ +# fractional power with an arbitrary branch cut +function pow(x,a,cut) + if(angle(x)/cut<=1) + return(abs(x)^a*exp(1im*angle(x)*a)) + else + return(abs(x)^a*exp(1im*(angle(x)-sign(cut)*2*pi)*a)) + end +end + +# asymptotic airy functions +# specify a branch cut for the fractional power +function airyai_asym(x,cut) + if(abs(real(pow(x,3/2,cut)))<airy_threshold) + return(exp(2/3*pow(x,3/2,cut))*airyai(x)) + else + ret=0 + for n in 0:airy_order + ret+=gamma(n+5/6)*gamma(n+1/6)*(-3/4)^n/(4*pi^(3/2)*factorial(n)*pow(x,3*n/2+1/4,cut)) + end + return ret + end +end +function airyaiprime_asym(x,cut) + if(abs(real(pow(x,3/2,cut)))<airy_threshold) + return(exp(2/3*pow(x,3/2,cut))*airyaiprime(x)) + else + ret=0 + for n in 0:airy_order + ret+=gamma(n+5/6)*gamma(n+1/6)*(-3/4)^n/(4*pi^(3/2)*factorial(n))*(-1/pow(x,3*n/2-1/4,cut)-(3/2*n+1/4)/pow(x,3*n/2+5/4,cut)) + end + return ret + end +end + +# solutions of (-\Delta+V-ip)phi=0 +# assume that p has an infinitesimal real part (and adjust the branch cuts appropriately) +function phi(p,x,E,V) + return(airyai_asym(exp(-1im*pi/3)*(E^(1/3)*x-E^(-2/3)*(V-1im*p)),pi)) +end +function dphi(p,x,E,V) + return(exp(-1im*pi/3)*E^(1/3)*airyaiprime_asym(exp(-1im*pi/3)*(E^(1/3)*x-E^(-2/3)*(V-1im*p)),pi)) +end +function eta(p,x,E,V) + return(exp(-1im*pi/3)*airyai_asym(-(E^(1/3)*x-E^(-2/3)*(V-1im*p)),pi/2)) +end +function deta(p,x,E,V) + return(-exp(-1im*pi/3)*E^(1/3)*airyaiprime_asym(-(E^(1/3)*x-E^(-2/3)*(V-1im*p)),pi/2)) +end + +# Laplace transform of psi +# assume that p has an infinitesimal real part (and adjust the branch cuts appropriately) +# for example, (1im*p-V)^(3/2) becomes pow(1im*p-V,3/2,-pi/2) because when 1im*p is real negative, its square root should be imaginary positive +function f(p,x,k0,E,V) + T=2im*k0/(1im*k0-sqrt(V-k0*k0)) + R=T-1 + + if x>=0 + C2=-1im*T/(pow(-1im*p,1/2,pi/2)*phi(p,0,E,V)-dphi(p,0,E,V))*((sqrt(V-k0*k0)+pow(-1im*p,1/2,pi/2))/(-1im*p+k0*k0)-2im*E^(-1/3)*pi*quadgk(y -> (pow(-1im*p,1/2,pi/2)*eta(p,0,E,V)-deta(p,0,E,V))*phi(p,y,E,V)*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2)-E^(-1)*pow(1im*p-V,3/2,-pi/2))),0,Inf)[1]) + FT=2*E^(-1/3)*pi*(quadgk(y -> phi(p,x,E,V)*eta(p,y,E,V)*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*x+E^(-2/3)*(1im*p-V),3/2,-pi/2)-pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2))),0,x)[1]+quadgk(y -> eta(p,x,E,V)*phi(p,y,E,V)*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2)-pow(E^(1/3)*x+E^(-2/3)*(1im*p-V),3/2,-pi/2))),x,Inf)[1]) + main=C2*phi(p,x,E,V)*exp(2im/3*(pow(E^(1/3)*x+E^(-2/3)*(1im*p-V),3/2,-pi/2)-E^(-1)*pow(1im*p-V,3/2,-pi/2)))+T*FT + + # subtract the contribution of the pole, which will be added back in after the integration + pole=psi_pole(x,k0,E,V)/(p+1im*k0*k0) + return(main-pole) + else + C1=-1im*T*((sqrt(V-k0*k0)*phi(p,0,E,V)+dphi(p,0,E,V))/(-1im*p+k0*k0)/(pow(-1im*p,1/2,pi/2)*phi(p,0,E,V)-dphi(p,0,E,V))+E^(-1/3)*quadgk(y -> phi(p,y,E,V)/(pow(-1im*p,1/2,pi/2)*phi(p,0,E,V)-dphi(p,0,E,V))*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2)-E^(-1)*pow(1im*p-V,3/2,-pi/2))),0,Inf)[1]) + FI=-1im*exp(1im*k0*x)/(-1im*p+k0*k0) + FR=-1im*exp(-1im*k0*x)/(-1im*p+k0*k0) + main=C1*exp(pow(-1im*p,1/2,pi/2)*x)+FI+R*FR + + # subtract the contribution of the pole, which will be added back in after the integration + pole=psi_pole(x,k0,E,V)/(p+1im*k0*k0) + return(main-pole) + end +end +# its derivative +function df(p,x,k0,E,V) + T=2im*k0/(1im*k0-sqrt(V-k0*k0)) + R=T-1 + + if x>=0 + C2=-1im*T/(pow(-1im*p,1/2,pi/2)*phi(p,0,E,V)-dphi(p,0,E,V))*((sqrt(V-k0*k0)+pow(-1im*p,1/2,pi/2))/(-1im*p+k0*k0)-2im*E^(-1/3)*pi*quadgk(y -> (pow(-1im*p,1/2,pi/2)*eta(p,0,E,V)-deta(p,0,E,V))*phi(p,y,E,V)*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2)-E^(-1)*pow(1im*p-V,3/2,-pi/2))),0,Inf)[1]) + dFT=2*E^(-1/3)*pi*(quadgk(y -> dphi(p,x,E,V)*eta(p,y,E,V)*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*x+E^(-2/3)*(1im*p-V),3/2,-pi/2)-pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2))),0,x)[1]+quadgk(y -> deta(p,x,E,V)*phi(p,y,E,V)*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2)-pow(E^(1/3)*x+E^(-2/3)*(1im*p-V),3/2,-pi/2))),x,Inf)[1]) + main=C2*dphi(p,x,E,V)*exp(2im/3*(pow(E^(1/3)*x+E^(-2/3)*(1im*p-V),3/2,-pi/2)-E^(-1)*pow(1im*p-V,3/2,-pi/2)))+T*dFT + + # subtract the contribution of the pole, which will be added back in after the integration + pole=dpsi_pole(x,k0,E,V)/(p+1im*k0*k0) + return(main-pole) + else + C1=-1im*T*((sqrt(V-k0*k0)*phi(p,0,E,V)+dphi(p,0,E,V))/(-1im*p+k0*k0)/(pow(-1im*p,1/2,pi/2)*phi(p,0,E,V)-dphi(p,0,E,V))+E^(-1/3)*quadgk(y -> phi(p,y,E,V)/(pow(-1im*p,1/2,pi/2)*phi(p,0,E,V)-dphi(p,0,E,V))*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2)-E^(-1)*pow(1im*p-V,3/2,-pi/2))),0,Inf)[1]) + dFI=k0*exp(1im*k0*x)/(-1im*p+k0*k0) + dFR=-k0*exp(-1im*k0*x)/(-1im*p+k0*k0) + main=C1*pow(-1im*p,1/2,pi/2)*exp(pow(-1im*p,1/2,pi/2)*x)+dFI+R*dFR + + # subtract the contribution of the pole, which will be added back in after the integration + pole=dpsi_pole(x,k0,E,V)/(p+1im*k0*k0) + return(main-pole) + end +end + +# psi (returns t,psi(x,t)) +function psi(x,k0,E,V,p_npoints,p_cutoff) + fft=fourier_fft(f,x,k0,E,V,p_npoints,p_cutoff) + # add the contribution of the pole + for i in 1:p_npoints + fft[2][i]=fft[2][i]+psi_pole(x,k0,E,V)*exp(-1im*k0*k0*fft[1][i]) + end + return(fft) +end +# its derivative +function dpsi(x,k0,E,V,p_npoints,p_cutoff) + fft=fourier_fft(df,x,k0,E,V,p_npoints,p_cutoff) + # add the contribution of the pole + for i in 1:p_npoints + fft[2][i]=fft[2][i]+dpsi_pole(x,k0,E,V)*exp(-1im*k0*k0*fft[1][i]) + end + return(fft) +end + +# compute Fourier transform by sampling and fft +function fourier_fft(A,x,k0,E,V,p_npoints,p_cutoff) + fun=zeros(Complex{Float64},p_npoints) + times=zeros(p_npoints) + + # prepare fft + for i in 1:p_npoints + fun[i]=p_cutoff/pi*A(1im*(-p_cutoff+2*p_cutoff*(i-1)/p_npoints),x,k0,E,V) + times[i]=(i-1)*pi/p_cutoff + end + + ifft!(fun) + + # correct the phase + for i in 2:2:p_npoints + fun[i]=-fun[i] + end + return([times,fun]) +end + +# asymptotic value of psi +function psi_pole(x,k0,E,V) + if x>=0 + return(1im*phi(-1im*k0*k0,x,E,V)*2*k0/(1im*k0*phi(-1im*k0*k0,0,E,V)+dphi(-1im*k0*k0,0,E,V))*exp(2im/3*(pow(E^(1/3)*x+E^(-2/3)*(k0*k0-V),3/2,-pi/2)-E^(-1)*pow(k0*k0-V,3/2,-pi/2)))) + else + return((1im*k0*phi(-1im*k0*k0,0,E,V)-dphi(-1im*k0*k0,0,E,V))/(1im*k0*phi(-1im*k0*k0,0,E,V)+dphi(-1im*k0*k0,0,E,V))*exp(-1im*k0*x)+exp(1im*k0*x)) + end +end +function dpsi_pole(x,k0,E,V) + if x>=0 + return(1im*dphi(-1im*k0*k0,x,E,V)*2*k0/(1im*k0*phi(-1im*k0*k0,0,E,V)+dphi(-1im*k0*k0,0,E,V))*exp(2im/3*(pow(E^(1/3)*x+E^(-2/3)*(k0*k0-V),3/2,-pi/2)-E^(-1)*pow(k0*k0-V,3/2,-pi/2)))) + else + return(-1im*k0*(1im*k0*phi(-1im*k0*k0,0,E,V)-dphi(-1im*k0*k0,0,E,V))/(1im*k0*phi(-1im*k0*k0,0,E,V)+dphi(-1im*k0*k0,0,E,V))*exp(-1im*k0*x)+1im*k0*exp(1im*k0*x)) + end +end + +# current +function J(ps,dps) + return(2*imag(conj(ps)*dps)) +end + +# complete computation of the current +function current(x,k0,E,V,p_npoints,p_cutoff) + ps=psi(x,k0,E,V,p_npoints,p_cutoff) + dps=dpsi(x,k0,E,V,p_npoints,p_cutoff) + Js=zeros(Complex{Float64},p_npoints) + for i in 1:p_npoints + Js[i]=J(ps[2][i],dps[2][i]) + end + return(Js) +end diff --git a/figs/fowler-nordheim.fig/Makefile b/figs/fowler-nordheim.fig/Makefile new file mode 100644 index 0000000..a6e5e53 --- /dev/null +++ b/figs/fowler-nordheim.fig/Makefile @@ -0,0 +1,29 @@ +PROJECTNAME=asymptotic + +PDFS=$(addsuffix .pdf, $(PROJECTNAME)) +TEXS=$(addsuffix .tikz.tex, $(PROJECTNAME)) + +all: $(PDFS) + +$(PDFS): $(addsuffix .dat, $(PROJECTNAME)) + gnuplot $(patsubst %.pdf, %.gnuplot, $@) > $(patsubst %.pdf, %.tikz.tex, $@) + pdflatex -jobname $(basename $@) -file-line-error $(patsubst %.pdf, %.tikz.tex, $@) + +asymptotic.dat: + julia asymptotic.jl > asymptotic.dat + +install: $(PDFS) + cp $^ $(INSTALLDIR)/ + +clean-aux: + rm -f $(addsuffix .aux, $(PROJECTNAME)) + rm -f $(addsuffix .log, $(PROJECTNAME)) + +clean-dat: + rm -f $(addsuffix .tikz.tex, $(PROJECTNAME)) + rm -f short-time.dat + +clean-tex: + rm -f $(PDFS) + +clean: clean-dat clean-aux clean-tex diff --git a/figs/fowler-nordheim.fig/asymptotic.gnuplot b/figs/fowler-nordheim.fig/asymptotic.gnuplot new file mode 100644 index 0000000..3296380 --- /dev/null +++ b/figs/fowler-nordheim.fig/asymptotic.gnuplot @@ -0,0 +1,59 @@ +datafile="asymptotic.dat" + +## can also set the following options +#set title "" +set ylabel "$|\\psi_{\\mathrm{FN}}|^2$" tc ls 1 #norotate +set y2label "$J_{\\mathrm{FN}}$" tc ls 2 #norotate +set xlabel "$x$" +# +#set xrange[:] +#set yrange [:] +set y2range [0:0.004] +# +## start ticks at 0, then every x +#set xtics 0,x +#set ytics 0,x +## puts 4 minor tics between tics (5 intervals, i.e. every 0.01) +set mxtics 5 +set mytics 5 +set my2tics 5 + +# default output canvas size: 12.5cm x 8.75cm +set term lua tikz size 12.5,8.75 standalone +# run +## gnuplot gnuplots && gnuplot_tikz out/latext/minimizer.tex + +set key off + +# 3=1+2 draw bottom and left sides of the box +#set border 3 +# don't show tics on opposite sides +set xtics nomirror +set ytics nomirror tc ls 1 +set y2tics nomirror tc ls 2 + +# Mathematica colors: +## 3f3d99 (dark blue) +## 9c4275 (dark pink) +## 9a8d3f (dark yellow) +## 3d9956 (dark green) +# My colors +## 4169E1 (pastel blue) +## DC143C (bright red) +## 32CD32 (bright green) +## 4B0082 (deep purple) +## DAA520 (ochre) + +# set linestyle +set style line 1 linetype rgbcolor "#4169E1" linewidth 3 +set style line 2 linetype rgbcolor "#DC143C" linewidth 3 +set style line 3 linetype rgbcolor "#32CD32" linewidth 3 +set style line 4 linetype rgbcolor "#4B0082" linewidth 3 +set style line 5 linetype rgbcolor "#DAA520" linewidth 3 + +set pointsize 0.6 + +set arrow to 0, graph 1 nohead lt 0 + +plot datafile using 1:2 with lines linestyle 1 ,\ + datafile using 1:3 with lines linestyle 2 axes x1y2 diff --git a/figs/fowler-nordheim.fig/asymptotic.jl b/figs/fowler-nordheim.fig/asymptotic.jl new file mode 100644 index 0000000..fd1d492 --- /dev/null +++ b/figs/fowler-nordheim.fig/asymptotic.jl @@ -0,0 +1,46 @@ +using QuadGK +using SpecialFunctions +using FFTW + +# numerical values +hbar=6.58e-16 # eV.s +m=9.11e-31 # kg +Vn=9 # eV +En=14e9 # V/m +Kn=4.5 # eV + +V=1 +E=En*hbar/(2*Vn^1.5*m^0.5)*sqrt(1.60e-19) +k0=sqrt(Kn/Vn) + +# rescale x to nm +nm_scale=sqrt(2*m*Vn)/hbar*1e9*sqrt(1.60e-19) + +# cutoffs +p_cutoff=20*k0 +p_npoints=256 + +# airy approximations +airy_threshold=30 +airy_order=5 + +# xbounds +xmax=10 +xmin=-10 +x_npoints=200 + +include("FN_base.jl") + +# compute asymptotic values +ps_asym=Array{Complex{Float64}}(undef,x_npoints) +dps_asym=Array{Complex{Float64}}(undef,x_npoints) +for i in 1:x_npoints + x=xmin+(xmax-xmin)*i/x_npoints + ps_asym[i]=psi_pole(x,k0,E,V) + dps_asym[i]=dpsi_pole(x,k0,E,V) +end + +# print asymptotic values +for i in 1:x_npoints + print((xmin+(xmax-xmin)*i/x_npoints)*nm_scale,' ',abs(ps_asym[i])^2,' ',J(ps_asym[i],dps_asym[i]),'\n') +end |