CoCalc Shared FilesHall1Dtunneling.sagewsOpen in CoCalc with one click!
Authors: David Cyganski, Bill Page
Views : 19
Description: From Poirier’s Bohmian Mechanics without Wavefunctions to Hall’s Many Interacting Worlds (1-D)

From Poirier's Bohmian Mechanics without Wavefunctions
To Hall's Many Interacting Worlds#

Ref:

  1. Quantum Phenomena Modeled by Interactions between Many Classical Worlds
    Michael J. W. Hall Dirk-André Deckert and Howard M. Wiseman,
    PHYSICAL REVIEW X 4, 041013 (23 October 2014)
    http://arxiv.org/abs/1402.6144

  2. Quantum Mechanics Without Wavefunctions
    Jeremy Schiff and Bill Poirier
    J. Chem. Phys. 136, 031102 (2012)
    http://arxiv.org/abs/1201.2382v1

  3. Verlet integration (Wikipedia)



%typeset_mode True from numpy import array, concatenate from mpmath import erfinv

1-D

Comparison of Poirier's many-D expression to Hall's 1-D "toy" expression.

vars= ['x']; d = 1 X = map(function,vars) # world index %var n ind=[n] # position of particle in world n P=map(lambda x:x(*ind),X);P
[x(n)]\left[x\left(n\right)\right]

Difference operators

For Hall's finite number of interacting worlds we need the following difference operators instead of derivatives.

def Dminus(x,i):return(x-x.subs(ind[i]==ind[i]-1)) def Dplus(x,i):return(x.subs(ind[i]==ind[i]+1)-x)
Dminus(P[0],0) Dminus(Dminus(P[0],0),0) Dplus(Dminus(Dminus(P[0],0),0),0) Dplus(Dplus(Dminus(Dminus(P[0],0),0),0),0)
x(n1)+x(n)-x\left(n - 1\right) + x\left(n\right)
2x(n1)+x(n2)+x(n)-2 \, x\left(n - 1\right) + x\left(n - 2\right) + x\left(n\right)
x(n+1)+3x(n1)x(n2)3x(n)x\left(n + 1\right) + 3 \, x\left(n - 1\right) - x\left(n - 2\right) - 3 \, x\left(n\right)
x(n+2)4x(n+1)4x(n1)+x(n2)+6x(n)x\left(n + 2\right) - 4 \, x\left(n + 1\right) - 4 \, x\left(n - 1\right) + x\left(n - 2\right) + 6 \, x\left(n\right)

etc.

Hall's "toy" expression

eqs. (24,25)

%var hbar,mm
def sif(n): return(1/((X[0](n))-(X[0](n-1)))) def sigma(n): return((sif(n)^2)*(sif(n+1)-2*sif(n)+sif(n-1))) def r(n): return(hbar^2/4/mm*(sigma(n+1)-sigma(n))) def U(n): return(hbar^2/8/mm*(sif(n+1)-sif(n))^2)

Poirier's expression

# Ensemble of Bohmian trajectories En = map(lambda x:x(*ind),X);En # Jacobian: How trajectories change with trajectory index (i.e. as a function of the initial position). Jn = matrix(map(lambda e:map(lambda c:Dminus(e,c),[i for i in range(d)]),En));Jn Kn = Jn^(-1);Kn bool(sif(n)==Kn[0,0])
[x(n)]\left[x\left(n\right)\right]
(x(n1)+x(n))\left(\begin{array}{r} -x\left(n - 1\right) + x\left(n\right) \end{array}\right)
(1x(n1)x(n))\left(\begin{array}{r} -\frac{1}{x\left(n - 1\right) - x\left(n\right)} \end{array}\right)
True\mathrm{True}

Quantum Force

Schiff and Poirier eq. (18)

eq18n = [ sum([ sum([ sum([ sum([ (hbar^2/4/mm)*Dplus( Kn[k,i]*Kn[p,j]*Dplus( Dminus( Kn[l,j],k ),l ),p ) for p in range(d)]) for k in range(d)]) for j in range(d)]) for l in range(d)]) for i in range(d)] eq18n[0]
hbar2(1x(n+2)x(n+1)2x(n+1)x(n)1x(n1)x(n)(x(n+1)x(n))21x(n+1)x(n)+1x(n1)x(n2)+2x(n1)x(n)(x(n1)x(n))2)4mm\frac{\mathit{hbar}^{2} {\left(\frac{\frac{1}{x\left(n + 2\right) - x\left(n + 1\right)} - \frac{2}{x\left(n + 1\right) - x\left(n\right)} - \frac{1}{x\left(n - 1\right) - x\left(n\right)}}{{\left(x\left(n + 1\right) - x\left(n\right)\right)}^{2}} - \frac{\frac{1}{x\left(n + 1\right) - x\left(n\right)} + \frac{1}{x\left(n - 1\right) - x\left(n - 2\right)} + \frac{2}{x\left(n - 1\right) - x\left(n\right)}}{{\left(x\left(n - 1\right) - x\left(n\right)\right)}^{2}}\right)}}{4 \, \mathit{mm}}

