Ian Jauslin
summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--Jauslin_Hagedornfest_2019.tex288
-rw-r--r--Makefile50
-rw-r--r--README50
-rw-r--r--figs/Millikan-Lauritsen_current.pngbin0 -> 29963 bytes
-rw-r--r--figs/anim_laser/Makefile12
-rw-r--r--figs/anim_laser/animate.py103
-rw-r--r--figs/anim_laser/laser.jl391
-rw-r--r--figs/animation/FN_base.jl170
-rw-r--r--figs/animation/Makefile12
-rw-r--r--figs/animation/animate.py127
-rw-r--r--figs/animation/animate_compute.jl66
-rw-r--r--figs/contour.fig/Makefile30
-rw-r--r--figs/contour.fig/contour.py82
-rw-r--r--figs/contour.fig/contour.tikz.tex34
-rw-r--r--figs/emitter.jpgbin0 -> 14324 bytes
-rw-r--r--figs/fowler-nordheim.fig/FN_base.jl170
-rw-r--r--figs/fowler-nordheim.fig/Makefile29
-rw-r--r--figs/fowler-nordheim.fig/asymptotic.gnuplot59
-rw-r--r--figs/fowler-nordheim.fig/asymptotic.jl46
-rw-r--r--figs/libs/Makefile28
l---------figs/potential.fig/Makefile1
-rw-r--r--figs/potential.fig/potential.tikz.tex17
-rw-r--r--figs/potential.fig/potential_square.tikz.tex17
-rw-r--r--figs/potential.fig/potential_square_photonic.tikz.tex22
-rw-r--r--figs/potential.fig/potential_square_thermal.tikz.tex17
-rw-r--r--libs/ian-presentation.cls187
26 files changed, 2008 insertions, 0 deletions
diff --git a/Jauslin_Hagedornfest_2019.tex b/Jauslin_Hagedornfest_2019.tex
new file mode 100644
index 0000000..c81be67
--- /dev/null
+++ b/Jauslin_Hagedornfest_2019.tex
@@ -0,0 +1,288 @@
+\documentclass{ian-presentation}
+
+\usepackage[hidelinks]{hyperref}
+\usepackage{graphicx}
+\usepackage{array}
+
+\begin{document}
+\pagestyle{empty}
+\hbox{}\vfil
+\bf\Large
+\hfil Time-evolution of electron emission\par
+\smallskip
+\hfil from a metal surface\par
+\vfil
+\large
+\hfil Ian Jauslin
+\normalsize
+\vfil
+\hfil\rm joint with {\bf Ovidiu Costin}, {\bf Rodica Costin}, and {\bf Joel L. Lebowitz}\par
+\vfil
+arXiv:{\tt \href{http://arxiv.org/abs/1808.00936}{1808.00936}}\hfill{\tt \href{http://ian.jauslin.org}{http://ian.jauslin.org}}
+\eject
+
+\setcounter{page}1
+\pagestyle{plain}
+
+\title{Field emission}
+\vfill
+\hfil\includegraphics[height=5cm]{emitter.jpg}
+\vfill
+\eject
+
+\title{Field emission}
+$$
+ V(x)=U\Theta(x)
+ ,\quad
+ E_{\mathrm F}=k_{\mathrm F}^2<U
+$$
+\hfil\includegraphics[height=5cm]{potential_square.pdf}
+\vfill
+\eject
+
+\title{Thermal emission}
+$$
+ V(x)=U\Theta(x)
+ ,\quad
+ k^2>U
+$$
+\hfil\includegraphics[height=5cm]{potential_square_thermal.pdf}
+\vfill
+\eject
+
+\title{Photonic emission}
+$$
+ V_t(x)=\Theta(x)(U-E_tx)
+ ,\quad
+ E_t=2\epsilon\omega\cos(\omega t)
+$$
+\hfil\includegraphics[height=5cm]{potential_square_photonic.pdf}
+\vfill
+\eject
+
+\title{Field emission}
+$$
+ V(x)=\Theta(x)(U-Ex)
+$$
+\hfil\includegraphics[height=5cm]{potential.pdf}
+\vfill
+\eject
+
+\title{Field emission}
+\begin{itemize}
+ \item \href{https://doi.org/10.1073\%2Fpnas.14.1.45}{[Millikan, Lauritsen, 1928]}: experimental plot of the logarithm of the current against $1/E$
+\end{itemize}
+\hfil\includegraphics[height=4.5cm]{Millikan-Lauritsen_current.png}
+\vfill
+\eject
+
+\title{Field emission through a triangular barrier}
+\vfill
+\begin{itemize}
+ \item \href{https://doi.org/10.1098/rspa.1928.0091}{[Fowler, Nordheim, 1928]}: predicted that the current is, for small $E$,
+ $$
+ J\approx CE^2e^{-\frac aE}
+ $$
+ \item (\href{https://doi.org/10.1088/1751-8113/44/5/05530@}{[Rokhlenko, 2011]}: studied the range of applicability of the approximation, and found more accurate approximations for larger fields.)
+\end{itemize}
+\vfill
+\eject
+
+\title{Fowler-Nordheim equation}
+\begin{itemize}
+ \item Schr\"odinger equation
+ $$
+ i\partial_t\psi=-\Delta\psi+\Theta(x)(U-Ex)\psi
+ $$
+ \item Fowler-Nordheim: stationary solution: $\psi_{\mathrm{FN}}(x,t)=e^{-ik^2t}\varphi_{\mathrm{FN}}(x)$
+ $$
+ \varphi_{\mathrm{FN}}(x)=
+ \left\{ \begin{array}{l@{\ }l}
+ e^{ikx}+R_Ee^{-ikx} & x<0\\
+ T_E\mathrm{Ai}(e^{-\frac{i\pi}3}(E^{\frac13}x-E^{-\frac23}(U-k^2)) & x>0
+ \end{array}\right.
+ $$
+ $R_E$ and $T_E$ are chosen so that $\varphi_{\mathrm{FN}}$ and $\partial\varphi_{\mathrm{FN}}$ are continuous at $x=0$.
+\end{itemize}
+\vfill
+\eject
+
+\title{Fowler-Nordheim equation}
+\vfill
+\hfil\includegraphics[height=5.5cm]{asymptotic.pdf}
+\vfill
+\eject
+
+\title{Initial value problem}
+\begin{itemize}
+ \item Initial condition:
+ $$
+ \psi(x,0)=
+ \left\{ \begin{array}{l@{\ }l}
+ e^{ikx}+R_0e^{-ikx} & x<0\\
+ T_0 e^{-\sqrt{U-k^2}x} & x>0
+ \end{array}\right.
+ $$
+ $R_0$ and $T_0$ ensure that $\psi$ and $\partial\psi$ are continuous.
+ \item {\bf Theorem}: $\psi(x,t)$ behaves asymptotically like $\psi_{\mathrm{FN}}$:
+ $$
+ \psi(x,t)
+ =\psi_{\mathrm{FN}}(x,t)+\left(\frac{t}{\tau_E(x)}\right)^{-\frac32}+O(t^{-\frac52})
+ .
+ $$
+\end{itemize}
+\vfill
+\eject
+
+\title{Photoemission}
+\begin{itemize}
+ \item Photoelectric effect: early observations: \href{https://doi.org/10.1002/andp.18872670827}{[Hertz, 1887]}, \href{https://doi.org/10.1002/andp.18882690206}{[Hallwachs, 1888]}, \href{https://doi.org/10.1002/andp.19003070611}{[Lenard, 1900]}.
+ \item When a metal is irradiated with ultra-violet light, electrons are ionized, with kinetic energies at integer multiples of $\hbar\omega$.
+ \item \href{https://doi.org/10.1002/andp.19053220607}{[Einstein, 1905]}: interpretation: quanta of light ({\it photons}) of energy $\hbar\omega$ are absorbed by the electrons, whose kinetic energy is raised by $n\hbar\omega$, and can escape the metal.
+\end{itemize}
+\vfill
+\eject
+
+\title{Photoemission}
+\begin{itemize}
+ \item Time dependent potential:
+ $$
+ V_t(x)=\Theta(x)(U-2\epsilon\omega\cos(\omega t)x)
+ $$
+ \item Magnetic gauge:
+ $$
+ \Psi(x,t)
+ :=\psi(x,t)e^{-ix\Theta(x)A(t)}
+ ,\quad
+ A(t):=\int_0^t ds\ 2\epsilon\omega\cos(\omega s)
+ =
+ 2\epsilon\sin(\omega t)
+ $$
+ satisfies
+ $$
+ i\partial_t\Psi(x,t)=\left((-i\nabla+\Theta(x)A(t))^2+\Theta(x)U\right)\Psi(x,t)
+ $$
+\end{itemize}
+\vfill
+\eject
+
+\title{Periodic solution}
+\vskip-10pt
+\begin{itemize}
+ \item \href{https://doi.org/10.1103/PhysRevA.72.023412}{[Faisal, Kami\'nski, Saczuk, 2005]}
+ $$
+ \Psi_{\mathrm{FKS}}(x,t)=\left\{\begin{array}{ll}
+ e^{ikx}\exp\left(-ik^2t\right)+\Psi_R(x,t)&\mathrm{\ if\ }x<0\\
+ \Psi_T(x,t)&\mathrm{\ if\ }x>0
+ \end{array}\right.
+ $$
+ $$
+ \Psi_R(x,t)=\sum_{M\in\mathbb Z}R_Me^{-iq_Mx}\exp\left(-iq_M^2t\right)
+ ,\quad
+ q_M=\sqrt{k^2+M\omega}
+ $$
+ $$
+ \Psi_T(x,t)=\sum_{M\in\mathbb Z}T_Me^{ip_Mx}\exp\left(-iUt-i\int_0^td\tau\ (p_M-A(\tau))^2\right)
+ $$
+ $$
+ p_M=\sqrt{k^2-U+M\omega-2\epsilon^2}
+ $$
+ \item $\Psi(x,t)$, $(-i\nabla+\Theta(x)A(t))\Psi(x,t)$ are continuous.
+\end{itemize}
+\vfill
+\eject
+
+\title{Initial value problem}
+\begin{itemize}
+ \item Initial condition:
+ $$
+ \Psi(x,0)=
+ \left\{ \begin{array}{l@{\ }l}
+ e^{ikx}+R_0e^{-ikx} & x<0\\
+ T_0 e^{-\sqrt{U-k^2}x} & x>0
+ \end{array}\right.
+ $$
+ $R_0$ and $T_0$ ensure that $\Psi$ and $\partial\Psi$ are continuous.
+ \item {\bf Conjecture} (in progress): $\Psi(x,t)$ behaves asymptotically like $\Psi_{\mathrm{FKS}}$:
+ $$
+ \psi(x,t)
+ =\psi_{\mathrm{FKS}}(x,t)+\left(\frac{t}{\tau_{\mathrm{FKS}}(x)}\right)^{-\frac32}+O(t^{-\frac52})
+ .
+ $$
+\end{itemize}
+\vfill
+\eject
+
+\title{Idea of the proof: field emission}
+\begin{itemize}
+ \item Laplace transform:
+ $$
+ \hat\psi_p(x):=\int_0^\infty dt\ e^{-pt}\psi(x,t)
+ $$
+ \item Schr\"odinger equation:
+ $$
+ (-\Delta+\Theta(x)V(x)-ip)\psi_p(x)=-i\psi(x,0)
+ ,\quad
+ V(x):=U-Ex
+ $$
+\end{itemize}
+\vfill
+\eject
+
+\title{Solution in Laplace space}
+\begin{itemize}
+ \item For simplicity, $R_0\equiv T_0\equiv0$.
+ \item Solution:
+ $$
+ \hat\psi_p(x)=
+ \left\{\begin{array}{>\displaystyle l@{\ }l}
+ c(p)e^{\sqrt{-ip}x}-\frac{ie^{ikx}}{-ip+k^2}
+ &\mathrm{if\ }x<0\\[0.5cm]
+ d(p)\varphi_p(x)
+ &\mathrm{if\ }x> 0
+ \end{array}\right.
+ $$
+ with
+ $$
+ (-\Delta+V(x)-ip)\varphi_p(x)=0
+ $$
+ $$
+ \varphi_p(x)=\mathrm{Ai}\left(e^{-\frac{i\pi}3}\left(E^{\frac13}x-E^{-\frac23}(U-ip)\right)\right)
+ $$
+\end{itemize}
+\vfill
+\eject
+
+\title{Solution in Laplace space}
+\begin{itemize}
+ \item $c$ and $d$ ensure that $\hat\psi_p(x)$ and $\partial\hat\psi_p(x)$ are continuous at $x=0$:
+ $$
+ c(p)=\frac{i(ik\varphi_p(0)-\partial\varphi_p(0))}{(-ip+k^2)(\sqrt{-ip}\varphi_p(0)-\partial\varphi_p(0))}
+ $$
+ $$
+ d(p)=-\frac{i}{(\sqrt{-ip}+ik)(\sqrt{-ip}\varphi_p(0)-\partial\varphi_p(0))}.
+ $$
+\end{itemize}
+\vfill
+\eject
+
+\title{Poles in Laplace plane}
+\vfill
+\hfil\includegraphics[height=5.5cm]{contour.pdf}
+\vfill
+\eject
+
+\title{Asymptotic behavior}
+\begin{itemize}
+ \item As $t\to\infty$:
+ $$
+ \psi(x,t)
+ =\psi_{\mathrm{FN}}(x,t)+\left(\frac{t}{\tau_E(x)}\right)^{-\frac32}+O(t^{-\frac52})
+ .
+ $$
+
+ \item If $k<0$ (reflected wave), then there is no pole on the imaginary axis, so there is no contribution as $t\to\infty$.
+ \item Similarly, the transmitted wave in the initial condition does not contribute.
+\end{itemize}
+
+\end{document}
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..4ed7b60
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,50 @@
+PROJECTNAME=$(basename $(wildcard *.tex))
+LIBS=$(notdir $(wildcard libs/*))
+FIGS=$(notdir $(wildcard figs/*.fig))
+
+PDFS=$(addsuffix .pdf, $(PROJECTNAME))
+SYNCTEXS=$(addsuffix .synctex.gz, $(PROJECTNAME))
+
+all: $(PROJECTNAME)
+
+$(PROJECTNAME): $(LIBS) $(FIGS) images
+ pdflatex -file-line-error $@.tex
+ pdflatex -synctex=1 $@.tex
+
+$(SYNCTEXS): $(LIBS) $(FIGS) images
+ pdflatex -synctex=1 $(patsubst %.synctex.gz, %.tex, $@)
+
+
+$(LIBS):
+ ln -fs libs/$@ ./
+
+
+$(FIGS):
+ make -C figs/$@
+ for pdf in $$(find figs/$@/ -name '*.pdf'); do ln -fs "$$pdf" ./ ; done
+
+images:
+ for jpg in figs/*.jpg; do ln -fs "$$jpg" ./ ; done
+ for png in figs/*.png; do ln -fs "$$png" ./ ; done
+
+
+clean-aux: clean-figs-aux
+ rm -f $(addsuffix .aux, $(PROJECTNAME))
+ rm -f $(addsuffix .log, $(PROJECTNAME))
+ rm -f $(addsuffix .out, $(PROJECTNAME))
+
+clean-libs:
+ rm -f $(LIBS)
+
+clean-figs:
+ $(foreach fig,$(addprefix figs/, $(FIGS)), make -C $(fig) clean; )
+ rm -f $(notdir $(wildcard figs/*.fig/*.pdf))
+
+clean-figs-aux:
+ $(foreach fig,$(addprefix figs/, $(FIGS)), make -C $(fig) clean-aux; )
+
+
+clean-tex:
+ rm -f $(PDFS) $(SYNCTEXS)
+
+clean: clean-aux clean-tex clean-libs clean-figs
diff --git a/README b/README
new file mode 100644
index 0000000..7976122
--- /dev/null
+++ b/README
@@ -0,0 +1,50 @@
+This directory contains the source files to typeset the presentation, and
+generate the figures. This can be accomplished by running
+ make
+
+Two animations were prepared for this talk. The first shows the time-evolution
+of the density for field emission. The source files are located in the
+'figs/animation' directory. To run the animation, run
+ cd figs/animation
+ make run
+The second shows the time-evolution of the density and current for
+photoemission. The source files are located in the 'figs/anim_laser' directory.
+To run the animation, run
+ cd figs/anim_laser
+ make run
+This can take a while the first time it is run.
+
+This document uses a custom class file, located in the 'libs' directory, which
+defines a number of commands.
+
+
+* Dependencies:
+
+ pdflatex
+ TeXlive packages:
+ amsfonts
+ array
+ graphics
+ hyperref
+ latex
+ pgf
+ standalone
+ GNU make
+ julia
+ python3
+ gnuplot
+This directory contains the source files to typeset the presentation, and
+generate the figures. This can be accomplished by running
+ make
+
+* Files:
+
+ Jauslin_Hagedornfest_2019.tex:
+ main LaTeX file
+
+ libs:
+ custom LaTeX class file
+
+ figs:
+ source code for the figures
+
diff --git a/figs/Millikan-Lauritsen_current.png b/figs/Millikan-Lauritsen_current.png
new file mode 100644
index 0000000..5645be8
--- /dev/null
+++ b/figs/Millikan-Lauritsen_current.png
Binary files differ
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
new file mode 100644
index 0000000..d84ad1f
--- /dev/null
+++ b/figs/emitter.jpg
Binary files differ
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}
+
diff --git a/libs/ian-presentation.cls b/libs/ian-presentation.cls
new file mode 100644
index 0000000..91bd487
--- /dev/null
+++ b/libs/ian-presentation.cls
@@ -0,0 +1,187 @@
+%%
+%% Ian's presentation class
+%%
+
+%% TeX format
+\NeedsTeXFormat{LaTeX2e}[1995/12/01]
+
+%% class name
+\ProvidesClass{ian-presentation}[2017/09/29]
+
+\def\ian@defaultoptions{
+ \pagestyle{plain}
+ \RequirePackage{color}
+ \RequirePackage{amssymb}
+}
+
+%% paper dimensions
+\setlength\paperheight{240pt}
+\setlength\paperwidth{427pt}
+
+%% fonts
+\input{size11.clo}
+\DeclareOldFontCommand{\rm}{\normalfont\rmfamily}{\mathrm}
+\DeclareOldFontCommand{\sf}{\normalfont\sffamily}{\mathsf}
+\DeclareOldFontCommand{\tt}{\normalfont\ttfamily}{\mathtt}
+\DeclareOldFontCommand{\bf}{\normalfont\bfseries}{\mathbf}
+\DeclareOldFontCommand{\it}{\normalfont\itshape}{\mathit}
+
+%% text dimensions
+\textheight=208pt
+\textwidth=370pt
+\hoffset=-1in
+\voffset=-1in
+\oddsidemargin=24pt
+\evensidemargin=24pt
+\topmargin=8pt
+\headheight=0pt
+\headsep=0pt
+\marginparsep=0pt
+\marginparwidth=0pt
+\footskip=16pt
+
+
+%% remove default skips
+\parindent=0pt
+\parskip=0pt
+\baselineskip=0pt
+
+%% something is wrong with \thepage, redefine it
+\gdef\thepage{\the\c@page}
+
+%% correct vertical alignment at the end of a document
+\AtEndDocument{
+ % save total slide count
+ \immediate\write\@auxout{\noexpand\gdef\noexpand\slidecount{\thepage}}
+ \vfill
+ \eject
+}
+
+
+%% footer
+\def\ps@plain{
+ \def\@oddhead{}
+ \def\@evenhead{\@oddhead}
+ \def\@oddfoot{\tiny\hfill\thepage/\safe\slidecount\hfill}
+ \def\@evenfoot{\@oddfoot}
+}
+\def\ps@empty{
+ \def\@oddhead{}
+ \def\@evenhead{\@oddhead}
+ \def\@oddfoot{}
+ \def\@evenfoot{\@oddfoot}
+}
+
+
+%% title of slide
+\def\title#1{
+ \hfil{\bf\large #1}\par
+ \hfil\vrule width0.75\textwidth height0.3pt\par
+ \vskip5pt
+}
+
+
+%% hyperlinks
+% hyperlinkcounter
+\newcounter{lncount}
+% hyperref anchor
+\def\hrefanchor{%
+\stepcounter{lncount}%
+\hypertarget{ln.\thelncount}{}%
+}
+
+%% define a command and write it to aux file
+\def\outdef#1#2{%
+ % define command%
+ \expandafter\xdef\csname #1\endcsname{#2}%
+ % hyperlink number%
+ \expandafter\xdef\csname #1@hl\endcsname{\thelncount}%
+ % write command to aux%
+ \immediate\write\@auxout{\noexpand\expandafter\noexpand\gdef\noexpand\csname #1\endcsname{\csname #1\endcsname}}%
+ \immediate\write\@auxout{\noexpand\expandafter\noexpand\gdef\noexpand\csname #1@hl\endcsname{\thelncount}}%
+}
+
+%% can call commands even when they are not defined
+\def\safe#1{%
+ \ifdefined#1%
+ #1%
+ \else%
+ {\color{red}\bf?}%
+ \fi%
+}
+
+
+%% itemize
+\newlength\itemizeskip
+% left margin for items
+\setlength\itemizeskip{20pt}
+\newlength\itemizeseparator
+% space between the item symbol and the text
+\setlength\itemizeseparator{5pt}
+% penalty preceding an itemize
+\def\itemizepenalty{0}
+% counter counting the itemize level
+\newcounter{itemizecount}
+
+% item symbol
+\def\itemizept#1{
+ \ifnum#1=1
+ \textbullet
+ \else
+ $\scriptstyle\blacktriangleright$
+ \fi
+}
+
+\newlength\current@itemizeskip
+\setlength\current@itemizeskip{0pt}
+\def\itemize{
+ \par\penalty\itemizepenalty\medskip\penalty\itemizepenalty
+ \addtocounter{itemizecount}{1}
+ \addtolength\current@itemizeskip{\itemizeskip}
+ \leftskip\current@itemizeskip
+}
+\def\enditemize{
+ \addtocounter{itemizecount}{-1}
+ \addtolength\current@itemizeskip{-\itemizeskip}
+ \par\leftskip\current@itemizeskip
+ \medskip
+}
+\newlength\itempt@total
+\def\item{
+ \settowidth\itempt@total{\itemizept\theitemizecount}
+ \addtolength\itempt@total{\itemizeseparator}
+ \par
+ \medskip
+ \hskip-\itempt@total\itemizept\theitemizecount\hskip\itemizeseparator
+}
+
+%% enumerate
+\newcounter{enumerate@count}
+\def\enumerate{
+ \setcounter{enumerate@count}0
+ \let\olditem\item
+ \let\olditemizept\itemizept
+ \def\item{
+ % counter
+ \stepcounter{enumerate@count}
+ % set header
+ \def\itemizept{\theenumerate@count.}
+ % hyperref anchor
+ \hrefanchor
+ % define tag (for \label)
+ \xdef\tag{\theenumerate@count}
+ \olditem
+ }
+ \itemize
+}
+\def\endenumerate{
+ \enditemize
+ \let\item\olditem
+ \let\itemizept\olditemizept
+}
+
+
+%% end
+\ian@defaultoptions
+
+\endinput