1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
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
|
