1
\documentclass[notitlepage]{article}
2
\usepackage{titling}
3
\usepackage[utf8]{inputenc}
4
\usepackage{titling}
5
\usepackage{fullpage}
6
\usepackage{mathtools}
7
\usepackage{etoolbox}
8
\cslet{lrbox*}\lrbox
9
\expandafter\patchcmd\csname lrbox*\endcsname{\setbox}{\global\setbox}{}{}
10
%\expandafter\show\csname lrbox*\endcsname % uncomment to see if it has worked
11
\cslet{endlrbox*}\endlrbox
12
\usepackage{sagetex}
13
\usepackage{url}
14
%\usepackage{amsmath}
15
\usepackage{cite}
16
\usepackage{xcolor}
17
\usepackage{breqn}
18
\usepackage{collect}
19
\setlength{\parskip}{1ex}
20
\usepackage[titletoc,title]{appendix}
21
\usepackage{float}
22
\usepackage{hyperref}
23
\hypersetup{
24
    colorlinks,
25
    linkcolor={red!50!black},
26
    citecolor={blue!50!black},
27
    urlcolor={blue!80!black}
28
}
29

30
\title{Discretization of the Quantum Dynamical Equations of Motion - One Dimension}
31
\author{Bill Page \small{(\href{mailto:bill.page@newsynthesis.org}{bill.page@newsynthesis.org})}}
32
\date{Draft \today}
33

34
\begin{document}
35

36
\maketitle
37
\thispagestyle{empty}
38

39
\begin{abstract}
40
\textbf{From Poirier’s Bohmian Mechanics without Wavefunctions, to Hall’s Many Interacting Worlds
41
in One Dimension}
42

43
In this note we show that Hall's ``toy'' model of quantum mechanics as many interacting interacting classical worlds in 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. These calculations are a preliminary exercise for the case of more than one dimension. The main emphasis is methodological with the main applications to follow in later articles.
44

45
We provide a numerical simulation that reproduces the results given previously by Poirier\cite{poirier2008}.
46
\end{abstract}
47

48
%\setcounter{section}{-1}
49

50
\section{Introduction}
51

52
In a general sense the quantum mechanical equations of motion originates with the work of Erwin Madelung\cite{madelung1927} on quantum hydrodynamics. David Bohm\cite{bohm1993} rediscovered and extended an interpretation of quantum mechanics originating with Louis de Broglie\cite{debroglie1960}, aka. Bohmian mechanics, which the notion of particles ``guided'' by the wavefunction. Peter Holland\cite{holland1995} and co-workers have continued to extend the de Broglie-Bohm theory in this direction, ultimately showing that it is possible to dispense with the wavefunction all together in favor of a continuous ensemble of quantum trajectories. Besides the conceptual metaphysical issues, this has lead to some new more powerful calculation methods that have become important for practical application of quantum mechanics to physical chemistry (e.g. Wyatt \cite{wyatt2005}). Poirier continues in this tradition while Hall, Deckert, and Wiseman\cite{hall2014quantum} emphasize the interpretation of the quantum trajectories as interaction between parallel worlds reminiscent but differing from Hugh Everett's many-worlds interpretation of quantum mechanics.
53

54
In this article the Hamiltonian formulaion of these equations of motion is important. Conservation of energy is an important indicator of the quality of the numerical integration. The Störmer–Verlet method implemented in this article is well known and preserves the symplectic structure on the system phase space.
55

56
From the point of view of Poirier\cite{schiff2012communication} this simulation is viewed as an approximation to the quantum dynamical field equations. One should think of this in hydrodynamic terms. But from the point of view of Hall, et al. these are discrete trajectories of the \textbf{same} particle in many different ``parallel'' worlds - but just one particle!
57

58
In the expressions which follow $\mathbf{C}$ is the co-ordinate which indexes these worlds. One should think of this as the same particle in many different worlds interacting with the copy of itself in other worlds. It is a an important aspect of the equations of motion that trajectories of these particles never cross. The co-ordinate function $\mathbf{C}$ is assumed to be ``uniformized''\footnote{ Schiff and Poirier\cite{schiff2012communication}} so that discretization of $\mathbf{C}$  assigns an equal ``density'' to the particle in each world.
59

60
In the example simulation included in this article the starting distribution of the particle in these different worlds is a Gaussian distribution. But the distribution evolves as it interacts with the classical potential. Some trajectories are scattered from the potential while other tunnel through. From the point of view of Poirier, an important goal of this simulation is to provide an efficient calculation of the scattering and tunneling amplitudes. But from the point of view of Hall this is a just purely "classical" mechanical problem with a rather unusual form of interaction. It is the form of this interaction between worlds that leads to statistics that approximate quantum mechanics in the limit of a very large number of interacting worlds.
61

62
\newcommand{\J}{\begin{vmatrix}\scriptstyle{J}\end{vmatrix}}
63