Compare to Hall's MIW quantum force

bool(eq18n[0]==r(n))
True\mathrm{True}

Fast Numeric Functions

Create a numeric function that depends on the nearest neighbors: x(n-2) x(n-1) x(n) x(n+1) x(n+2)

xsk=[X[i](n+j) for i in range(d) for j in [-2,-1,0,1,2]];xsk xsv=[var(vars[i]+str(2+j)) for i in range(d) for j in [-2,-1,0,1,2]]; xsv
[x(n2),x(n1),x(n),x(n+1),x(n+2)]\left[x\left(n - 2\right), x\left(n - 1\right), x\left(n\right), x\left(n + 1\right), x\left(n + 2\right)\right]
[x0,x1,x2,x3,x4]\left[x_{0}, x_{1}, x_{2}, x_{3}, x_{4}\right]
eq18n3=map(lambda ex: ex.subs(dict(zip(xsk,xsv))),eq18n); eq18n3
[hbar2(1x0x12x1x2+1x2x3(x1x2)21x1x22x2x3+1x3x4(x2x3)2)4mm]\left[\frac{\mathit{hbar}^{2} {\left(\frac{\frac{1}{x_{0} - x_{1}} - \frac{2}{x_{1} - x_{2}} + \frac{1}{x_{2} - x_{3}}}{{\left(x_{1} - x_{2}\right)}^{2}} - \frac{\frac{1}{x_{1} - x_{2}} - \frac{2}{x_{2} - x_{3}} + \frac{1}{x_{3} - x_{4}}}{{\left(x_{2} - x_{3}\right)}^{2}}\right)}}{4 \, \mathit{mm}}\right]
rf=fast_callable(eq18n3[0].subs({hbar:2/10,mm:1}), vars=xsv, domain=float)
p1=[random() for i in xsv];p1
[0.333791178387,0.000756870815695,0.458751929544,0.949630476402,0.678418378242]\left[0.333791178387, 0.000756870815695, 0.458751929544, 0.949630476402, 0.678418378242\right]
rf(*p1)
0.02272340235780.0227234023578
%timeit [rf(*p1)]
625 loops, best of 3: 413 ns per loop

Quantum Potential###

Why does Poirier eq. (20) include a second term?

eq20n = sum([ sum([ sum([ (hbar^2/4/mm)*(Kn[k,j]*Dplus( Dplus( Kn[l,j],l ),k )+1/2*Dplus(Kn[l,j],l)*Dplus(Kn[k,j],k)) for k in range(d)]) for j in range(d)]) for l in range(d)]) eq20n
((1x(n+1)x(n)+1x(n1)x(n))22(1x(n+2)x(n+1)2x(n+1)x(n)1x(n1)x(n))x(n1)x(n))hbar28mm\frac{{\left({\left(\frac{1}{x\left(n + 1\right) - x\left(n\right)} + \frac{1}{x\left(n - 1\right) - x\left(n\right)}\right)}^{2} - \frac{2 \, {\left(\frac{1}{x\left(n + 2\right) - x\left(n + 1\right)} - \frac{2}{x\left(n + 1\right) - x\left(n\right)} - \frac{1}{x\left(n - 1\right) - x\left(n\right)}\right)}}{x\left(n - 1\right) - x\left(n\right)}\right)} \mathit{hbar}^{2}}{8 \, \mathit{mm}}

Only the first term corresponds to Hall's expression.

eq20n1 = sum([ sum([ sum([ (hbar^2/4/mm)*(1/2*Dplus(Kn[l,j],l)*Dplus(Kn[k,j],k)) for k in range(d)]) for j in range(d)]) for l in range(d)]) eq20n1
hbar2(1x(n+1)x(n)+1x(n1)x(n))28mm\frac{\mathit{hbar}^{2} {\left(\frac{1}{x\left(n + 1\right) - x\left(n\right)} + \frac{1}{x\left(n - 1\right) - x\left(n\right)}\right)}^{2}}{8 \, \mathit{mm}}
bool(eq20n1==U(n))
True\mathrm{True}
eq20n3=eq20n.subs(dict(zip(xsk,xsv))); eq20n3
((1x1x21x2x3)2+2(1x1x22x2x3+1x3x4)x1x2)hbar28mm\frac{{\left({\left(\frac{1}{x_{1} - x_{2}} - \frac{1}{x_{2} - x_{3}}\right)}^{2} + \frac{2 \, {\left(\frac{1}{x_{1} - x_{2}} - \frac{2}{x_{2} - x_{3}} + \frac{1}{x_{3} - x_{4}}\right)}}{x_{1} - x_{2}}\right)} \mathit{hbar}^{2}}{8 \, \mathit{mm}}
ru=fast_callable(eq20n3.subs({hbar:2/10,mm:1}), vars=xsv, domain=float)
ru(*p1)
0.7889198529670.788919852967

