1
\documentclass[notitlepage]{article}
2
\usepackage{titling}
3
\usepackage[utf8]{inputenc}
4
\usepackage{titling}
5
\usepackage{fullpage}
6
\usepackage{mathtools}
7
\usepackage{sagetex}
8
\usepackage{url}
9
%\usepackage{amsmath}
10
\usepackage{cite}
11
\usepackage{hyperref}
12
\usepackage{xcolor}
13
\hypersetup{
14
    colorlinks,
15
    linkcolor={red!50!black},
16
    citecolor={blue!50!black},
17
    urlcolor={blue!80!black}
18
}
19
\usepackage{breqn}
20
\usepackage{collect}
21
\setlength{\parskip}{1ex}
22

23
\title{Discretization of Poirier's Quantum Dynamical Equations of Motion}
24
\author{Bill Page \small{(\href{mailto:bill.page@newsynthesis.org}{bill.page@newsynthesis.org})}}
25
\date{\today}
26

27
\begin{document}
28

29
\maketitle
30
\thispagestyle{empty}
31

32
\begin{abstract}
33
\textbf{From Poirier’s Bohmian Mechanics without Wavefunctions, to Hall’s Many Interacting Worlds in More Than One Dimension}
34

35
In this note we show that Hall's "toy" model of quantum mechanics as many interacting interacting classical worlds one dimension follows directly from the approximation of Poirier's quantum mechanical equations of motion as difference equations with a judicious choice of difference operators.  Using the same approach we can construct similar models for the two- and three-dimensional case. We observe that the discretization of Poirier's expression for the quantum force is identical to the gradient of the discretization of the quantum potential in both the 1-D case and in the 2 and 3-D cases when using a symmetric difference operator.
36
\end{abstract}
37

38
\setcounter{section}{-1}
39
\section{Preliminary}
40
\begin{sagesilent}
41
import numpy
42
import mpmath
43
\end{sagesilent}
44
In this document we have used the computer algebra system
45
Sage\footnote{\sagestr{version()}\cite{SageMath2016} on SageMathCloud\cite{SageMathCloud2016}} with components of
46
SciPy\footnote{numpy Version \sagestr{numpy.version.version} \cite{SciPy2016}}
47
and
48
mpmath\footnote{mpmath Version \sagestr{mpmath.__version__} \cite{mpmath2016}}
49
for both symbolic and numeric calculations. The source code for this document can be found
50
online\footnote{\url{https://cloud.sagemath.com/projects/b04b5777-e269-4c8f-a4b8-b21dbe1c93c6/files/discretization.tex}}. The inset bolded text below and many of the formulas in this paper are input and output produced directly by the Sage system.
51
\begin{sageblock}
52
from numpy import array,concatenate,isnan
53
from mpmath import erfinv
54
\end{sageblock}
55
It is convenient to define a compact notation for functions and some special symbols:
56

57
\definecollection{prelim}
58
\begin{sagesilent}
59
begincol=chr(92)+"begin { collect* } { prelim }"
60
endcol=chr(92)+"end { collect* }"
61
\end{sagesilent}
62

63
ok { \sagestr{begincol} }
64
\begin{sageblock}
65
def fun(*y): # fun(x,y)(f)(x,y)=f
66
    cache={}
67
    def upd(f,x,y): # memoized
68
        try: return cache[x]
69
        except KeyError:
70
            try: cache[x] = f.subs(dict([[y[i],x[i]] for i in range(len(y))]))
71
            except AttributeError: # ask forgiveness
72
                def sub(f): return f.subs(dict([[y[i],x[i]] for i in range(len(y))]))
73
                cache[x] = map(sub,f)
74
        return cache[x]
75
    return lambda f: lambda *x: upd(f,x,y)
76
def argscript(self, *args): return "%s_{%s}"%(self.name(),','.join(map(repr, args)))
77
def functions(*names):
78
    return map(lambda nam:function(nam, print_latex_func=argscript),names)
79
def detJscript(self, *args):
80
    return "\\begin{vmatrix}\\scriptstyle{J}\\end{vmatrix}_{%s}"% (
81
        ','.join(map(repr, args)))
82
detJ = function('detJ', print_latex_func=detJscript)
83
def squared(x): return x.dot(x)  # vector (numpy array) x^2
84
hbar = var('hbar',latex_name='\\hbar')
85
mu = var('mu',latex_name='\\mu')
86
\end{sageblock}
87
\sagestr{endcol}
88

89

90
The following is a work-round for a bug\footnote{Some commands for manipulating and simplifying expressions do not respect the \textbf{\rm print\_latex\_func} option. \url{https://trac.sagemath.org/ticket/19151}} in Sage.
91
\begin{sageblock}
92
def fix(it):
93
    return it.subs(dict([[function('x')(n+i,m+j),X[0](n+i,m+j)]
94
                            for i in range(-1,2) for j in range(-1,2)] +
95
                        [[function('y')(n+i,m+j),X[1](n+i,m+j)]
96
                            for i in range(-1,2) for j in range(-1,2)]))
97
\end{sageblock}
98

99
\section{2-Dimensional System}
100
In Schiff and Poirier \cite{schiff2012communication} eq. (18), the quantum dynamical equation of motion
101

102
\begin{equation}
103
\mu \ddot{x^i}+\frac{\partial V(x)}{\partial x^i}- \frac{\hbar^2}{4 m}\frac{\partial}{\partial C^m}\left(K^k_i K^m_j \frac{\partial^2 K^l_j}{\partial C^k \partial C^l}\right)=0
104
\end{equation}
105
when specialized to 2-D and discretized by replacing derivatives with difference operators is the 2-D analog of Hall’s \cite{hall2014quantum} “toy” expression eq. (24).
106

107
By discretized we mean that the "uniformizing co-ordinate" $\mathbf{C}$ will be replaced with discrete values representing a finite number of "worlds". Derivatives with respect to $\mathbf{C}$ will be replaced with difference operators in the expressions for the quantum potential and quantum force.
108

109
Let $E =[x_{n,m},y_{n,m}]$ represent the configuration space of one particle in 2-D space or two particles in 1-D space in the world with co-ordinates $n,m$. In Sage this can be represented as:
110
\begin{sageblock}
111
X = functions('x','y')
112
d = len(X)
113
C = var('n,m') # C=(n,m)
114
E = tuple([x(*C) for x in X])
115
\end{sageblock}
116
\begin{equation} \sage{E} \end{equation}
117

118
Forward and backward difference approximations to the derivatives are
119
\begin{sageblock}
120
def DM(x,n):return x-x.subs(n==n-1)
121
def DP(x,n):return x.subs(n==n+1)-x
122
\end{sageblock}
123
Similarly the central difference approximation to this derivative is
124
\begin{sageblock}
125
def D(x,n):return x.subs(n==n+1)/2-x.subs(n==n-1)/2
126
\end{sageblock}
127

128
To derive Hall's equations for a finite number of interacting worlds we replace derivatives with difference operators so the Jacobian becomes:
129
\begin{sageblock}
130
J=fun(n,m)(matrix([[DM(E[j],C[i]) for i in range(d)] for j in range(d)]))
131
\end{sageblock}
132
\begin{equation} \sage{J(n,m)} \end{equation}
133
It will be convenient for our calculations to represent the inverse of the Jacobian in terms of its adjugate and delay evaluation of the determinant until needed. The adjugate matrix is
134
\begin{sageblock}
135
A=fun(n,m)(fix(J(n,m).adjoint()))
136
\end{sageblock}
137
\begin{equation} \sage{A(n,m)} \end{equation}
138
The inverse is
139
\begin{sageblock}
140
K=fun(n,m)(A(n,m)/detJ(n,m))
141
\end{sageblock}
142
\begin{equation} \sage{K(n,m)} \end{equation}
143
The determinant of the Jacobian is
144
\begin{sageblock}
145
Jdet=fun(n,m)(fix(J(n,m).determinant().expand()))
146
dets = {detJ(n+i,m+j):Jdet(n+i,m+j) for i in range(-2,3) for j in range(-2,3)}
147
\end{sageblock}
148
%\begin{equation} \sage{Jdet(n,m)} \end{equation}
149

150
It is $\sage{fix(factor(J(n,m)^(-1)))==K(n,m).subs(detJ(n,m)==Jdet(n,m))}$ that $J(n,m)^{-1}=K(n,m) \text{ where }{\sage{detJ(n,m)}=\operatorname{Jdet}(n,m)}$.
151

152
\subsection{Quantum Force}
153

154
The Schiff and Poirier expression eq. (18) in 2-D has many more terms than in the 1-D case!
155

156
\begin{sageblock}
157
R=fun(n,m) (vector([ (hbar^2/4/mu)*sum( sum( sum( sum(
158
 DP( K(n,m)[k,i] * K(n,m)[p,j] * DP( DM( K(n,m)[l,j],C[l] ),C[k] ),C[p] )
159
     for p in range(d)) for k in range(d)) for j in range(d)) for l in range(d))
160
         for i in range(d)])
161
\end{sageblock}
162

163
\subsection{Quantum Potential}
164

165
\begin{sageblock}
166
U=fun(n,m)(
167
  sum( sum( sum(
168
    -(hbar^2/4/mu) *
169
    (K(n,m)[k,j] * DM( DP( K(n,m)[l,j],C[l] ),C[k] ) +
170
      1/2 * DP(K(n,m)[l,j],D[l]) * DP(K(n,m)[k,j],C[k]))
171
        for k in range(d)) for j in range(d)) for l in range(d)) )
172
\end{sageblock}
173

174
\begin{dmath}
175
\sage{ eq20.subs({hbar:1,mu:1}) }
176
\end{dmath}
177

178
\begin{sageblock}
179
U2=fun(n,m)(
180
  sum( sum( sum( -(hbar^2/4/mu)*(1/2 * DP(K(n,m)[l,j],C[l]) * DP(K(n,m)[k,j],C[k]))
181
    for k in range(d)) for j in range(d)) for l in range(d)) )
182
\end{sageblock}
183

184
\begin{dmath}
185
\sage{ eq20n2 }
186
\end{dmath}
187

188
\begin{sageblock}
189
U3=fun(n,m)(
190
  sum( sum( sum( -(hbar^2/4/mu)*(K(n,m)[k,j] * DM( DP( K(n,m)[l,j],C[l] ),C[k] ))
191
    for k in range(d)) for j in range(d)) for l in range(d)) )
192
\end{sageblock}
193

194
\begin{dmath}
195
\sage{ eq20n3 }
196
\end{dmath}
197

198
It is $\sage{bool(eq20 == eq20n2+eq20n3)}$ that $eq20 = eq20n2+eq20n3$.
199

200
In the discretized expressions it will be convenient to replace the $X$ co-ordinates with scalar quantities. The forward and backward difference approximations will then require $2x5x5 = 50$ scalar quantities.
201

202
\begin{sageblock}
203
xyk = [X[i](n+j,m+k)
204
          for j in range(-2,3) for k in range(-2,3) for i in range(d)]
205
xyv = [var(X[i].name()+str(2+j)+str(2+k))
206
          for j in range(-2,3) for k in range(-2,3) for i in range(d)]
207
\end{sageblock}
208

209
\begin{sageblock}
210
xnm = X[0](n,m).subs(dict(zip(xyk,xyv)))
211
ynm = X[1](n,m).subs(dict(zip(xyk,xyv)))
212
\end{sageblock}
213
\begin{dmath}
214
\sage{[xnm,ynm]}
215
\end{dmath}
216
\bibliographystyle{bibstyles/utphys}
217
\bibliography{bibliography}
218
\end{document}
219

220
We can compute the discretized quantum force as the gradient of the discretized quantum potential symbolically. We would like to compare this with the discretization of Poirier's expression for the force.
221

222
Let $U(n,m)$ be the potential at $n,m$ with $U2(n,m)$ and $U3(n,m)$ being just the 2nd and 3rd order terms respectively.
223
\begin{sageblock}
224
def U(n,m): return eq20.subs({var('n'):n,var('m'):m})
225
def U2(n,m): return eq20n2.subs({var('n'):n,var('m'):m})
226
def U3(n,m): return eq20n3.subs({var('n'):n,var('m'):m})
227
\end{sageblock}
228
Gradients:
229
\begin{sageblock}
230
dets = {detJ(n+i,m+j):Jdet(n+i,m+j) for i in range(-2,3) for j in range(-1,3)}
231
rrx = -sum(U(n+i,m+j) for i in range(-2,3) for j in range(-2,3)
232
          ).subs(dets).subs(dict(zip(xyk,xyv))).diff(xnm)
233
rry = -sum([U(n+i,m+j) for i in range(-2,3) for j in range(-2,3)]
234
          ).subs(dets).subs(dict(zip(xyk,xyv))).diff(ynm)
235
\end{sageblock}
236
Or computing the gradients of the 2nd and 3rd order terms separately:
237
\begin{sageblock}
238
rr2x = -sum([U2(n+i,m+j) for i in range(-2,3) for j in range(-2,3)]
239
           ).subs(dets).subs(dict(zip(xyk,xyv))).diff(xnm)
240
rr2y = -sum([U2(n+i,m+j) for i in range(-2,3) for j in range(-2,3)]
241
           ).subs(dets).subs(dict(zip(xyk,xyv))).diff(ynm)
242
rr3x = -sum([U3(n+i,m+j) for i in range(-2,3) for j in range(-2,3)]
243
           ).subs(dets).subs(dict(zip(xyk,xyv))).diff(xnm)
244
rr3y = -sum([U3(n+i,m+j) for i in range(-2,3) for j in range(-2,3)]
245
           ).subs(dets).subs(dict(zip(xyk,xyv))).diff(ynm)
246
\end{sageblock}
247

248
It is $\sage{bool(rr2x.expand()+rr3x.expand() == rrx.expand())}$ that $rrx = rr2x+rr3x$.
249

250
It is $\sage{bool(rr2y.expand()+rr3y.expand() == rry.expand())}$ that $rrx = rr2y+rr3y$.
251

252
Surprisingly, just like in the 1-D case, the force term from the 3rd order term is exactly -2 times the force term obtained from the 2nd order term so the difference between the force from just the 2nd order term so the force from the full expression for the quantum potential is differs from the 2nd term by only a sign.
253

254
It is $\sage{bool(rr3x.expand() == -2*rr2x.expand())}$ that $rr3x = -2 rr2x$.
255

256
It is $ \sage{bool(rr3y.expand() == -2*rr2y.expand())}$ that $rr3y = -2 rr2y$.
257

258
But unlike the 1-D case, the result for the full expression is not the same as the expression we obtained directly by substituting difference operators for derivatives in Poirier's expression (eq18n) for the quantum force!
259

260
\begin{sageblock}
261
test=dict(zip(xyv,[random() for i in xyv]))
262
\end{sageblock}
263

264
It is $\sage{bool(eq18[0].subs(dets).subs(dict(zip(xyk,xyv))).subs(test)==rrx.subs(test))}$ that $eq18[0] = rrx$.
265

266
It is $\sage{bool(eq18[1].subs(dets).subs(dict(zip(xyk,xyv))).subs(test)==rry.subs(test))}$ that $eq18[1] = rry$.
267

268
Repeating all of the above using a symmetric difference approximation to the derivatives we find however that the two expressions for the quantum force turn out to be the same.
269

270
\begin{sageblock}
271
@CachedFunction
272
def J(n,m):
273
    return matrix([[D(E[j],i).subs({var('n'):n,var('m'):m})
274
             for i in range(d)]
275
               for j in range(d)])
276
@CachedFunction
277
def A(n,m): return fix(J(n,m).adjoint()).subs({var('n'):n,var('m'):m})
278
@CachedFunction
279
def K(n,m): return (A(n,m)/detJ(n,m)).subs({var('n'):n,var('m'):m})
280
@CachedFunction
281
def Jdet(n,m): return J(n,m).determinant().expand()
282
\end{sageblock}
283
\begin{equation} \sage{J(n,m)} \end{equation}
284
\begin{equation} \sage{A(n,m)} \end{equation}
285
\begin{equation} \sage{K(n,m)} \end{equation}
286
\begin{sageblock}
287
eq18 = [ (hbar^2/4/mu)*sum( sum( sum( sum(
288
 D( K(n,m)[k,i] * K(n,m)[p,j] * D( D( K(n,m)[l,j],l ),k ),p )
289
   for p in range(d))
290
     for k in range(d))
291
       for j in range(d))
292
         for l in range(d)) for i in range(d)]
293
eq20 = sum( sum( sum(
294
 -(hbar^2/4/mu) *
295
   (K(n,m)[k,j] * D( D( K(n,m)[l,j],l ),k ) +
296
     1/2 * D(K(n,m)[l,j],l) * D(K(n,m)[k,j],k))
297
       for k in range(d)) for j in range(d)) for l in range(d))
298
eq20n2 = sum( sum( sum( -(hbar^2/4/mu)*(1/2 * DP(K(n,m)[l,j],l) * DP(K(n,m)[k,j],k))
299
   for k in range(d)) for j in range(d)) for l in range(d))
300
eq20n3 = sum( sum( sum( -(hbar^2/4/mu)*(K(n,m)[k,j] * DM( DP( K(n,m)[l,j],l ),k ))
301
   for k in range(d)) for j in range(d)) for l in range(d))
302
\end{sageblock}
303

304
\begin{dmath} \sage{eq20} \end{dmath}
305

306
\begin{sageblock}
307
xyk = [X[i](n+j,m+k)
308
          for j in range(-4,5) for k in range(-4,5) for i in range(d)]
309
xyv = [var(X[i].name()+str(5+j)+str(5+k))
310
          for j in range(-4,5) for k in range(-4,5) for i in range(d)]
311
xnm = X[0](n,m).subs(dict(zip(xyk,xyv)))
312
ynm = X[1](n,m).subs(dict(zip(xyk,xyv)))
313
\end{sageblock}
314

315
\begin{dmath}
316
\sage{[xnm,ynm]}
317
\end{dmath}
318

319
\begin{sageblock}
320
def U(n,m): return eq20.subs({var('n'):n,var('m'):m})
321
def U2(n,m): return eq20n2.subs({var('n'):n,var('m'):m})
322
def U3(n,m): return eq20n3.subs({var('n'):n,var('m'):m})
323
dets = {detJ(n+i,m+j):Jdet(n+i,m+j) for i in range(-6,7) for j in range(-6,7)}
324
rrx = -sum(U(n+i,m+j) for i in range(-3,4) for j in range(-3,4)
325
          ).subs(dets).subs(dict(zip(xyk,xyv))).diff(xnm)
326
rry = -sum([U(n+i,m+j) for i in range(-3,4) for j in range(-3,4)]
327
          ).subs(dets).subs(dict(zip(xyk,xyv))).diff(ynm)
328
rr2x = -sum([U2(n+i,m+j) for i in range(-3,4) for j in range(-3,4)]
329
           ).subs(dets).subs(dict(zip(xyk,xyv))).diff(xnm)
330
rr2y = -sum([U2(n+i,m+j) for i in range(-3,4) for j in range(-3,4)]
331
           ).subs(dets).subs(dict(zip(xyk,xyv))).diff(ynm)
332
rr3x = -sum([U3(n+i,m+j) for i in range(-3,4) for j in range(-3,4)]
333
           ).subs(dets).subs(dict(zip(xyk,xyv))).diff(xnm)
334
rr3y = -sum([U3(n+i,m+j) for i in range(-3,4) for j in range(-3,4)]
335
           ).subs(dets).subs(dict(zip(xyk,xyv))).diff(ynm)
336

337
\end{sageblock}
338

339
\begin{sageblock}
340
vals=[RDF(random()) for i in xyv]
341
test=dict(zip(xyv,vals))
342
constants = {hbar:0.2,mu:1.0}
343
\end{sageblock}
344

345

346
It is $\sage{squared(map(lambda x:x.subs(test).subs(constants),array([rr2x,rr2y]) - array([rrx,rry])))<1e-15}$ that $[rrx,rry] = [rr2x,rr2y]+[rr3x,rr3y]$.
347

348
The force term from the 3rd order term is still exactly -2 times the force term obtained from the 2nd order term.
349

350
It is $\sage{squared(map(lambda x:x.subs(test).subs(constants),array([rr3x,rr3y]) - -2*array([rr2x,rr2y])))<1e-15}$ that $[rr3x,rr3y] = -2 [rr2x,rr2y]$.
351

352
But now the full expression is the same as the expression we obtained directly by substituting difference operators for derivatives in Poirier's expression (eq18n) for the quantum force. It is very time consuming to show this symbolically but numerical serve as a practical demonstration.
353

354
\begin{sageblock}
355
ff0=fast_callable(rrx.subs(constants), vars=xyv, domain=RDF)
356
ff1=fast_callable(rry.subs(constants), vars=xyv, domain=RDF)
357
def ff(p): return array([ff0(*p),ff1(*p)])
358
gg0=fast_callable(eq18[0].subs(dets).subs(dict(zip(xyk,xyv))).subs(constants), vars=xyv, domain=RDF)
359
gg1=fast_callable(eq18[1].subs(dets).subs(dict(zip(xyk,xyv))).subs(constants), vars=xyv, domain=RDF)
360
def gg(p): return array([gg0(*p),gg1(*p)])
361
\end{sageblock}
362

363
It is $\sage{ abs(ff0(*vals)-gg0(*vals))  < 1.0e-12 and abs(ff1(*vals)-gg1(*vals)) < 1.0e-12 }$ that $eq18 = [rrx,rry]$.
364

365
\section{Simulation}
366

367
Using the expressions derived above we can generate efficient numerical functions for use in the solution of the equations of motion.
368

369
Classical Potential
370
\begin{sageblock}
371
    P0(x)=-9/(2*cosh(2*x)^2)
372
\end{sageblock}
373
\sageplot[width=.5\textwidth]{
374
	plot(P0,(-3,3),axes_labels=['distance','energy'],
375
	axes_labels_size=1,legend_label='potential')
376
}
377

378
Classical Force
379
\begin{sageblock}
380
	F0(x)= -(18*sinh(2*x))/(cosh(2*x)^3)
381
\end{sageblock}
382
\sageplot[width=.5\textwidth]{
383
    plot(F0,(-3,3),axes_labels=['distance','force'],axes_labels_size=1)
384
}
385

386
Quantum Potential
387

388
$ \sage{U(n,m).subs(dets).subs(dict(zip(xyk,xyv))).subs(constants).variables()} $
389

390
%\begin{sageblock}
391
%ru0=fast_callable(U(n,m).subs(dets).subs(dict(zip(xyk,xyv))).subs(constants), 
392
%vars=xyv, domain=RDF)
393
%def ru(x): return ru0(*x)
394
%\end{sageblock}
395

396
\bibliographystyle{bibstyles/utphys}
397
\bibliography{bibliography}
398
\end{document}
399
%sagemathcloud={"latex_command":"pdflatex -synctex=1 -interact=nonstopmode 'discretization.tex'"}
400

401