64
\definecollection{prelim}
65
\begin{sagesilent}[prelim]
66
import numpy
67
from numpy import array,concatenate,vectorize,isnan
68
rdf=numpy.array([RDF.pi()]).dtype # ensure compatible types
69
import mpmath
70
from mpmath import erfinv
71

72
def fun(*y): # fun(x,y)(f)(x,y)=f
73
    cache={}
74
    def upd(f,x,y): # memoized
75
        try: return cache[x]
76
        except KeyError: # not found
77
            try: cache[x] = f.subs(dict([[y[i],x[i]]
78
                                        for i in range(len(y))]))
79
            except AttributeError: # ask forgiveness
80
                def sub(f): return f.subs(dict([[y[i],x[i]]
81
                                        for i in range(len(y))]))
82
                cache[x] = map(sub,f)
83
        return cache[x]
84
    return lambda f: lambda *x: upd(f,x,y)
85
def argscript(self, *args):
86
    return "%s_{%s}"%(self.name(),','.join(map(repr, args)))
87
def functions(*names):
88
    return map(lambda nam:function(nam, print_latex_func=argscript),names)
89
def detJscript(self, *args):
90
    return '\J_{%s}'% (','.join(map(repr, args)))
91
J__= function('J_') # determinant of J
92
J_ = function('J_', print_latex_func=detJscript)
93
def squared(x): # x^2
94
    try: return x.dot(x) # numpy array
95
    except: return power(x,2)
96
\end{sagesilent}
97

98
\definecollection{shift}
99
\begin{sagesilent}[shift]
100
def shift(x,n):
101
    try: # deduce slicing
102
        start = (x.__array_interface__['data'][0] -
103
                 x.base.__array_interface__['data'][0])/x.itemsize
104
        stop = start + x.shape[0]
105
        if start+n<0 or stop+n>x.base.shape[0]:
106
            raise RuntimeError('bounds')
107
        return x.base[start+n:stop+n]
108
    except AttributeError:
109
        raise RuntimeError('expecting view')
110
\end{sagesilent}
111

112

