| Download
Discretization of the Quantum Dynamical Equations of Motion - One Dimension
Project: Many Interacting Worlds
Path: miw/discretization1.tex
Views: 86\documentclass[notitlepage]{article}1\usepackage{titling}2\usepackage[utf8]{inputenc}3\usepackage{titling}4\usepackage{fullpage}5\usepackage{mathtools}6\usepackage{etoolbox}7\cslet{lrbox*}\lrbox8\expandafter\patchcmd\csname lrbox*\endcsname{\setbox}{\global\setbox}{}{}9%\expandafter\show\csname lrbox*\endcsname % uncomment to see if it has worked10\cslet{endlrbox*}\endlrbox11\usepackage{sagetex}12\usepackage{url}13%\usepackage{amsmath}14\usepackage{cite}15\usepackage{xcolor}16\usepackage{breqn}17\usepackage{collect}18\setlength{\parskip}{1ex}19\usepackage[titletoc,title]{appendix}20\usepackage{float}21\usepackage{hyperref}22\hypersetup{23colorlinks,24linkcolor={red!50!black},25citecolor={blue!50!black},26urlcolor={blue!80!black}27}2829\title{Discretization of the Quantum Dynamical Equations of Motion - One Dimension}30\author{Bill Page \small{(\href{mailto:bill.page@newsynthesis.org}{bill.page@newsynthesis.org})}}31\date{Draft \today}3233\begin{document}3435\maketitle36\thispagestyle{empty}3738\begin{abstract}39\textbf{From Poirier’s Bohmian Mechanics without Wavefunctions, to Hall’s Many Interacting Worlds40in One Dimension}4142In 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.4344We provide a numerical simulation that reproduces the results given previously by Poirier\cite{poirier2008}.45\end{abstract}4647%\setcounter{section}{-1}4849\section{Introduction}5051In 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.5253In 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.5455From 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!5657In 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.5859In 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.6061\newcommand{\J}{\begin{vmatrix}\scriptstyle{J}\end{vmatrix}}6263\definecollection{prelim}64\begin{sagesilent}[prelim]65import numpy66from numpy import array,concatenate,vectorize,isnan67rdf=numpy.array([RDF.pi()]).dtype # ensure compatible types68import mpmath69from mpmath import erfinv7071def fun(*y): # fun(x,y)(f)(x,y)=f72cache={}73def upd(f,x,y): # memoized74try: return cache[x]75except KeyError: # not found76try: cache[x] = f.subs(dict([[y[i],x[i]]77for i in range(len(y))]))78except AttributeError: # ask forgiveness79def sub(f): return f.subs(dict([[y[i],x[i]]80for i in range(len(y))]))81cache[x] = map(sub,f)82return cache[x]83return lambda f: lambda *x: upd(f,x,y)84def argscript(self, *args):85return "%s_{%s}"%(self.name(),','.join(map(repr, args)))86def functions(*names):87return map(lambda nam:function(nam, print_latex_func=argscript),names)88def detJscript(self, *args):89return '\J_{%s}'% (','.join(map(repr, args)))90J__= function('J_') # determinant of J91J_ = function('J_', print_latex_func=detJscript)92def squared(x): # x^293try: return x.dot(x) # numpy array94except: return power(x,2)95\end{sagesilent}9697\definecollection{shift}98\begin{sagesilent}[shift]99def shift(x,n):100try: # deduce slicing101start = (x.__array_interface__['data'][0] -102x.base.__array_interface__['data'][0])/x.itemsize103stop = start + x.shape[0]104if start+n<0 or stop+n>x.base.shape[0]:105raise RuntimeError('bounds')106return x.base[start+n:stop+n]107except AttributeError:108raise RuntimeError('expecting view')109\end{sagesilent}110111112In this document we have used the computer algebra system113Sage\footnote{\sagestr{version()}\cite{SageMath2016} on SageMathCloud\cite{SageMathCloud2016}} with components of114SciPy\footnote{numpy Version \sagestr{numpy.version.version} \cite{SciPy2016}}115and116mpmath\footnote{mpmath Version \sagestr{mpmath.__version__} \cite{mpmath2016}}117for both symbolic and numeric calculations. This document and the source code for this work can be found118online\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}.119120For example we define some special symbols for Planck's constant and mass:121\begin{sageblock}122hbar = var('hbar',latex_name=r'\hbar')123m = var('m')124\end{sageblock}125126\definecollection{fancy}127\definecollection{prelim1}128\begin{collect}{prelim1}{\noindent}{\par}129The 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.130\end{collect}131\definecollection{prelim2}132\begin{sagecollect}{prelim2}133def fix(it):134ret = it.subs(dict([[function('x')(n+i), X[0](n+i)] for i in range(-2,3)] +135[[J__(n+i), J_(n+i)] for i in range(-2,3)]))136return ret137def factor(x):138try: return fix(x.factor())139except: return vector(map(lambda y:fix(y.factor()),x))140\end{sagecollect}141142\section{One Dimensional System}143144We show below that in one dimension Schiff and Poirier's quantum dynamical equation of motion, \cite{schiff2012communication} eqs. (18-20),145\begin{equation} \label{eq18}146m \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)=0147\end{equation}148that follows from a quantum potential of the form149\begin{equation} \label{eq20}150Q = -\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)151\end{equation}152153when 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.154155\begin{align}156r_N(x_n;X) &= \frac{\hbar^2}{4 m}[\sigma_{n+1}(X)-\sigma_n(X)] \label{eq24} \\157\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] \nonumber158\end{align}159160By 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$.161162In 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.163164Let $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:165\begin{sageblock}166X = functions('x')167d = len(X)168C = (var('n'),)169E = tuple([x(*C) for x in X])170\end{sageblock}171\begin{equation} \sage{E} \end{equation}172173\section{Forward/Backward Difference Approximation}174175In this approximation to the derivative each world interacts only locally but with the176\emph{two} nearest neighbors to the left and right.177178The forward and backward difference operators179\begin{sageblock}180def DM(x,n):return x-x.subs(n==n-1)181def DP(x,n):return x.subs(n==n+1)-x182\end{sageblock}183approximate derivatives.184185To derive Hall's equations for a finite number of interacting worlds we replace derivatives with difference operators so the Jacobian becomes a $1x1$ matrix186\begin{sageblock}187J = fun(n) (matrix([[DM(E[j],C[i]) for i in range(d)] for j in range(d)]))188\end{sageblock}189\begin{equation} \sage{J(n)} \end{equation}190It 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 just191\begin{sageblock}192A = fun(n) (J(n).adjoint())193\end{sageblock}194\begin{equation} \sage{A(n)} \end{equation}195so the inverse ``matrix'' is196\begin{sageblock}197K = fun(n) (A(n)/J_(n))198\end{sageblock}199\begin{equation} \sage{K(n)} \end{equation}200where the symbol $\J_n$ represents the determinant. We will eventually need to replace the symbolic determinants with then explicit expressions201\begin{sageblock}202_J = fun(n) (factor(J(n).determinant()))203dets = {J_(n+i):_J(n+i) for i in range(-2,3)}204\end{sageblock}205\begin{equation} \sage{_J(n)} \end{equation}206It 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)}$.207208\subsection{Quantum Potential}209210The 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.211\begin{sageblock}212Q = fun(n) ( -(hbar^2/4/m) * sum( sum( sum(213(K(n)[k,j] * DM( DP( K(n)[l,j],C[l] ),C[k] ) +2141/2 * DP(K(n)[l,j],C[l]) * DP(K(n)[k,j],C[k]))215for k in range(d)) for j in range(d)) for l in range(d)) )216\end{sageblock}217\begin{dmath}218\sage{ expand(factor(Q(n))) }219\end{dmath}220221\begin{sageblock}222Q2 = fun(n) ( -(hbar^2/4/m) * sum( sum( sum(223(1/2 * DP(K(n)[l,j],C[l]) * DP(K(n)[k,j],C[k]) )224for k in range(d)) for j in range(d)) for l in range(d)) )225\end{sageblock}226\begin{dmath}227\sage{ expand(factor(Q2(n))) }228\end{dmath}229230\begin{sageblock}231Q3 = fun(n) ( -(hbar^2/4/m) * sum( sum( sum(232(K(n)[k,j] * DM( DP( K(n)[l,j],C[l] ),C[k]) )233for k in range(d)) for j in range(d)) for l in range(d)) )234\end{sageblock}235\begin{dmath}236\sage{ expand(factor(Q3(n))) }237\end{dmath}238239In 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.240241\begin{sageblock}242xyk = [X[i](n+j)243for j in range(-2,3) for i in range(d)]244xyv = [var(X[i].name()+str(2+j))245for j in range(-2,3) for i in range(d)]246xykv = dict(zip(xyk,xyv))247xyvk = dict(zip(xyv,xyk))248xn = X[0](n).subs(xykv)249\end{sageblock}250\begin{dmath}251\sage{(n==xn,)}252\end{dmath}253254\subsection{Quantum Force}255256The 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 follows257\begin{sageblock}258R = fun(n) (vector([ (hbar^2/4/m) * sum( sum( sum( sum(259DP( K(n)[k,i] * K(n)[p,j] * DP( DM( K(n)[l,j],C[l] ),C[k] ),C[p] )260for p in range(d)) for k in range(d)) for j in range(d)) for l in range(d))261for i in range(d)]) )262\end{sageblock}263\begin{dmath}264\sage{ expand(R(n)[0]).subs(dets) }265\end{dmath}266267We 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.268269To 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$.270271\noindent272Gradients:273\begin{sageblock}274r = fun(n) (vector([ -sum(Q(n+i) for i in range(-2,3)275).subs(dets).subs(xykv).diff(xn).expand()276for k in range(d)]))277\end{sageblock}278\begin{dmath}279\sage{ expand(r(n).subs(xyvk)) }280\end{dmath}281with $r_2$ and $r_3$ the gradients of the 2nd and 3rd order terms of the potential, separately.282\begin{sageblock}283r2 = fun(n) (vector([ -sum(Q2(n+i) for i in range(-2,3)284).subs(dets).subs(xykv).diff(xn).expand()285for k in range(d)]))286\end{sageblock}287\begin{dmath}288\sage{ expand(r2(n).subs(xyvk)) }289\end{dmath}290291\begin{sageblock}292r3 = fun(n) (vector([ -sum([Q3(n+i) for i in range(-2,3)]293).subs(dets).subs(xykv).diff(xn).expand()294for k in range(d)]))295\end{sageblock}296\begin{dmath}297\sage{ expand(r3(n).subs(xyvk)) }298\end{dmath}299300Of course $r2(n)+r3(n)=r(n)$ is $\sage{r2(n)+r3(n)==r(n)}$.301302Surprisingly, 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.303304It is $\sage{r3(n)==-2*r2(n)}$ that $r3(n)= -2\ r2(n)$.305306This 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.307308It is $ \sage{bool(factor(r(n))==factor(R(n).subs(dets).subs(xykv)))} $ that $r(n)==R(n)$309310Hall's equation \eqref{eq24} can be written311312\begin{sageblock}313sigma = fun(n) ( (1/(X[0](n)-X[0](n-1))^2) *314(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))) )315r_hall = fun(n) (hbar^2/4/m * (sigma(n+1) - sigma(n)))316\end{sageblock}317\begin{dmath}318\sage{ expand(r_hall(n)) }319\end{dmath}320321It is $ \sage{bool(r(n)[0]==r_hall(n).subs(xykv))} $ that $rr=r_{hall}$.322323Thus 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.324325Save these calculations for later.326\begin{sageblock}327Q_hall = Q(n).subs(dets)328R_hall = R(n).subs(dets)329xyv_hall = xyv330xykv_hall = xykv331\end{sageblock}332333\section{Symmetric Difference Approximation}334335The central difference approximation to this derivative is336\begin{sageblock}337def D(x,n):return x.subs(n==n+1)/2-x.subs(n==n-1)/2338\end{sageblock}339340In 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.341342Repeating 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.343344\begin{sageblock}345J = fun(n) (matrix([[D(E[j],C[i]) for i in range(d)] for j in range(d)]))346\end{sageblock}347\begin{equation} \sage{J(n)} \end{equation}348\begin{sageblock}349A = fun(n) (fix(J(n).adjoint()))350K = fun(n) (A(n)/J_(n))351\end{sageblock}352\begin{equation} \sage{K(n)} \end{equation}353\begin{sageblock}354_J = fun(n) (fix(J(n).determinant().expand()))355dets = {J_(n+i):_J(n+i) for i in range(-6,7)}356\end{sageblock}357\begin{equation} \sage{_J(n)} \end{equation}358359\begin{sageblock}360Q = fun(n) ( -(hbar^2/4/m) *361sum( sum( sum(362(K(n)[k,j] * D( D( K(n)[l,j],C[l] ),C[k] ) +3631/2 * D(K(n)[l,j],C[l]) * D(K(n)[k,j],C[k]))364for k in range(d)) for j in range(d)) for l in range(d)) )365Q2 = fun(n) ( -(hbar^2/4/m) *366sum( sum( sum( (1/2 * D(K(n)[l,j],C[l]) * D(K(n)[k,j],C[k]))367for k in range(d)) for j in range(d)) for l in range(d)) )368Q3 = fun(n) ( -(hbar^2/4/m) *369sum( sum( sum( (K(n)[k,j] * D( D( K(n)[l,j],C[l] ),C[k] ))370for k in range(d)) for j in range(d)) for l in range(d)) )371\end{sageblock}372\begin{dmath}373\sage{ expand(factor(Q(n))) }374\end{dmath}375\begin{sageblock}376R = fun(n) (vector([ (hbar^2/4/m)*sum( sum( sum( sum(377D( K(n)[k,i] * K(n)[p,j] * D( D( K(n)[l,j],C[l] ),C[k] ),C[p] )378for p in range(d)) for k in range(d)) for j in range(d)) for l in range(d))379for i in range(d)]) )380\end{sageblock}381\begin{dmath}382\sage{ expand(R(n)[0].subs(dets)) }383\end{dmath}384385\begin{sageblock}386xyk = [X[i](n+j) for j in range(-4,5) for i in range(d)]387xyv = [var(X[i].name()+str(5+j)) for j in range(-4,5) for i in range(d)]388xykv = dict(zip(xyk,xyv))389xn = X[0](n).subs(xykv)390\end{sageblock}391392\begin{dmath}393\sage{(n==xn)}394\end{dmath}395396The neighborhood for the symmetric difference operators is larger than previous case.397\begin{sageblock}398r = fun(n) (vector([-sum(Q(n+i) for i in range(-3,4)399).subs(dets).subs(xykv).diff(xn).expand()400for k in range(d)]))401r2 = fun(n) (vector([-sum(Q2(n+i) for i in range(-3,4)402).subs(dets).subs(xykv).diff(xn).expand()403for k in range(d)]))404r3 = fun(n) (vector([-sum([Q3(n+i) for i in range(-3,4)]405).subs(dets).subs(xykv).diff(xn).expand()406for k in range(d)]))407\end{sageblock}408409It is $\sage{r3(n)==-2*r2(n)}$ that $r3(n)= -2\ r2(n)$.410411It is $ \sage{bool(factor(r(n))==factor(R(n).subs(dets).subs(xykv)))} $ that $r(n)=R(n)$412413This 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.414415Save these calculations for later.416\begin{sageblock}417Q_poirier = Q(n).subs(dets)418R_poirier = R(n).subs(dets)419xyv_poirer = xyv420xykv_poirier = xykv421\end{sageblock}422423424\section{Simulation}425426Using 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.427428We 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.429430For efficiency we expression the following operations in ``vectorized'' form. First we need the difference operator431\begin{sageblock}432def dM(x,n):return shift(x,n) - shift(x,n-1)433\end{sageblock}434435436\subsection{Classical Potential}437Poirier's problem involves scattering a proton-mass particle from an Eckart potential slightly less than the particles incident energy.438\begin{sageblock}439alf = 2.5 ; x0=0 ; B = 0.0024440V_(x) = B /cosh(alf *(x-x0))^2441#V_(x) = 2/cosh(2*x)^2442#V_(x) = 0443V = vectorize(fast_callable( V_(x), vars=[x], domain=RDF))444_ = plot(V,(-3,3),axes_labels=['distance','energy'],axes_labels_size=1,445legend_label='potential')446\end{sageblock}447\begin{figure}[H] \centering448\sageplot[width=.9\textwidth]{_}449\caption{Classical Potential}450\end{figure}451452\subsection{Classical Force}453\begin{sageblock}454F = vectorize(fast_callable( -diff(V_(x),x), vars=[x], domain=RDF))455_ =plot(F,(-3,3),axes_labels=['distance','force'],axes_labels_size=1,456legend_label='force')457\end{sageblock}458\begin{figure}[H] \centering459\sageplot[width=.9\textwidth]{_}460\caption{Classical Force}461\end{figure}462463\subsection{Quantum Potential}464465\begin{sagesilent}[fancy]466def qq(hbar,m,x):467dx0 = dM(x,0)468dxp1 = dM(x,+1)469dxm1 = dM(x,-1)470return hbar^2/4/m*(-1/2/dxp1^2 + 3/2/dx0^2 - 1/dxm1/dx0)471\end{sagesilent}472\includesagecollection{fancy}473474\subsection{Quantum Force}475In general the force is a vector but in the following we treat it as a single scalar quantity.476477Vectorized Quantum Force478\begin{sagesilent}[fancy]479def ff(hbar,m,x):480dx0 = dM(x,0)481dxp1 = dM(x,+1)482dxp2 = dM(x,+2)483dxm1 = dM(x,-1)484return 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)485\end{sagesilent}]486\includesagecollection{fancy}487488\subsection{Initial Data}489490\subsubsection{Spatial Distribution}491492In 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.493\begin{sagesilent}[fancy]494def boundary(x): # returns x as a view over a base with boundary values495return concatenate((496array([-1e18,-1e14,-1e10,-1e6]),497x,498array([1e6,1e10,1e14,1e18])))[4:-4]499\end{sagesilent}500\includesagecollection{fancy}501502For 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$.503\begin{sagesilent}[fancy]504N = 2000505var('alf,x0,m,v,x1')506model1={m:RDF(2000),alf:RDF(0.70),v:RDF(0.00164317),x0:RDF(-7.0)}507W(x)=abs((alf/pi)^(1/4)*exp(-(alf/2)*(x-x0)^2)*exp(I*m*v*x))^2508\end{sagesilent}509\includesagecollection{fancy}510\begin{equation}511\sage{W(x)}512\end{equation}513We may compute a uniform mapping into this distribution by inverting the cummulative distribution.514\begin{sagesilent}[fancy]515assume(alf>0)516erf0=integrate(W(x),x,-oo,x1).subs(model1)517y1=var('y1');inverse_erf=function('inverse_erf')518erfinv0=solve(519integrate(W(x),x,-oo,x1).subs(model1)==y1,x1)[0].rhs()520def erfinv1(x): return erfinv0.subs(y1==x).substitute_function(inverse_erf,erfinv)521\end{sagesilent}522\includesagecollection{fancy}523After adding the open boundary conditions, we can check the accuracy by comparing the discrete representation to the original distribution.524\begin{sageblock}525p0 = array([erfinv1((i+0.001)/(N-1+0.002)) for i in range(N)], dtype=rdf)526b0 = boundary(p0)527_ = (list_plot([[b0[i],2/N/(b0[i+1]-b0[i-1])] for i in range(4,len(b0)-4) ],528axes_labels=['location','density'], color='red',529legend_label='discrete')+530plot(W(x).subs(model1),(x,-15,5),531legend_label='continuous'))532\end{sageblock}533\begin{figure}[H] \centering534\sageplot[width=.9\textwidth]{_}535\caption{Initial Spatial Density}536\end{figure}537538\subsubsection{Distribution of Velocity}539540The particles have an identical initial velocity in each world of $\sage{v.subs(model1)}$.541\begin{sageblock}542v0=array([v.subs(model1) for i in range(N)],dtype=rdf)543\end{sageblock}544545\subsection{Energy of the Ensemble}546\noindent547Total classical potential548\begin{sageblock}549def PU(m,x): return sum(V(i) for i in x)550\end{sageblock}551\noindent552Total kinetic energy553\begin{sageblock}554def KU(m,v): return RDF(0.5)*m*(v.dot(v))555\end{sageblock}556\noindent557Total quantum potential558\begin{sageblock}559def QU(m,hbar,x): return qq(m,hbar,x).sum()560\end{sageblock}561562\subsection{Acceleration}563564Acceleration is the combined effect of both classical and quantum forces.565\begin{sageblock}566def A(hbar,m,x): return (F(x) + ff(hbar,m,x))/m567\end{sageblock}568569\subsection{Step Size Controller}570571Continuous integrating step size controller \cite{hairer2005}572\begin{sagesilent}[fancy]573def G(a,v):574alpha = RDF(-100.0) # sensitivity575eps = RDF(10.0^-10)576return -alpha*(a.dot(v))/max(v.dot(v),eps)577\end{sagesilent}578\includesagecollection{fancy}579580\subsection{Integration}581Integrate using the Störmer–Verlet algorithm with adaptive step size.582583\begin{sagesilent}[fancy]584t_end = 11600; t_samples = 100; x_samples = 100585t_sample = t_end/t_samples586x_sample = max(1,int(N)/int(x_samples))587x_start = max(0,int(N-x_samples*x_sample)/int(2))588XS = range(x_start,N-x_start-1,x_sample)589ds0 = RDF(0.007) # initial step size590dsn = 1 # initial step size divider591xb = boundary(p0)592while True:593try: # step size594ds = ds0/dsn595rho = RDF(1.0) # step density596dt = ds/rho597m0 = RDF(2000.0)598hbar0 = RDF(1.0)599# initial values600x = boundary(p0); v = v0.copy(); a = A(hbar0,m0,x)601rho = rho - RDF(0.5)*G(a,v)*dt602t = 0; T = [t]603XX = {t:x[XS].copy()}; XV = {t:v[XS].copy()}604KK = KU(m0,v); PP = PU(m0,x); QQ = QU(hbar0,m0,x)605E1 = KK+PP+QQ # total energy606Rho = {t:rho}; TK = {t:KK}; TP = {t:PP}; TQ = {t:QQ}; E = {t:E1}607try:608while t<t_end:609print "t =", t, "rho =", rho610t1 = t + t_sample611while t<t1:612rho = rho + G(a,v)*dt613dt = ds/rho; t += dt614x += v*dt + 0.5*a*dt^2615v += 0.5*a*dt; a = A(hbar0,m0,x); v += 0.5*a*dt616if isnan(rho) or rho>100 or rho<0.01:617raise ValueError(618"Step control failed at %s. rho=%s"%(t,rho))619# Check no-crossing620if (x[1:]-x[:N-1]).min()<0:621raise ValueError("crossing")622# Check total energy conservation623KK = KU(m0,v); PP = PU(m0,x); QQ = QU(hbar0,m0,x)624E2 = KK+PP+QQ625if abs(E1-E2)>0.01:626raise ValueError(627"Energy conservation bound failed at %s. Delta E:"\628"|%s| > 0.01."%(t,E1-E2))629E1 = E2630T += [t]; XX[t] = x[XS].copy(); XV[t] = v[XS].copy()631Rho[t] = rho; TK[t] = KK; TP[t] = PP; TQ[t] = QQ; E[t] = E1632except KeyboardInterrupt:633print "Interrupted at %s ..."%(t)634T += [t]; XX[t] = x[XS].copy(); XV[t] = v[XS].copy()635Rho[t] = rho; TK[t] = KK; TP[t] = PP; TQ[t] = QQ; E[t] = E1636tmax = t637break638except ValueError as msg:639print msg640dsn = dsn + 1641print "Trying a shorter initial step size: %s."%(ds0/dsn)642continue643print "t =", tmax, "rho =",rho644\end{sagesilent}645\includesagecollection{fancy}646647The adaptation of the integration stepsize during the simulation is given by \textbf{\rm rho}. A value greater than 1 indicates a smaller stepsize.648\begin{sageblock}649_ = line2d([[t,Rho[t]] for t in T],axes_labels=['time','rho'])650\end{sageblock}651\begin{figure}[H] \centering652\sageplot[width=.9\textwidth]{_}653\caption{Step Density}654\end{figure}655656Energy 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.657\begin{sageblock}658_ = (line2d([[t,TK[t]] for t in T], axes_labels=['time','energy'],659legend_label='kinetic', color='green') +660line2d([[t,TQ[t]] for t in T], axes_labels=['time','energy'],661legend_label='quantum potential',color='red') +662line2d([[t,TP[t]] for t in T], axes_labels=['time','energy'],663legend_label='classical potential',color='blue'))664\end{sageblock}665\begin{figure}[H] \centering666\sageplot[width=.9\textwidth]{_}667\caption{Quantum Potential, Classical Potential and Kinetic Energy}668\end{figure}669670Conservation of energy is an important indicator of the quality of the numerical integration.671Note in particular how the changes in total energy over time remain well-bounded.672\begin{sageblock}673_ = line2d([[t,E[t]] for t in T],axes_labels=['time','energy'],674legend_label='total energy',color='black')675\end{sageblock}676\begin{figure}[H] \centering677\sageplot[width=.9\textwidth]{_}678\caption{Energy Conservation}679\end{figure}680681The 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.682\begin{sageblock}683_ = (line([[t,XX[t].mean()] for t in T], color='black',684thickness=2, legend_label="average",685axes_labels=['time','location'], axes_labels_size=1)+686sum(line([[t,XX[t][i]] for t in T], color=hue((0.5*i)/len(XS)), alpha=0.5,687axes_labels=['time','location'], axes_labels_size=1)688for i in range(len(XS))))689\end{sageblock}690\begin{figure}[H] \centering691\sageplot[width=.9\textwidth]{_}692\caption{Trajectories}693\end{figure}694695The distribution at the end of the simulation shows the reflected and transmitted "wave packets" separating from each other.696\begin{sageblock}697_ = list_plot([[x[i-1],2/N/(x[i+1]-x[i-1])] for i in range(1,N-1) ],698axes_labels=['location','density'])699\end{sageblock}700\begin{figure}[H] \centering701\sageplot[width=.9\textwidth]{_}702\caption{Final Spatial Density}703\end{figure}704705Transmission probability706\begin{sageblock}707_ = line([[t,(XV[t]>0).mean()] for t in T ],708#_ = line([[t,(XX[t]>0).mean()] for t in T ],709axes_labels=['time','transmission'])710\end{sageblock}711\begin{figure}[H] \centering712\sageplot[width=.9\textwidth]{_}713\caption{Evolution of Spatial Density}714\end{figure}715Final $\sage{(v>0).mean()}$.716717\newpage718\begin{appendices}719\section{Sage Preliminaries}\label{SageAppendix}720721It is convenient to define a compact notation for functions:722\includesagecollection{prelim}723\includecollection{prelim1}724\includesagecollection{prelim2}725This routine shifts a numpy array view left or right by $n$ cells.726\includesagecollection{shift}727728\end{appendices}729\newpage730\listoffigures731\bibliographystyle{../bibstyles/utphys}732\noindent733\bibliography{../bibliography}734735\end{document}736737738%sagemathcloud={"latex_command":"pdflatex -synctex=1 -interact=nonstopmode 'poirier-problem1.tex'"}739740