Computation of total energy below shows that Hall's expression is correct.

eq20n4=U(n).subs(dict(zip(xsk,xsv))); eq20n4
hbar2(1x1x21x2x3)28mm\frac{\mathit{hbar}^{2} {\left(\frac{1}{x_{1} - x_{2}} - \frac{1}{x_{2} - x_{3}}\right)}^{2}}{8 \, \mathit{mm}}
rU=fast_callable(eq20n4.subs({hbar:2/10,mm:1}), vars=xsv, domain=float)
rU(*p1)
3.266732526583.26673252658

Classical force

a = 0.25 c = 5 #F(x)=((x-c)/(a^3*(3.14)))*exp(-(x-c)^2/a^2)*0.7 F(x)= 9*sech(5*(x-c))^2*tanh(5*(x-c))
FF=fast_callable(F(x), vars=[x], domain=float)
plot(FF,(0,10),axes_labels=['location','force'],axes_labels_size=1)

Classical potential

#PO=sage_eval("lambda x:"+fricas(-F(x)).integrate(x).simplify().unparsed_input_form());PO(x) PO(x)=9/(10*cosh(5*x+-5*c)^2)
plot(PO,(0,10),axes_labels=['location','energy'],axes_labels_size=1)

Simulation

Simulation of a single particle in one dimension over 100 worlds. The initial spacial distribution is a Gaussian "wave packet". The "uniformizing function" is the cummulative distribution, this case the error function. We can use the inverse of the uniformizing function to create a gaussian packet from a uniform range.

Initial conditions

N=100 r0=[1.95*i/N for i in range(-N/2,(N+1)/2)] p0=map(lambda x:float(2.0*erfinv(x)),r0)
line2d(zip(r0,p0),axes_labels=['range','location'],axes_labels_size=1)
# density import numpy def den(p): return map(lambda x:1/x,p[1:]-p[0:N-1]) line2d(zip(p0[1:N],den(array(p0))),axes_labels=['location','density of worlds'],axes_labels_size=1)
# initial velocity of particle in each world v0=[1.2 for i in p0]

Acceleration

# (quantum force + classical force)/mass def rf1(X,i): global N lm2=X[i-2] if i>1 else -(10.0^10) lm1=X[i-1] if i>0 else -(10.0^6) rp1=X[i+1] if i+1<N else 10.0^6 rp2=X[i+2] if i+2<N else 10.0^10 return rf(lm2,lm1,X[i],rp1,rp2)+FF(X[i])
%md Total classical potential

Total classical potential

def PU(x): return sum([PO(x[i]) for i in range(N)])

Total quantum Potential

def ru1(X,i): global N lm2=X[i-2] if i>1 else -(10.0^10) lm1=X[i-1] if i>0 else -(10.0^6) rp1=X[i+1] if i+1<N else 10.0^6 rp2=X[i+2] if i+2<N else 10.0^10 # (quantum potential return rU(lm2,lm1,X[i],rp1,rp2) def QU(x): return sum([ru1(x,i) for i in range(N)])

Total kinetic energy

def K(x): return sum(x^2)/2

Solution

Run the integration

Velocity Verlet

Störmer–Verlet discretization with continuous integrating stepsize controller

Ref: Explicit, Time Reversible, Adaptive Step Size Control by Ernst Hairer and Gustaf Söderlind