113
In this document we have used the computer algebra system
114
Sage\footnote{\sagestr{version()}\cite{SageMath2016} on SageMathCloud\cite{SageMathCloud2016}} with components of
115
SciPy\footnote{numpy Version \sagestr{numpy.version.version} \cite{SciPy2016}}
116
and
117
mpmath\footnote{mpmath Version \sagestr{mpmath.__version__} \cite{mpmath2016}}
118
for both symbolic and numeric calculations. This document and the source code for this work can be found
119
online\footnote{\LaTeX\  \href{https://cloud.sagemath.com/projects/b04b5777-e269-4c8f-a4b8-b21dbe1c93c6/files/miw/\jobname.tex}{source} \href{https://cloud.sagemath.com/projects/b04b5777-e269-4c8f-a4b8-b21dbe1c93c6/files/miw/\jobname.pdf}{pdf} \href{https://cloud.sagemath.com/projects/b04b5777-e269-4c8f-a4b8-b21dbe1c93c6/files/miw/\jobname.sagews}{SageWorksheet}. The online versions may be more up to date.}. The inset bold text below and many of the formulas in this paper are input and output produced directly by the Sage system. For more details please refer to Appendix \ref{SageAppendix}.
120

121
For example we define some special symbols for Planck's constant and mass:
122
\begin{sageblock}
123
hbar = var('hbar',latex_name=r'\hbar')
124
m = var('m')
125
\end{sageblock}
126

127
\definecollection{fancy}
128
\definecollection{prelim1}
129
\begin{collect}{prelim1}{\noindent}{\par}
130
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.
131
\end{collect}
132
\definecollection{prelim2}
133
\begin{sagecollect}{prelim2}
134
def fix(it):
135
    ret = it.subs(dict([[function('x')(n+i), X[0](n+i)] for i in range(-2,3)] +
136
                       [[J__(n+i), J_(n+i)] for i in range(-2,3)]))
137
    return ret
138
def factor(x):
139
    try: return fix(x.factor())
140
    except: return vector(map(lambda y:fix(y.factor()),x))
141
\end{sagecollect}
142

143
\section{One Dimensional System}
144

145
We show below that in one dimension Schiff and Poirier's quantum dynamical equation of motion,  \cite{schiff2012communication} eqs. (18-20),
146
\begin{equation} \label{eq18}
147
m \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
148
\end{equation}
149
that follows from a quantum potential of the form
150
\begin{equation} \label{eq20}
151
Q = -\frac{\hbar^2}{4 m} \left(K^k_j \frac{\partial^2 K^l_j}{\partial C^l \partial C^k} + \frac{1}{2}\frac{\partial K^l_j}{\partial C^l}\frac{\partial K^k_j}{\partial C^k} \right)
152
\end{equation}
153

154
when discretized by replacing derivatives with difference operators, is equivalent to Hall’s \cite{hall2014quantum} ``toy'' expression eq. (24) for a finite number of classically interacting ``parallel'' worlds.
155

156
\begin{align}
157
r_N(x_n;X) &= \frac{\hbar^2}{4 m}[\sigma_{n+1}(X)-\sigma_n(X)] \label{eq24} \\
158
\sigma_n(X) &= \frac{1}{(x_n-x_{n-1})^2} \left[ \frac{1}{x_{n+1}-x_n} - \frac{2}{x_n-x_{n-1}} + \frac{1}{x_{n-1}-x_{n-2}}\right]  \nonumber
159
\end{align}
160

161
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. $\mathbf{C}$ becomes a discrete index $n$. 
162

163
In the following we write most expressions in their general vector form even though in one dimension considerable simplification is possible to connect with related articles where we use these expressions in more than one dimension.
164

165
Let $E =[x_{n}]$ represent the configuration space of one particle in 1-D space with world co-ordinate $n$. In Sage this can be represented as:
166
\begin{sageblock}
167
X = functions('x')
168
d = len(X)
169
C = (var('n'),)
170
E = tuple([x(*C) for x in X])
171
\end{sageblock}
172
\begin{equation} \sage{E} \end{equation}
173

174
\section{Forward/Backward Difference Approximation}
175

176
In this approximation to the derivative each world interacts only locally but with the
177
\emph{two} nearest neighbors to the left and right.
178

179
The forward and backward difference operators
180
\begin{sageblock}
181
def DM(x,n):return x-x.subs(n==n-1)
182
def DP(x,n):return x.subs(n==n+1)-x
183
\end{sageblock}
184
 approximate derivatives.
185

186
To derive Hall's equations for a finite number of interacting worlds we replace derivatives with difference operators so the Jacobian becomes a $1x1$ matrix
187
\begin{sageblock}
188
J = fun(n) (matrix([[DM(E[j],C[i]) for i in range(d)] for j in range(d)]))
189
\end{sageblock}
190
\begin{equation} \sage{J(n)} \end{equation}
191
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 just
192
\begin{sageblock}
193
A = fun(n) (J(n).adjoint())
194
\end{sageblock}
195
\begin{equation} \sage{A(n)} \end{equation}
196
so the inverse ``matrix'' is
197
\begin{sageblock}
198
K = fun(n) (A(n)/J_(n))
199
\end{sageblock}
200
\begin{equation} \sage{K(n)} \end{equation}
201
where the symbol $\J_n$ represents the determinant. We will eventually need to replace the symbolic determinants with then explicit expressions
202
\begin{sageblock}
203
_J = fun(n) (factor(J(n).determinant()))
204
dets = {J_(n+i):_J(n+i) for i in range(-2,3)}
205
\end{sageblock}
206
\begin{equation} \sage{_J(n)} \end{equation}
207
It is $\sage{factor(J(n)^(-1))==K(n).subs(J_(n)==_J(n))}$ that $J(n)^{-1}=K(n) \text{ where }{\sage{J_(n)}=\operatorname{_J}(n)}$.
208

209
\subsection{Quantum Potential}
210

211
The Schiff and Poirier expression for the quantum potential \eqref{eq20} can be approximated in terms of difference operations. Let $Q(n)=Q_2(n)+Q_3(n)$ be the 2nd and 3rd order terms, respectively.
212
\begin{sageblock}
213
Q = fun(n) ( -(hbar^2/4/m) * sum( sum( sum(
214
    (K(n)[k,j] * DM( DP( K(n)[l,j],C[l] ),C[k] ) +
215
        1/2 * DP(K(n)[l,j],C[l]) * DP(K(n)[k,j],C[k]))
216
            for k in range(d)) for j in range(d)) for l in range(d)) )
217
\end{sageblock}
218
\begin{dmath}
219
\sage{ expand(factor(Q(n))) }
220
\end{dmath}
221

222
\begin{sageblock}
223
Q2 = fun(n) ( -(hbar^2/4/m) * sum( sum( sum(
224
    (1/2 * DP(K(n)[l,j],C[l]) * DP(K(n)[k,j],C[k]) )
225
        for k in range(d)) for j in range(d)) for l in range(d)) )
226
\end{sageblock}
227
\begin{dmath}
228
\sage{ expand(factor(Q2(n))) }
229
\end{dmath}
230

231
\begin{sageblock}
232
Q3 = fun(n) ( -(hbar^2/4/m) * sum( sum( sum(
233
    (K(n)[k,j] * DM( DP( K(n)[l,j],C[l] ),C[k]) )
234
        for k in range(d)) for j in range(d)) for l in range(d)) )
235
\end{sageblock}
236
\begin{dmath}
237
\sage{ expand(factor(Q3(n))) }
238
\end{dmath}
239

240
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.
241

242
\begin{sageblock}
243
xyk = [X[i](n+j)
244
          for j in range(-2,3) for i in range(d)]
245
xyv = [var(X[i].name()+str(2+j))
246
          for j in range(-2,3) for i in range(d)]
247
xykv = dict(zip(xyk,xyv))
248
xyvk = dict(zip(xyv,xyk))
249
xn = X[0](n).subs(xykv)
250
\end{sageblock}
251
\begin{dmath}
252
\sage{(n==xn,)}
253
\end{dmath}
254

255
\subsection{Quantum Force}
256

257
The Schiff and Poirier expression for the quantum force is the last term on the left of \eqref{eq18}. It can be rewritten as terms of difference operations as follows
258
\begin{sageblock}
259
R = fun(n) (vector([ (hbar^2/4/m) * sum( sum( sum( sum(
260
   DP( K(n)[k,i] * K(n)[p,j] * DP( DM( K(n)[l,j],C[l] ),C[k] ),C[p] )
261
       for p in range(d)) for k in range(d)) for j in range(d)) for l in range(d))
262
           for i in range(d)]) )
263
\end{sageblock}
264
\begin{dmath}
265
\sage{ expand(R(n)[0]).subs(dets) }
266
\end{dmath}
267

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

270
To get the corresponding expression for the force we sum the potential symbolically over a sufficiently large area and then differentiate by $x(n)$. Since the potential is local, involving only nearby worlds, the result is independent of the size of the area so long as the area is large enough to include all the neighbors of a given world at $n$.
271

272
\noindent
273
Gradients:
274
\begin{sageblock}
275
r = fun(n) (vector([ -sum(Q(n+i) for i in range(-2,3)
276
                  ).subs(dets).subs(xykv).diff(xn).expand()
277
                      for k in range(d)]))
278
\end{sageblock}
279
\begin{dmath}
280
\sage{ expand(r(n).subs(xyvk)) }
281
\end{dmath}
282
with $r_2$ and $r_3$ the gradients of the 2nd and 3rd order terms of the potential, separately.
283
\begin{sageblock}
284
r2 = fun(n) (vector([ -sum(Q2(n+i) for i in range(-2,3)
285
          ).subs(dets).subs(xykv).diff(xn).expand()
286
              for k in range(d)]))
287
\end{sageblock}
288
\begin{dmath}
289
\sage{ expand(r2(n).subs(xyvk)) }
290
\end{dmath}
291

292
\begin{sageblock}
293
r3 = fun(n) (vector([ -sum([Q3(n+i) for i in range(-2,3)]
294
         ).subs(dets).subs(xykv).diff(xn).expand()
295
             for k in range(d)]))
296
\end{sageblock}
297
\begin{dmath}
298
\sage{ expand(r3(n).subs(xyvk)) }
299
\end{dmath}
300

301
Of course $r2(n)+r3(n)=r(n)$ is $\sage{r2(n)+r3(n)==r(n)}$.
302

303
Surprisingly, 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 differs from the 2nd term by only a sign.
304

305
It is $\sage{r3(n)==-2*r2(n)}$ that $r3(n)= -2\ r2(n)$.
306

307
This result is the same as the expression we obtained directly by substituting difference operators for derivatives in Poirier's expression eq. (18) for the quantum force. We will see later that this is not the case in more than one dimension.
308

309
It is $ \sage{bool(factor(r(n))==factor(R(n).subs(dets).subs(xykv)))} $ that $r(n)==R(n)$
310

311
Hall's equation \eqref{eq24} can be written
312

313
\begin{sageblock}
314
sigma = fun(n) ( (1/(X[0](n)-X[0](n-1))^2) *
315
    (1/(X[0](n+1)-X[0](n)) - 2/(X[0](n)-X[0](n-1)) + 1/(X[0](n-1)-X[0](n-2))) )
316
r_hall = fun(n) (hbar^2/4/m * (sigma(n+1) - sigma(n)))
317
\end{sageblock}
318
\begin{dmath}
319
\sage{ expand(r_hall(n)) }
320
\end{dmath}
321

322
It is $ \sage{bool(r(n)[0]==r_hall(n).subs(xykv))} $ that $rr=r_{hall}$.
323

324
Thus Hall's model of quantum mechanics as many interacting interacting classical worlds in one dimension follows directly from the approximation of Poirier's quantum mechanical equations of motion as difference equations with a certain choice of operators. Obviously this is not the only way to approximate Poirier's equations as difference equations.
325

326
Save these calculations for later.
327
\begin{sageblock}
328
Q_hall = Q(n).subs(dets)
329
R_hall = R(n).subs(dets)
330
xyv_hall = xyv
331
xykv_hall = xykv
332
\end{sageblock}
333

334
\section{Symmetric Difference Approximation}
335

336
The central difference approximation to this derivative is
337
\begin{sageblock}
338
def D(x,n):return x.subs(n==n+1)/2-x.subs(n==n-1)/2
339
\end{sageblock}
340

341
In this approximation we will see that each world interacts with its four neighbors to the left and right.  In fact we may continue this process of approximation to higher orders with more neighbors but in all these Hall-like models this non-locality in phase space is always bounded. Perhaps most interesting form the point of view of efficient computation is the lowest order approximation given previously. But the surprise is that all of these approximations gives quantum mechanics in the limit of a very large number of worlds (i.e. a continuum) and are good approximations even for a very small number of worlds.
342

343
Repeating the computations above using a symmetric difference approximation to the derivatives we again find that the quantum force computed as the gradient of the discretized quantum potential is identical to the discretization of Poirier's expression for the force.
344

345
\begin{sageblock}
346
J = fun(n) (matrix([[D(E[j],C[i]) for i in range(d)] for j in range(d)]))
347
\end{sageblock}
348
\begin{equation} \sage{J(n)} \end{equation}
349
\begin{sageblock}
350
A = fun(n) (fix(J(n).adjoint()))
351
K = fun(n) (A(n)/J_(n))
352
\end{sageblock}
353
\begin{equation} \sage{K(n)} \end{equation}
354
\begin{sageblock}
355
_J = fun(n) (fix(J(n).determinant().expand()))
356
dets = {J_(n+i):_J(n+i) for i in range(-6,7)}
357
\end{sageblock}
358
\begin{equation} \sage{_J(n)} \end{equation}
359

360
\begin{sageblock}
361
Q = fun(n) ( -(hbar^2/4/m) *
362
    sum( sum( sum(
363
        (K(n)[k,j] * D( D( K(n)[l,j],C[l] ),C[k] ) +
364
            1/2 * D(K(n)[l,j],C[l]) * D(K(n)[k,j],C[k]))
365
                for k in range(d)) for j in range(d)) for l in range(d)) )
366
Q2 = fun(n) ( -(hbar^2/4/m) *
367
    sum( sum( sum( (1/2 * D(K(n)[l,j],C[l]) * D(K(n)[k,j],C[k]))
368
        for k in range(d)) for j in range(d)) for l in range(d)) )
369
Q3 = fun(n) ( -(hbar^2/4/m) *
370
    sum( sum( sum( (K(n)[k,j] * D( D( K(n)[l,j],C[l] ),C[k] ))
371
        for k in range(d)) for j in range(d)) for l in range(d)) )
372
\end{sageblock}
373
\begin{dmath}
374
\sage{ expand(factor(Q(n))) }
375
\end{dmath}
376
\begin{sageblock}
377
R = fun(n) (vector([ (hbar^2/4/m)*sum( sum( sum( sum(
378
    D( K(n)[k,i] * K(n)[p,j] * D( D( K(n)[l,j],C[l] ),C[k] ),C[p] )
379
        for p in range(d)) for k in range(d)) for j in range(d)) for l in range(d))
380
            for i in range(d)]) )
381
\end{sageblock}
382
\begin{dmath}
383
\sage{ expand(R(n)[0].subs(dets)) }
384
\end{dmath}
385

386
\begin{sageblock}
387
xyk = [X[i](n+j) for j in range(-4,5) for i in range(d)]
388
xyv = [var(X[i].name()+str(5+j)) for j in range(-4,5) for i in range(d)]
389
xykv = dict(zip(xyk,xyv))
390
xn = X[0](n).subs(xykv)
391
\end{sageblock}
392

393
\begin{dmath}
394
\sage{(n==xn)}
395
\end{dmath}
396

397
The neighborhood for the symmetric difference operators is larger than previous case.
398
\begin{sageblock}
399
r = fun(n) (vector([-sum(Q(n+i) for i in range(-3,4)
400
            ).subs(dets).subs(xykv).diff(xn).expand()
401
                for k in range(d)]))
402
r2 = fun(n) (vector([-sum(Q2(n+i) for i in range(-3,4)
403
            ).subs(dets).subs(xykv).diff(xn).expand()
404
                for k in range(d)]))
405
r3 = fun(n) (vector([-sum([Q3(n+i) for i in range(-3,4)]
406
            ).subs(dets).subs(xykv).diff(xn).expand()
407
               for k in range(d)]))
408
\end{sageblock}
409

410
It is $\sage{r3(n)==-2*r2(n)}$ that $r3(n)= -2\ r2(n)$.
411

412
It is $ \sage{bool(factor(r(n))==factor(R(n).subs(dets).subs(xykv)))} $ that $r(n)=R(n)$
413

414
This expression is the same as the expression we obtained directly by substituting difference operators for derivatives in Poirier's expression \eqref{eq18} for the quantum force.
415

416
Save these calculations for later.
417
\begin{sageblock}
418
Q_poirier = Q(n).subs(dets)
419
R_poirier = R(n).subs(dets)
420
xyv_poirer = xyv
421
xykv_poirier = xykv
422
\end{sageblock}
423

424

425
\section{Simulation}
426

427
Using the expressions derived above we can generate efficient numerical functions which we use later in determining the motion of a particle in many worlds simultaneously.
428

429
We will use Hall's difference approximation to Poirier's equations above for a numerical simulation of quantum motion in a small number of ``parallel'' worlds.
430

431
For efficiency we expression the following operations in  ``vectorized'' form. First we need the difference operator
432
\begin{sageblock}
433
def dM(x,n):return shift(x,n) - shift(x,n-1)
434
\end{sageblock}
435

436

437
\subsection{Classical Potential}
438
Poirier's problem involves scattering a proton-mass particle from an Eckart potential slightly less than the particles incident energy.
439
\begin{sageblock}
440
alf = 2.5 ; x0=0 ; B = 0.0024
441
V_(x) = B /cosh(alf *(x-x0))^2
442
#V_(x) = 2/cosh(2*x)^2
443
#V_(x) = 0
444
V = vectorize(fast_callable( V_(x), vars=[x], domain=RDF))
445
_ = plot(V,(-3,3),axes_labels=['distance','energy'],axes_labels_size=1,
446
         legend_label='potential')
447
\end{sageblock}
448
\begin{figure}[H] \centering
449
\sageplot[width=.9\textwidth]{_}
450
\caption{Classical Potential}
451
\end{figure}
452

453
\subsection{Classical Force}
454
\begin{sageblock}
455
F = vectorize(fast_callable( -diff(V_(x),x), vars=[x], domain=RDF))
456
_ =plot(F,(-3,3),axes_labels=['distance','force'],axes_labels_size=1,
457
        legend_label='force')
458
\end{sageblock}
459
\begin{figure}[H] \centering
460
\sageplot[width=.9\textwidth]{_}
461
\caption{Classical Force}
462
\end{figure}
463

464
\subsection{Quantum Potential}
465

466
\begin{sagesilent}[fancy]
467
def qq(hbar,m,x):
468
    dx0  = dM(x,0)
469
    dxp1 = dM(x,+1)
470
    dxm1 = dM(x,-1)
471
    return hbar^2/4/m*(-1/2/dxp1^2 + 3/2/dx0^2 - 1/dxm1/dx0)
472
\end{sagesilent}
473
\includesagecollection{fancy}
474

475
\subsection{Quantum Force}
476
In general the force is a vector but in the following we treat it as a single scalar quantity.
477

478
Vectorized Quantum Force
479
\begin{sagesilent}[fancy]
480
def ff(hbar,m,x):
481
    dx0  = dM(x,0)
482
    dxp1 = dM(x,+1)
483
    dxp2 = dM(x,+2)
484
    dxm1 = dM(x,-1)
485
    return hbar^2/2/m*(-1/dxp1^3 + 1/2/dxp2/dxp1^2 + 1/dx0^3 - 1/2/dxp1/dx0^2 - 1/2/dxm1/dx0^2 + 1/2/dxp1^2/dx0)
486
\end{sagesilent}]
487
\includesagecollection{fancy}
488

489
\subsection{Initial Data}
490

491
\subsubsection{Spatial Distribution}
492

493
In order to compute the quantum force on the particles in world $n$ we introduce four widely separated fictitious particles in worlds to the left and right to represent open boundary conditions.
494
\begin{sagesilent}[fancy]
495
def boundary(x): # returns x as a view over a base with boundary values
496
    return concatenate((
497
        array([-1e18,-1e14,-1e10,-1e6]),
498
        x,
499
        array([1e6,1e10,1e14,1e18])))[4:-4]
500
\end{sagesilent}
501
\includesagecollection{fancy}
502

503
For example consider the spatial distribution of a Gaussian ensemble of one particle in  $N$ ``parallel'' worlds. Everything is in atomic units below, so $\hbar=1$.
504
\begin{sagesilent}[fancy]
505
N = 2000
506
var('alf,x0,m,v,x1')
507
model1={m:RDF(2000),alf:RDF(0.70),v:RDF(0.00164317),x0:RDF(-7.0)}
508
W(x)=abs((alf/pi)^(1/4)*exp(-(alf/2)*(x-x0)^2)*exp(I*m*v*x))^2
509
\end{sagesilent}
510
\includesagecollection{fancy}
511
\begin{equation}
512
\sage{W(x)}
513
\end{equation}
514
We may compute a uniform mapping into this distribution by inverting the cummulative distribution.
515
\begin{sagesilent}[fancy]
516
assume(alf>0)
517
erf0=integrate(W(x),x,-oo,x1).subs(model1)
518
y1=var('y1');inverse_erf=function('inverse_erf')
519
erfinv0=solve(
520
    integrate(W(x),x,-oo,x1).subs(model1)==y1,x1)[0].rhs()
521
def erfinv1(x): return erfinv0.subs(y1==x).substitute_function(inverse_erf,erfinv)
522
\end{sagesilent}
523
\includesagecollection{fancy}
524
After adding the open boundary conditions, we can check the accuracy by comparing the discrete representation to the original distribution.
525
\begin{sageblock}
526
p0 = array([erfinv1((i+0.001)/(N-1+0.002)) for i in range(N)], dtype=rdf)
527
b0 = boundary(p0)
528
_ = (list_plot([[b0[i],2/N/(b0[i+1]-b0[i-1])] for i in range(4,len(b0)-4) ],
529
              axes_labels=['location','density'], color='red',
530
              legend_label='discrete')+
531
     plot(W(x).subs(model1),(x,-15,5),
532
              legend_label='continuous'))
533
\end{sageblock}
534
\begin{figure}[H] \centering
535
\sageplot[width=.9\textwidth]{_}
536
\caption{Initial Spatial Density}
537
\end{figure}
538

539
\subsubsection{Distribution of Velocity}
540

541
The particles have an identical initial velocity in each world of $\sage{v.subs(model1)}$.
542
\begin{sageblock}
543
v0=array([v.subs(model1) for i in range(N)],dtype=rdf)
544
\end{sageblock}
545

546
\subsection{Energy of the Ensemble}
547
\noindent
548
Total classical potential
549
\begin{sageblock}
550
def PU(m,x): return sum(V(i) for i in x)
551
\end{sageblock}
552
\noindent
553
Total kinetic energy
554
\begin{sageblock}
555
def KU(m,v): return RDF(0.5)*m*(v.dot(v))
556
\end{sageblock}
557
\noindent
558
Total quantum potential
559
\begin{sageblock}
560
def QU(m,hbar,x): return qq(m,hbar,x).sum()
561
\end{sageblock}
562

563
\subsection{Acceleration}
564

565
Acceleration is the combined effect of both classical and quantum forces.
566
\begin{sageblock}
567
def A(hbar,m,x): return (F(x) + ff(hbar,m,x))/m
568
\end{sageblock}
569

570
\subsection{Step Size Controller}
571

572
Continuous integrating step size controller \cite{hairer2005}
573
\begin{sagesilent}[fancy]
574
def G(a,v):
575
    alpha = RDF(-100.0) # sensitivity
576
    eps = RDF(10.0^-10)
577
    return -alpha*(a.dot(v))/max(v.dot(v),eps)
578
\end{sagesilent}
579
\includesagecollection{fancy}
580

581
\subsection{Integration}
582
Integrate using the Störmer–Verlet algorithm with adaptive step size.
583

584
\begin{sagesilent}[fancy]
585
t_end = 11600; t_samples = 100; x_samples = 100
586
t_sample = t_end/t_samples
587
x_sample = max(1,int(N)/int(x_samples))
588
x_start = max(0,int(N-x_samples*x_sample)/int(2))
589
XS = range(x_start,N-x_start-1,x_sample)
590
ds0 = RDF(0.007) # initial step size
591
dsn = 1 # initial step size divider
592
xb = boundary(p0)
593
while True:
594
    try: # step size
595
        ds = ds0/dsn
596
        rho = RDF(1.0) # step density
597
        dt = ds/rho
598
        m0 = RDF(2000.0)
599
        hbar0 = RDF(1.0)
600
        # initial values
601
        x = boundary(p0); v = v0.copy(); a = A(hbar0,m0,x)
602
        rho = rho - RDF(0.5)*G(a,v)*dt
603
        t = 0;  T = [t]
604
        XX = {t:x[XS].copy()}; XV = {t:v[XS].copy()}
605
        KK = KU(m0,v); PP = PU(m0,x); QQ = QU(hbar0,m0,x)
606
        E1 = KK+PP+QQ  # total energy
607
        Rho = {t:rho}; TK = {t:KK}; TP = {t:PP}; TQ = {t:QQ}; E = {t:E1}
608
        try:
609
            while t<t_end:
610
                print "t =", t, "rho =", rho
611
                t1 = t + t_sample
612
                while t<t1:
613
                    rho = rho + G(a,v)*dt
614
                    dt = ds/rho; t += dt
615
                    x += v*dt + 0.5*a*dt^2
616
                    v += 0.5*a*dt; a = A(hbar0,m0,x); v += 0.5*a*dt
617
                    if isnan(rho) or rho>100 or rho<0.01:
618
                        raise ValueError(
619
                            "Step control failed at %s. rho=%s"%(t,rho))
620
                # Check no-crossing
621
                if (x[1:]-x[:N-1]).min()<0:
622
                    raise ValueError("crossing")
623
                # Check total energy conservation
624
                KK = KU(m0,v); PP = PU(m0,x); QQ = QU(hbar0,m0,x)
625
                E2 = KK+PP+QQ
626
                if abs(E1-E2)>0.01:
627
                    raise ValueError(
628
                        "Energy conservation bound failed at %s. Delta E:"\
629
                        "|%s| > 0.01."%(t,E1-E2))
630
                E1 = E2
631
                T += [t]; XX[t] = x[XS].copy(); XV[t] = v[XS].copy()
632
                Rho[t] = rho; TK[t] = KK; TP[t] = PP; TQ[t] = QQ; E[t] = E1
633
        except KeyboardInterrupt:
634
            print "Interrupted at %s ..."%(t)
635
            T += [t]; XX[t] = x[XS].copy(); XV[t] = v[XS].copy()
636
            Rho[t] = rho; TK[t] = KK; TP[t] = PP; TQ[t] = QQ; E[t] = E1
637
        tmax = t
638
        break
639
    except ValueError as msg:
640
        print msg
641
        dsn = dsn + 1
642
        print "Trying a shorter initial step size: %s."%(ds0/dsn)
643
        continue
644
print "t =", tmax, "rho =",rho
645
\end{sagesilent}
646
\includesagecollection{fancy}
647

648
The adaptation of the integration stepsize during the simulation is given by \textbf{\rm rho}. A value greater than 1 indicates a smaller stepsize.
649
\begin{sageblock}
650
_ = line2d([[t,Rho[t]] for t in T],axes_labels=['time','rho'])
651
\end{sageblock}
652
\begin{figure}[H] \centering
653
\sageplot[width=.9\textwidth]{_}
654
\caption{Step Density}
655
\end{figure}
656

657
Energy is exchanged between classical, quantum mechanical and kinetic forms due to the interaction of the particle with the potential barrier (in each world) as well as the interaction between worlds.
658
\begin{sageblock}
659
_ = (line2d([[t,TK[t]] for t in T], axes_labels=['time','energy'],
660
        legend_label='kinetic', color='green') +
661
    line2d([[t,TQ[t]] for t in T], axes_labels=['time','energy'],
662
        legend_label='quantum potential',color='red') +
663
    line2d([[t,TP[t]] for t in T], axes_labels=['time','energy'],
664
        legend_label='classical potential',color='blue'))
665
\end{sageblock}
666
\begin{figure}[H] \centering
667
\sageplot[width=.9\textwidth]{_}
668
\caption{Quantum Potential, Classical Potential and Kinetic Energy}
669
\end{figure}
670

671
Conservation of energy is an important indicator of the quality of the numerical integration.
672
Note in particular how the changes in total energy over time remain well-bounded.
673
\begin{sageblock}
674
_ = line2d([[t,E[t]] for t in T],axes_labels=['time','energy'],
675
        legend_label='total energy',color='black')
676
\end{sageblock}
677
\begin{figure}[H] \centering
678
\sageplot[width=.9\textwidth]{_}
679
\caption{Energy Conservation}
680
\end{figure}
681

682
The trajectory of the particle in each world is shown below.  Notice how the trajectories do not cross, yet in some worlds the particle is scattered from the barrier while in others the particle tunnels across the barrier ``pushed'' by its counterpart in other worlds.
683
\begin{sageblock}
684
_ = (line([[t,XX[t].mean()] for t in T], color='black',
685
          thickness=2, legend_label="average",
686
          axes_labels=['time','location'], axes_labels_size=1)+
687
     sum(line([[t,XX[t][i]] for t in T], color=hue((0.5*i)/len(XS)), alpha=0.5,
688
              axes_labels=['time','location'], axes_labels_size=1)
689
                  for i in range(len(XS))))
690
\end{sageblock}
691
\begin{figure}[H] \centering
692
\sageplot[width=.9\textwidth]{_}
693
\caption{Trajectories}
694
\end{figure}
695

696
The distribution at the end of the simulation shows the reflected and transmitted "wave packets" separating from each other.
697
\begin{sageblock}
698
_ = list_plot([[x[i-1],2/N/(x[i+1]-x[i-1])] for i in range(1,N-1) ],
699
              axes_labels=['location','density'])
700
\end{sageblock}
701
\begin{figure}[H] \centering
702
\sageplot[width=.9\textwidth]{_}
703
\caption{Final Spatial Density}
704
\end{figure}
705

706
Transmission probability
707
\begin{sageblock}
708
_ = line([[t,(XV[t]>0).mean()] for t in T ],
709
#_ = line([[t,(XX[t]>0).mean()] for t in T ],
710
         axes_labels=['time','transmission'])
711
\end{sageblock}
712
\begin{figure}[H] \centering
713
\sageplot[width=.9\textwidth]{_}
714
\caption{Evolution of Spatial Density}
715
\end{figure}
716
Final $\sage{(v>0).mean()}$.
717

718
\newpage
719
\begin{appendices}
720
\section{Sage Preliminaries}\label{SageAppendix}
721

722
It is convenient to define a compact notation for functions:
723
\includesagecollection{prelim}
724
\includecollection{prelim1}
725
\includesagecollection{prelim2}
726
This routine shifts a numpy array view left or right by $n$ cells.
727
\includesagecollection{shift}
728

729
\end{appendices}
730
\newpage
731
\listoffigures
732
\bibliographystyle{../bibstyles/utphys}
733
\noindent
734
\bibliography{../bibliography}
735

736
\end{document}
737

738

739
%sagemathcloud={"latex_command":"pdflatex -synctex=1 -interact=nonstopmode 'poirier-problem1.tex'"}
740