diff options
Diffstat (limited to 'figs')
22 files changed, 1433 insertions, 0 deletions
diff --git a/figs/Millikan-Lauritsen_current.png b/figs/Millikan-Lauritsen_current.png Binary files differnew file mode 100644 index 0000000..5645be8 --- /dev/null +++ b/figs/Millikan-Lauritsen_current.png diff --git a/figs/anim_laser/Makefile b/figs/anim_laser/Makefile new file mode 100644 index 0000000..9b66e37 --- /dev/null +++ b/figs/anim_laser/Makefile @@ -0,0 +1,12 @@ +PROJECTNAME=animate + +all: animate.dat + +run: animate.dat + python3 animate.py animate.dat & + +animate.dat: + julia laser.jl > animate.dat + +clean: + rm -f animate.dat diff --git a/figs/anim_laser/animate.py b/figs/anim_laser/animate.py new file mode 100644 index 0000000..1e6995a --- /dev/null +++ b/figs/anim_laser/animate.py @@ -0,0 +1,103 @@ +from matplotlib import pyplot as pl +from matplotlib import animation +import sys + +# read data +# time dependent data +frames=[] +## asymptotic data (located in the first block) +#asym=[] +infile=open(sys.argv[1],'r') +row=[] +for line in infile: + if line=='\n': + frames.append(row) + row=[] + else: + dat=[] + for n in line.split(): + dat.append(float(n)) + row.append(dat) +infile.close() + + +# set up plot +fig = pl.figure() +pl.subplot(211) +axr=fig.gca() +asym_rho, = axr.plot([],[],linewidth=3.5,color='#00FF00') +cn_rho, = axr.plot([],[],color='#FF0000') + +pl.subplot(212) +axJ=fig.gca() +asym_J, = axJ.plot([],[],linewidth=3.5,color='#00FF00') +cn_J, = axJ.plot([],[],color='#FF0000') + +# plot ranges +xmax=0 +maxyr=0 +maxyJ=0 +for frame in frames: + for i in range(len(frame)): + if frame[i][1]>xmax: + xmax=frame[i][1] + if frame[i][2]>maxyr: + maxyr=frame[i][2] + if frame[i][3]>maxyr: + maxyr=frame[i][3] + if frame[i][4]>maxyJ: + maxyJ=frame[i][4] + if frame[i][5]>maxyJ: + maxyJ=frame[i][5] +xmin=0 +minyr=0 +minyJ=0 +for frame in frames: + for i in range(len(frame)): + if frame[i][1]<xmin: + xmin=frame[i][1] + if frame[i][2]<minyr: + minyr=frame[i][2] + if frame[i][3]<minyr: + minyr=frame[i][3] + if frame[i][4]<minyJ: + minyJ=frame[i][4] + if frame[i][5]<minyJ: + minyJ=frame[i][5] + +# animate +def init_plot(): + axr.set_ylim(minyr,maxyr) + axr.set_xlim(xmin,xmax) + axJ.set_ylim(minyJ,maxyJ) + axJ.set_xlim(xmin,xmax) + + axr.vlines(0,minyr,maxyr,linestyles="dotted") + axJ.vlines(0,minyJ,maxyJ,linestyles="dotted") + axJ.hlines(0,xmin,xmax,linestyles="dotted") +def update(frame): + axr.set_title("t=% .3f fs" % (frame[0][0])) + xdata=[] + ydata=[] + asym_data=[] + cn_data=[] + for i in range(len(frame)): + xdata.append(frame[i][1]) + asym_data.append(frame[i][2]) + cn_data.append(frame[i][3]) + asym_rho.set_data(xdata,asym_data) + cn_rho.set_data(xdata,cn_data) + + xdata=[] + ydata=[] + asym_data=[] + cn_data=[] + for i in range(len(frame)): + xdata.append(frame[i][1]) + asym_data.append(frame[i][4]) + cn_data.append(frame[i][5]) + asym_J.set_data(xdata,asym_data) + cn_J.set_data(xdata,cn_data) +anim = animation.FuncAnimation(fig, update, frames=frames, blit=False, interval=100, repeat=True, init_func=init_plot) +#anim.save('laser_schrodinger.m4v', fps=10, extra_args=['-vcodec', 'libx264']) +pl.show() diff --git a/figs/anim_laser/laser.jl b/figs/anim_laser/laser.jl new file mode 100644 index 0000000..78cabff --- /dev/null +++ b/figs/anim_laser/laser.jl @@ -0,0 +1,391 @@ +using FFTW +using Printf +using SpecialFunctions +using LinearAlgebra + +# numerical values +hbar=6.58e-16 # eV.s +m=9.11e-31 # kg +V_=6.6 # eV +I=1e17 # W/m^2 +EF=3.3 # eV +omega_=1.5498*4 # eV + +# dimensionless quantities +V=1 +k=sqrt(EF/V_) +epsilon=sqrt(hbar^3*I/(8*m*V_*omega_^2)) +omega=omega_/V_ + +#epsilon=8*epsilon +#omega=4*omega + +nu=60 +ell=20 + +# for Crank-Nicolson +xmin_cn=-100 +xmax_cn=100 +nx_cn=4000 +tmax=10 +# mu=dt/dx^2 +mu=0.1 +nt=floor(Int,tmax/(mu*((xmax_cn-xmin_cn)/nx_cn)^2)) +# only save up to memt time steps +memt=10000 + + +# asymptotic value of \varphi_t(x) +function phi_asym(x,T,gkz,k,epsilon,omega,V,nu,ell) + # compute residues + res=residue(x,T,gkz,k,epsilon,omega,V,nu) + times=compute_times(omega,nu,ell) + N=(2*nu+1)*(2*ell+1) + fun=zeros(Complex{Float64},N) + for m in 1:N + for n in -nu:nu + fun[m]+=-1im*res[n+nu+1]*exp(-1im*(k*k+n*omega)*times[m]) + end + end + return fun +end +# asymptotic value of \partial\varphi_t(x) +function dphi_asym(x,T,gkz,k,epsilon,omega,V,nu,ell) + # compute residues + res=dresidue(x,T,gkz,k,epsilon,omega,V,nu) + times=compute_times(omega,nu,ell) + N=(2*nu+1)*(2*ell+1) + fun=zeros(Complex{Float64},N) + for m in 1:N + for n in -nu:nu + fun[m]+=-1im*res[n+nu+1]*exp(-1im*(k*k+n*omega)*times[m]) + end + end + return fun +end + +# return times +function compute_times(omega,nu,ell) + N=(2*nu+1)*(2*ell+1) + times=zeros(N) + for j in 1:N + times[j]=pi/(omega*(nu+0.5))*(j-1) + end + return times +end + +# \kappa_m^{(\sigma)} +function K(m,sigma,V,k,epsilon,omega) + return pow(V+2*epsilon*epsilon-k*k-sigma-m*omega,0.5,pi/2) +end +# mathfrak q_k +function Q(k,V) + return sqrt(V-k*k) +end + +# T_0 +function T0(k,V) + return 2im*k/(1im*k-sqrt(V-k*k)) +end + +# compute g +function g(kappa,epsilon,omega,N) + B=zeros(Complex{Float64},N) + # prepare vector of B's + for i in 1:N + t=2*pi/(omega*N)*(i-1) + B[i]=exp(1im*epsilon*epsilon/omega*sin(2*omega*t)-kappa*4*epsilon/omega*cos(omega*t)) + end + # take the fft + fft!(B) + # correct by a factor 1/N + for i in 1:N + B[i]=B[i]/N + end + return(B) +end + +# index of m-n +function mmn(n,nu) + return (n>=0 ? n+1 : n+4*nu+2) +end + +# compute residues at -ik^2-in\omega +function residue(x,T,gkz,k,epsilon,omega,V,nu) + # compute residue + out=zeros(Complex{Float64},2*nu+1) + if x>=0 + for n in -nu:nu + for m in -nu:nu + out[n+nu+1]+=1im*gkz[m+nu+1][mmn(n-m,nu)]*exp(-K(m,0,V,k,epsilon,omega)*x)*T[m+nu+1] + end + end + else + for n in -nu:nu + out[n+nu+1]=(n==0 ? 1im*(exp(1im*k*x)-exp(-1im*k*x)) : 0) + for m in -nu:nu + out[n+nu+1]+=1im*T[m+nu+1]*gkz[m+nu+1][mmn(n-m,nu)]*exp(-1im*x*pow(k*k+n*omega,0.5,-pi/2)) + end + end + end + + return out +end +# compute \partial residues at -ik^2-in\omega +function dresidue(x,T,gkz,k,epsilon,omega,V,nu) + # compute residue + out=zeros(Complex{Float64},2*nu+1) + if x>=0 + for n in -nu:nu + for m in -nu:nu + out[n+nu+1]+=-K(m,0,V,k,epsilon,omega)*1im*gkz[m+nu+1][mmn(n-m,nu)]*exp(-K(m,0,V,k,epsilon,omega)*x)*T[m+nu+1] + end + end + else + for n in -nu:nu + out[n+nu+1]=(n==0 ? -k*(exp(1im*k*x)+exp(-1im*k*x)) : 0) + for m in -nu:nu + out[n+nu+1]+=pow(k*k+n*omega,0.5,-pi/2)*T[m+nu+1]*gkz[m+nu+1][mmn(n-m,nu)]*exp(-1im*x*pow(k*k+n*omega,0.5,-pi/2)) + end + end + end + + return out +end + + +# matching condition for residue +function residue_matching(k,epsilon,omega,V,nu) + # preliminary: g^{(\kappa_m^{(0)})} + # g^{(\kappa_m^{(0)})}_n=gkz[m+nu+1][n+1] if n>=0 and gkz[m+nu+1][n+4*nu+2] if n<0 + gkz=Array{Array{Complex{Float64}},1}(undef,2*nu+1) + for m in -nu:nu + gkz[m+nu+1]=g(K(m,0,V,k,epsilon,omega),epsilon,omega,4*nu+1) + end + + # solve matching condition + G=zeros(Complex{Float64},2*nu+1,2*nu+1) + v=zeros(Complex{Float64},2*nu+1) + for n in -nu:nu + for m in -nu:nu + G[n+nu+1,m+nu+1]=((K(m,0,V,k,epsilon,omega)-1im*pow(k*k+n*omega,0.5,-pi/2))*gkz[m+nu+1][mmn(n-m,nu)]-epsilon*(gkz[m+nu+1][mmn(n-m+1,nu)]-gkz[m+nu+1][mmn(n-m-1,nu)])) + end + end + v[nu+1]=-2im*k + + return (G\v,gkz) +end + + +# 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 + +# current +function J(phi,dphi,x,epsilon,omega,t) + if x>=0 + return(2*imag(conj(phi)*(dphi+2im*epsilon*sin(omega*t)*phi))) + else + return(2*imag(conj(phi)*dphi)) + end +end + +# compute \psi_t(x) using Crank-Nicolson +function psi_cn(V,omega,epsilon,k,xs,ts,memt) + nx=length(xs) + nt=length(ts) + + # length of output vector + len=min(nt,memt) + # how often to write to output + freq=ceil(Int,nt/memt) + + # initialize + psis=zeros(Complex{Float64},nx,len) + times=zeros(Float64,len) + + # previous time step + psi_prev=zeros(Complex{Float64},nx) + + # init + for i in 1:nx + times[1]=0 + psis[i,1]=(xs[i]<0 ? exp(1im*k*xs[i])+(T0(k,V)-1)*exp(-1im*k*xs[i]) : T0(k,V)*exp(-Q(k,V)*xs[i])) + psi_prev[i]=psis[i,1] + end + # the boundary condition is that psi at the boundary is constant (nothing needs to be done for this) + + # matrix structure of the Crank-Nicolson algorithm + # diagonals of M + M0=zeros(Complex{Float64},nx-2) + Mp=zeros(Complex{Float64},nx-3) + Mm=zeros(Complex{Float64},nx-3) + v=zeros(Complex{Float64},nx-2) + + for j in 1:nt-1 + # print progress + progress(j,nt-1,1000) + + # t_j + t=ts[j] + # t_{j+1} + tp=ts[j+1] + dt=tp-t + for mu in 2:nx-1 + # x_mu + x=xs[mu] + dx2=(xs[mu+1]-xs[mu])*(xs[mu]-xs[mu-1]) + # tridiagonal matrix + M0[mu-1]=1im/dt-1/dx2-(x<0 ? 0 : (V-2*epsilon*omega*cos(omega*tp)*x)/2) + if mu<nx-1 + Mp[mu-1]=1/(2*dx2) + end + if mu>2 + Mm[mu-2]=1/(2*dx2) + end + # right side + v[mu-1]=1im*psi_prev[mu]/dt-(psi_prev[mu+1]+psi_prev[mu-1]-2*psi_prev[mu])/(2*dx2)+(x<0 ? 0 : (V-2*epsilon*omega*cos(omega*t)*x)*psi_prev[mu]/2) + # correct for boundary conditions + # assumes the boundary condition is constant in time! + if mu==2 + v[mu-1]-=psi_prev[1]/(2*dx2) + end + if mu==nx-1 + v[mu-1]-=psi_prev[nx]/(2*dx2) + end + end + + M=Tridiagonal(Mm,M0,Mp) + + # copy to psi_prev + psi_prev[2:nx-1]=M\v + + # write to output + if j%freq==0 + psis[:,Int(j/freq)+1]=psi_prev + times[Int(j/freq)+1]=tp + end + end + + return (times,psis) +end + +# returns the times in b that are closest to those in a +# assumes the vectors are ordered +function nearest_times(a,b) + out=zeros(Int,length(a)) + pointer=1 + for i in 1:length(a) + if pointer==length(b) + out[i]=length(b) + end + for j in pointer:length(b)-1 + if b[j+1]>a[i] + out[i]=(abs(b[j+1]-a[i])<abs(b[j]-a[i]) ? j+1 : j) + pointer=out[i] + break + end + # all remaining points are beyond the range + if j==length(b)-1 + out[i]=length(b) + pointer=out[i] + end + end + end + return out +end + +# space derivative using Crank Nicolson +function deriv_cn(phi,i,xs) + if i==1 + return (phi[i+1]-phi[i])/(xs[i+1]-xs[i]) + end + if i==length(phi) + return (phi[i]-phi[i-1])/(xs[i]-xs[i-1]) + end + return (phi[i+1]-phi[i])/(xs[i+1]-xs[i])/2 + (phi[i]-phi[i-1])/(xs[i]-xs[i-1])/2 +end + +################################################ + + +# print progress +function progress(j,tot,freq) + if (j-1)%ceil(Int,tot/freq)==0 + if j>1 + @printf(stderr,"\r") + end + @printf(stderr,"%d/%d",j,tot) + end + if j==tot + @printf(stderr,"%d/%d\n",j,tot) + end +end + +# print animation data using Crank Nicolson and compare to Laplace +function print_anim_cn_cmp() + @printf(stderr,"epsilon=% .8e omega=% .8e k=% .8e V=% .8e z=% .8e T=% .8e gamma=% .8e\n",epsilon,omega,k,V,4*sqrt(nu)*epsilon/sqrt(omega),2*pi/(omega_/hbar)*1e15,sqrt((V-k^2)/(4*epsilon^2))) + xmin=-10 + xmax=10 + nx=200 + + # array of times for the asymptote + times=compute_times(omega,nu,ell) + + # positions at which to compute phi + xs=Array(0:nx_cn)*(xmax_cn-xmin_cn)/nx_cn.+xmin_cn + # times at which to compute phi + ts=Array(0:nt-1)*tmax/nt + + # compute matching condition + (T,gkz)=residue_matching(k,epsilon,omega,V,nu) + + # compute wave function + (t_psi,psi)=psi_cn(V,omega,epsilon,k,xs,ts,memt) + + # nearest times using Crank Nicolson + indices_t=nearest_times(times,t_psi) + indices_x=nearest_times(Array(0:nx-1)*(xmax-xmin)/nx.+xmin,xs) + + ps_cn=Array{Array{Complex{Float64},1}}(undef,nx) + ps_asym=Array{Array{Complex{Float64},1}}(undef,nx) + dps_cn=Array{Array{Complex{Float64},1}}(undef,nx) + dps_asym=Array{Array{Complex{Float64},1}}(undef,nx) + for i in 1:nx + progress(i,nx,1000) + x=xmin+(xmax-xmin)*(i-1)/nx + ps_cn[i]=Array{Complex{Float64}}(undef,length(times)) + dps_cn[i]=Array{Complex{Float64}}(undef,length(times)) + for j in 1:length(times) + # phase change in magnetic gauge + phase=(x<0 ? 1 : exp(-2im*epsilon*sin(omega*times[j])*x)) + ps_cn[i][j]=phase*psi[indices_x[i],indices_t[j]] + dps_cn[i][j]=phase*deriv_cn(psi[:,indices_t[j]],indices_x[i],xs)-(xs[indices_x[i]]<0 ? 0 : 2im*epsilon*sin(omega*times[j])*psi[indices_x[i],indices_t[j]]) + end + ps_asym[i]=phi_asym(x,T,gkz,k,epsilon,omega,V,nu,ell) + dps_asym[i]=dphi_asym(x,T,gkz,k,epsilon,omega,V,nu,ell) + end + + # print values at each time + for j in 1:length(times) + if times[j]>tmax + break + end + for i in 1:nx + x=(xmin+(xmax-xmin)*i/nx) + # remove point at 0 + if i!=Int(nx/2) + @printf("% .8e % .8e % .8e % .8e % .8e % 8e\n",times[j]*hbar/V_*1e15,x*sqrt(hbar^2*1.6e-19/(2*m*V_))*1e9,abs(ps_asym[i][j])^2,abs(ps_cn[i][j])^2,J(ps_asym[i][j],dps_asym[i][j],x,epsilon,omega,times[j]),J(ps_cn[i][j],dps_cn[i][j],x,epsilon,omega,times[j])) + end + end + print('\n') + end +end + + +print_anim_cn_cmp() diff --git a/figs/animation/FN_base.jl b/figs/animation/FN_base.jl new file mode 100644 index 0000000..af2a1ee --- /dev/null +++ b/figs/animation/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/animation/Makefile b/figs/animation/Makefile new file mode 100644 index 0000000..d582a40 --- /dev/null +++ b/figs/animation/Makefile @@ -0,0 +1,12 @@ +PROJECTNAME=animate + +all: animate.dat + +run: animate.dat + python3 animate.py animate.dat & + +animate.dat: + julia animate_compute.jl > animate.dat + +clean: + rm -f animate.dat diff --git a/figs/animation/animate.py b/figs/animation/animate.py new file mode 100644 index 0000000..aaa3a80 --- /dev/null +++ b/figs/animation/animate.py @@ -0,0 +1,127 @@ +from matplotlib import pyplot as pl +from matplotlib import animation +import sys + +# read data +# time dependent data +frames=[] +# asymptotic data (located in the first block) +asym=[] +infile=open(sys.argv[1],'r') +row=[] +for line in infile: + # read first block + if len(asym)==0: + if line=='\n': + asym=row + row=[] + else: + dat=[] + for n in line.split(): + dat.append(float(n)) + row.append(dat) + # read other blocks + else: + if line=='\n': + frames.append(row) + row=[] + else: + dat=[] + for n in line.split(): + dat.append(float(n)) + row.append(dat) +infile.close() + + +# set up plot +fig = pl.figure() +#pl.subplot(211) +axr=fig.gca() +asym_rho, = axr.plot([],[],linewidth=3.5,color='#00FF00') +rho, = axr.plot([],[],color='red') + +#pl.subplot(212) +#axJ=fig.gca() +#asym_J, = axJ.plot([],[],linewidth=3.5,color='#00FF00') +#J, = axJ.plot([],[],color='red') + +# plot ranges +xmax=0 +maxyr=0 +maxyJ=0 +for frame in frames: + for i in range(len(frame)): + if frame[i][1]>xmax: + xmax=frame[i][1] + if frame[i][2]>maxyr: + maxyr=frame[i][2] + if frame[i][3]>maxyJ: + maxyJ=frame[i][3] +for i in range(len(asym)): + if asym[i][0]>xmax: + xmax=asym[i][0] + if asym[i][1]>maxyr: + maxyr=asym[i][1] + if asym[i][2]>maxyJ: + maxyJ=asym[i][2] +xmin=0 +minyr=0 +minyJ=0 +for frame in frames: + for i in range(len(frame)): + if frame[i][1]<xmin: + xmin=frame[i][1] + if frame[i][2]<minyr: + minyr=frame[i][2] + if frame[i][3]<minyJ: + minyJ=frame[i][3] +for i in range(len(asym)): + if asym[i][0]<xmin: + xmin=asym[i][0] + if asym[i][1]<minyr: + minyr=asym[i][1] + if asym[i][2]<minyJ: + minyJ=asym[i][2] + + +# plot asymptotes +asym_rho_datax=[] +asym_rho_datay=[] +for i in range(len(asym)): + asym_rho_datax.append(asym[i][0]) + asym_rho_datay.append(asym[i][1]) +asym_rho.set_data(asym_rho_datax,asym_rho_datay) +#asym_J_datax=[] +#asym_J_datay=[] +#for i in range(len(asym)): +# asym_J_datax.append(asym[i][0]) +# asym_J_datay.append(asym[i][2]) +#asym_J.set_data(asym_J_datax,asym_J_datay) + +# animate +def init_plot(): + axr.set_ylim(minyr,maxyr) + axr.set_xlim(xmin,xmax) + #axJ.set_ylim(minyJ,maxyJ) + #axJ.set_xlim(xmin,xmax) + + axr.vlines(0,minyr,maxyr,linestyles="dotted") + #axJ.vlines(0,minyJ,maxyJ,linestyles="dotted") +def update(frame): + axr.set_title("t=% .3f fs" % (frame[0][0])) + xdata=[] + ydata=[] + for i in range(len(frame)): + xdata.append(frame[i][1]) + ydata.append(frame[i][2]) + rho.set_data(xdata,ydata) + + #xdata=[] + #ydata=[] + #for i in range(len(frame)): + # xdata.append(frame[i][1]) + # ydata.append(frame[i][3]) + #J.set_data(xdata,ydata) +anim = animation.FuncAnimation(fig, update, frames=frames, blit=False, interval=100, repeat=True, init_func=init_plot) +#anim.save('schrodinger_barrier.mp4', fps=15, extra_args=['-vcodec', 'libx264']) +pl.show() diff --git a/figs/animation/animate_compute.jl b/figs/animation/animate_compute.jl new file mode 100644 index 0000000..452320e --- /dev/null +++ b/figs/animation/animate_compute.jl @@ -0,0 +1,66 @@ +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 +En=10e9 # 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=hbar*sqrt(1.6e-19)/sqrt(2*m*Vn)*1e9 + +# 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 wave function +ps=Array{Array{Array{Complex{Float64},1},1}}(undef,x_npoints) +dps=Array{Array{Array{Complex{Float64},1},1}}(undef,x_npoints) +for i in 1:x_npoints + print(stderr,i,'/',x_npoints,'\n') + x=xmin+(xmax-xmin)*i/x_npoints + ps[i]=psi(x,k0,E,V,p_npoints,p_cutoff) + dps[i]=dpsi(x,k0,E,V,p_npoints,p_cutoff) +end + +# 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])/(2*k0),'\n') +end +print('\n') + +# print values at each time +for j in 1:p_npoints + for i in 1:x_npoints + print(real(ps[i][1][j])*hbar/Vn*1e15,' ',(xmin+(xmax-xmin)*i/x_npoints)*nm_scale,' ',abs(ps[i][2][j])^2,' ',J(ps[i][2][j],dps[i][2][j])/(2*k0),'\n') + end + print('\n') +end diff --git a/figs/contour.fig/Makefile b/figs/contour.fig/Makefile new file mode 100644 index 0000000..bade7d0 --- /dev/null +++ b/figs/contour.fig/Makefile @@ -0,0 +1,30 @@ +PROJECTNAME=contour + +PDFS=$(addsuffix .pdf, $(PROJECTNAME)) + +all: $(PDFS) + +$(PDFS): poles.tikz.tex + pdflatex -jobname $(basename $@) -file-line-error $(patsubst %.pdf, %.tikz.tex, $@) + +poles.tikz.tex: + python3 contour.py > poles.tikz.tex + +install: $(PDFS) + cp $^ $(INSTALLDIR)/ + +$(LIBS): + ln -fs libs/$@ ./ + +clean-libs: + rm -f $(LIBS) + +clean-aux: + rm -f $(addsuffix .aux, $(PROJECTNAME)) + rm -f $(addsuffix .log, $(PROJECTNAME)) + rm -f poles.tikz.tex + +clean-tex: + rm -f $(PDFS) + +clean: clean-libs clean-aux clean-tex diff --git a/figs/contour.fig/contour.py b/figs/contour.fig/contour.py new file mode 100644 index 0000000..102df86 --- /dev/null +++ b/figs/contour.fig/contour.py @@ -0,0 +1,82 @@ +#!/usr/bin/env python3 + +import cmath +import math +import scipy.special as sp +from scipy import optimize as op +import random +import sys + +# number of roots +nr_roots=4 + +# size of plot +plotsize_x=3 +plotsize_y=3 +# rescale plot (so roots are not too close together) +plotsize_scale_x=1 +plotsize_scale_y=6 + +# numerical values +hbar=6.58e-16 # eV.s +m=9.11e-31 # kg +Vn=9 # eV +En=20e9 # V/m + +V=1 +E=En*hbar/(Vn**1.5*m**0.5)*math.sqrt(1.60e-19) + +# sqrt with branch cut along iR_+ +def sqrt_p(x): + r,p=cmath.polar(x) + if(p<cmath.pi/2): + return(cmath.rect(math.sqrt(r),p/2)) + else: + return(cmath.rect(math.sqrt(r),(p-2*math.pi)/2)) + +# solution of (-\Delta+V-ip)phi=0 +def phi(p,x,E,V): + return(sp.airy(cmath.exp(-1j*math.pi/3)*(E**(1/3)*x-E**(-2/3)*(V-1j*p)))[0]) +# its derivative +def dphi(p,x,E,V): + return(cmath.exp(-1j*math.pi/3)*E**(1/3)*sp.airy(cmath.exp(-1j*math.pi/3)*(E**(1/3)*x-E**(-2/3)*(V-1j*p)))[1]) + + +# the function whose roots are to be computed and its derivative +def f(p): + return(sqrt_p(-1j*p)*phi(p,0,E,V)-dphi(p,0,E,V)) +def df(p): + return(-1j/(2*sqrt_p(-1j*p))*phi(p,0,E,V)+sqrt_p(-1j*p)*dp_phi(p,0,E,V)-dp_dphi(p,0,E,V)) + +# derivatives of phi with respect to p +def dp_phi(p,x,E,V): + return(1j*cmath.exp(-1j*math.pi/3)*E**(-2/3)*sp.airy(cmath.exp(-1j*math.pi/3)*(E**(1/3)*x-E**(-2/3)*(V-1j*p)))[1]) +def dp_dphi(p,x,E,V): + return(-1j*(x-(V-1j*p)/E)*phi(p,x,E,V)) + +# check whether the root was already found +def check(root,roots): + # reject if the root doesn't fit on the plot + if(plotsize_scale_x*root.real<-plotsize_x or abs(plotsize_scale_y*root.imag)>plotsize_y): + return(False) + for x in roots: + if(abs(root-x)<1e-12): + return(False) + return(True) + +# find roots +roots=[] +while len(roots)<nr_roots: + try: + root=op.newton(f, -random.random()*plotsize_x/plotsize_scale_x+(2*random.random()-1)*plotsize_y/plotsize_scale_y*1j, fprime=df) + if(check(root,roots)): + roots.append(root) + print("found "+str(len(roots))+" roots of "+str(nr_roots), file=sys.stderr) + except RuntimeError: + root=0 + except: + break + +# print result +for root in roots: + print("\\pole{(% .3f,% .3f)}" % (plotsize_scale_x*root.real,plotsize_scale_y*root.imag)) diff --git a/figs/contour.fig/contour.tikz.tex b/figs/contour.fig/contour.tikz.tex new file mode 100644 index 0000000..953d8f0 --- /dev/null +++ b/figs/contour.fig/contour.tikz.tex @@ -0,0 +1,34 @@ +\documentclass{standalone} + +\usepackage{tikz} + +\def\pole#1{ + \draw#1++(0.1,0.1)--++(-0.2,-0.2); + \draw#1++(-0.1,0.1)--++(0.2,-0.2); + \draw[color=red, line width=1pt]#1circle(0.3); +} + +\begin{document} +\begin{tikzpicture} + +% branch cut +\fill[color=gray](-3.3,-0.1)--++(3.3,0)--++(0,0.2)--++(-3.3,0)--cycle; + +% axes +\draw(-3.3,0)--++(6.6,0); +\draw(0,-3.3)--++(0,6.6); + +% -ik0^2 +\pole{(0,-2)} + +% other poles +\input{poles.tikz} + +% contour +\draw[color=red, line width=1pt](-3.3,-0.3)--++(3.3,0); +\draw[color=red, line width=1pt](0,-0.3)..controls(0.3,-0.3)and(0.3,0.3)..(0,0.3); +\draw[color=red, line width=1pt](-3.3,0.3)--++(3.3,0); + + +\end{tikzpicture} +\end{document} diff --git a/figs/emitter.jpg b/figs/emitter.jpg Binary files differnew file mode 100644 index 0000000..d84ad1f --- /dev/null +++ b/figs/emitter.jpg 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 diff --git a/figs/libs/Makefile b/figs/libs/Makefile new file mode 100644 index 0000000..33b81e2 --- /dev/null +++ b/figs/libs/Makefile @@ -0,0 +1,28 @@ +PROJECTNAME=$(basename $(basename $(wildcard *.tikz.tex))) +LIBS=$(notdir $(wildcard libs/*)) + +PDFS=$(addsuffix .pdf, $(PROJECTNAME)) + +all: $(PDFS) + +$(PDFS): $(LIBS) + echo $(LIBS) + pdflatex -jobname $(basename $@) -file-line-error $(patsubst %.pdf, %.tikz.tex, $@) + +install: $(PDFS) + cp $^ $(INSTALLDIR)/ + +$(LIBS): + ln -fs libs/$@ ./ + +clean-libs: + rm -f $(LIBS) + +clean-aux: + rm -f $(addsuffix .aux, $(PROJECTNAME)) + rm -f $(addsuffix .log, $(PROJECTNAME)) + +clean-tex: + rm -f $(PDFS) + +clean: clean-libs clean-aux clean-tex diff --git a/figs/potential.fig/Makefile b/figs/potential.fig/Makefile new file mode 120000 index 0000000..704310e --- /dev/null +++ b/figs/potential.fig/Makefile @@ -0,0 +1 @@ +../libs/Makefile
\ No newline at end of file diff --git a/figs/potential.fig/potential.tikz.tex b/figs/potential.fig/potential.tikz.tex new file mode 100644 index 0000000..ac1ba36 --- /dev/null +++ b/figs/potential.fig/potential.tikz.tex @@ -0,0 +1,17 @@ +\documentclass{standalone} + +\usepackage{tikz} + +\begin{document} +\begin{tikzpicture} + +\draw(-3,0)--(3,0); +\draw(0,-1.5)--(0,3); + +\draw[color=gray,line width=1pt, densely dotted](-2.5,1.25)--++(2.5,0); + +\draw[line width=1.5pt](-2.5,0)--(0,0)--(0,2.5)--(1.975,-1.25); + +\end{tikzpicture} +\end{document} + diff --git a/figs/potential.fig/potential_square.tikz.tex b/figs/potential.fig/potential_square.tikz.tex new file mode 100644 index 0000000..bd13467 --- /dev/null +++ b/figs/potential.fig/potential_square.tikz.tex @@ -0,0 +1,17 @@ +\documentclass{standalone} + +\usepackage{tikz} + +\begin{document} +\begin{tikzpicture} + +\draw(-3,0)--(3,0); +\draw(0,-1.5)--(0,3); + +\draw[color=gray,line width=1pt, densely dotted](-2.5,1.25)--++(2.5,0); + +\draw[line width=1.5pt](-2.5,0)--(0,0)--(0,2.5)--(2.5,2.5); + +\end{tikzpicture} +\end{document} + diff --git a/figs/potential.fig/potential_square_photonic.tikz.tex b/figs/potential.fig/potential_square_photonic.tikz.tex new file mode 100644 index 0000000..8a35fb4 --- /dev/null +++ b/figs/potential.fig/potential_square_photonic.tikz.tex @@ -0,0 +1,22 @@ +\documentclass{standalone} + +\usepackage{tikz} + +\begin{document} +\begin{tikzpicture} + +\draw(-3,0)--(3,0); +\draw(0,-1.5)--(0,3); + +\draw[color=gray,line width=1pt, densely dotted](-2.5,1.25)--++(2.5,0); +\draw[color=gray,line width=1pt, densely dotted](-1,2.00)--++(1,0); +\draw[color=gray,line width=1pt, densely dotted](-1,2.75)--++(3.5,0); + +\draw[->,line width=0.75pt,color=red](-0.5,1.25)--++(0,0.75); +\draw[->,line width=0.75pt,color=red](-0.5,2.00)--++(0,0.75); + +\draw[line width=1.5pt](-2.5,0)--(0,0)--(0,2.5)--(2.5,2.5); + +\end{tikzpicture} +\end{document} + diff --git a/figs/potential.fig/potential_square_thermal.tikz.tex b/figs/potential.fig/potential_square_thermal.tikz.tex new file mode 100644 index 0000000..3e8fd22 --- /dev/null +++ b/figs/potential.fig/potential_square_thermal.tikz.tex @@ -0,0 +1,17 @@ +\documentclass{standalone} + +\usepackage{tikz} + +\begin{document} +\begin{tikzpicture} + +\draw(-3,0)--(3,0); +\draw(0,-1.5)--(0,3); + +\draw[color=gray,line width=1pt, densely dotted](-2.5,2.75)--++(5,0); + +\draw[line width=1.5pt](-2.5,0)--(0,0)--(0,2.5)--(2.5,2.5); + +\end{tikzpicture} +\end{document} + |