%time t1 = 10; step = 0.1; def A(x): return array([rf1(x,i) for i in range(N)]) XX = array(p0) VV = array(v0) AX = A(XX) # minimum step size ds = 0.002 # controller def G(p,q): return -alpha*sum(map(prod,zip(p,q)))/sum(map(prod,zip(q,q)))+beta*(1-ds/dt) alpha = 30.0 beta = 0.5 # step density rho = 1; dt=ds/rho; rho = rho - G(AX,VV)*dt/2 t = 0 XV = [concatenate((XX,VV))]; tt = [t] rhos = [rho] while t<t1: # Check total energy conservation T1=K(VV)+PU(XX)+QU(XX) for i in range(int(step/dt)): rho = rho + G(AX,VV)*dt dt = ds/rho; t += dt XX += VV*dt + AX/2*dt^2 VV += AX/2*dt; AX = A(XX); VV += AX/2*dt T2=K(VV)+PU(XX)+QU(XX) if abs(T1-T2)>0.05: print "Integration failed at %s."%t,"Try a shorter time step: %s."%ds/2 break XV.append(concatenate((XX,VV))); tt += [t]; rhos.append(rho)
CPU time: 23.67 s, Wall time: 23.67 s
len(rhos) line2d(zip(tt,rhos),axes_labels=['time','rho'],legend_label='step density')
102102
show(sum([line([[tt[i],XV[i][j]] for i in range(0,len(tt))],rgbcolor = hue(float(j)/N)) for j in range(N)]),figsize=[10,7],axes_labels=['time','location'],axes_labels_size=1)
animate([plot(PO,(-5,20),axes_labels=['location','classical potential'],axes_labels_size=1, figsize=[10,4]) + text("time:"+("%4.1f"%X[0]).rjust(5), (20,1),color='black') + list_plot(zip(X[1][0:N],map(PO,X[1][0:N])), size=15, xmin=-5, xmax=20,ymax=1) for X in zip(tt,XV)]).show(gif=true,delay=20)
animate([line2d(zip(X[1][0:N],den(X[1][0:N])),xmin=-5, xmax=20,ymin=0,ymax=75, axes_labels=['location','density'],axes_labels_size=1)+ text("time:"+("%4.1f"%X[0]).rjust(5), (20,75),color='black') for X in zip(tt,XV)]).show(gif=true,delay=20)
# momentum density animate([line2d(zip(mu[1],den(mu[1])),xmin=-3,xmax=3,ymin=0,ymax=10^3, axes_labels=['velocity','density'],axes_labels_size=1)+ text("time:"+("%4.1f"%mu[0]).rjust(5), (3,1000),color='black') for mu in zip(tt,map(lambda x: x[N:2*N].sort() or x[N:2*N],XV))]).show(gif=true,delay=20)
<string>:1: RuntimeWarning: divide by zero encountered in double_scalars

Velocity

def average_velocity(t): return sum(map(abs,XV[t][N:2*N]))/N def num_left(t): return sum(map(lambda x:1 if x<0 else 0,XV[t][N:2*N])) def average_left(t): return sum(map(lambda x:-x if x<0 else 0,XV[t][N:2*N]))/(0.000001+num_left(t)) def num_right(t): return sum(map(lambda x:1 if x>=0 else 0,XV[t][N:2*N])) def average_right(t): return sum(map(lambda x: x if x>0 else 0,XV[t][N:2*N]))/(0.000001+num_right(t)) tot=line2d([[tt[i],average_velocity(i)] for i in range(len(tt))],color='blue',legend_label='average',legend_color='blue',ymax=2) left=line2d([[tt[i],average_left(i)] for i in range(len(tt))],color='red',legend_label='reflection',legend_color='red',ymax=2) right=line2d([[tt[i],average_right(i)] for i in range(len(tt))],color='green',legend_label='transmission',legend_color='green',ymax=2) show(tot+left+right,figsize=[9,5],axes_labels=['time','velocity'],axes_labels_size=1)
sum(map(abs,XV[len(tt)-1][N:2*N]))-sum(map(abs,XV[0][N:2*N]))
19.475570972-19.475570972

Final state

print "# of worlds with particles moving left: ", num_left(len(tt)-1), "average velocity: ",average_left(len(tt)-1) print "# of worlds with particles moving right: ", num_right(len(tt)-1), "average velocity: ",average_right(len(tt)-1) print "overall average: ",average_velocity(len(tt)-1)
# of worlds with particles moving left: 60 average velocity: 0.586271573569 # of worlds with particles moving right: 40 average velocity: 1.63370330985 overall average: 1.00524429028

Conservation of energy

U1=line2d([[tt[i],QU(XV[i][0:N])] for i in range(len(tt))], color='blue',legend_color='blue',legend_label="quantum potential") U2=line2d([[tt[i],PU(XV[i][0:N])] for i in range(len(tt))], color='green', legend_color='green', legend_label="classical potential") K1=line2d([[tt[i],K(XV[i][N:2*N])] for i in range(len(tt))], legend_color='red',color='red',legend_label="kinetic energy") show(K1+U1+U2,figsize=[9,5],axes_labels=['time','energy'],axes_labels_size=1)
line([[tt[i],K(XV[i][N:2*N])+PU(XV[i][0:N])+QU(XV[i][0:N])] for i in range(len(tt))], legend_label="total energy",axes_labels=['time','energy'],axes_labels_size=1)
K(XV[len(tt)-1][N:2*N])+QU(XV[len(tt)-1][0:N])+PU(XV[len(tt)-1][0:N])-K(XV[1][N:2*N])-QU(XV[1][0:N])-PU(XV[1][0:N])
0.0114684504031-0.0114684504031