Authors: David Cyganski, Bill Page
Views : 45
Description: From Poirier’s Bohmian Mechanics without Wavefunctions to Hall’s Many Interacting Worlds (1-D)
Compute Environment: Ubuntu 18.04 (Deprecated)

# 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.


︠e5ce94ab-7772-486b-b9b6-6c041b1c308e︠
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

$\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)
Dminus(Dminus(P,0),0)
Dplus(Dminus(Dminus(P,0),0),0)
Dplus(Dplus(Dminus(Dminus(P,0),0),0),0)

$-x\left(n - 1\right) + x\left(n\right)$
$-2 \, x\left(n - 1\right) + x\left(n - 2\right) + x\left(n\right)$
$x\left(n + 1\right) + 3 \, x\left(n - 1\right) - x\left(n - 2\right) - 3 \, x\left(n\right)$
$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(n))-(X(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])

$\left[x\left(n\right)\right]$
$\left(\begin{array}{r} -x\left(n - 1\right) + x\left(n\right) \end{array}\right)$
$\left(\begin{array}{r} -\frac{1}{x\left(n - 1\right) - x\left(n\right)} \end{array}\right)$
$\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

$\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==r(n))

$\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

$\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]$
$\left[x_{0}, x_{1}, x_{2}, x_{3}, x_{4}\right]$
eq18n3=map(lambda ex: ex.subs(dict(zip(xsk,xsv))),eq18n); eq18n3

$\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.subs({hbar:2/10,mm:1}), vars=xsv, domain=float)

p1=[random() for i in xsv];p1

$\left[0.826055631769, 0.287873495586, 0.0578043754268, 0.0243675702724, 0.674135524924\right]$
rf(*p1)

$514.24975008$
%timeit [rf(*p1)]

625 loops, best of 3: 387 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

$\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.

U(n)

$\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(eq20n==U(n))

$\mathrm{False}$
eq20n3=eq20n.subs(dict(zip(xsk,xsv))); eq20n3

$\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.788919852967$

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

eq20n4=U(n).subs(dict(zip(xsk,xsv))); eq20n4

$\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.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])


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')

$102$ 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).rjust(5), (20,1),color='black') +
list_plot(zip(X[0:N],map(PO,X[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[0:N],den(X[0:N])),xmin=-5, xmax=20,ymin=0,ymax=75,
axes_labels=['location','density'],axes_labels_size=1)+
text("time:"+("%4.1f"%X).rjust(5), (20,75),color='black')
for X in zip(tt,XV)]).show(gif=true,delay=20) # momentum density
animate([line2d(zip(mu,den(mu)),xmin=-3,xmax=3,ymin=0,ymax=10^3,
axes_labels=['velocity','density'],axes_labels_size=1)+
text("time:"+("%4.1f"%mu).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[N:2*N]))

$-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[N:2*N])-QU(XV[0:N])-PU(XV[0:N])


$-0.0114684504031$