Ian Jauslin
summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorIan Jauslin <jauslin@ias.edu>2019-05-22 22:59:03 -0400
committerIan Jauslin <jauslin@ias.edu>2019-05-22 22:59:03 -0400
commitdf7449e4a29ec8d3182cf7b2aebcb86f7ac596c2 (patch)
tree182c5c85f57264b38c60a02ec268e5c00d937b34 /figs/fowler-nordheim.fig/FN_base.jl
As presented at VirginiaTech on 2019-05-24HEADv1.0master
Diffstat (limited to 'figs/fowler-nordheim.fig/FN_base.jl')
-rw-r--r--figs/fowler-nordheim.fig/FN_base.jl170
1 files changed, 170 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