CoCalc Shared Filesmiw / discretization1.texOpen in CoCalc with one click!
Author: Bill Page
Views : 13
Description: Discretization of the Quantum Dynamical Equations of Motion - One Dimension
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'"}