CoCalc Public FilesHall2Dentangle1.sagewsOpen with one click!
Authors: David Cyganski, Bill Page
Views : 44
Description: Each moving point in the following diagram represents the position of two particles in a single world (four worlds are shown). As the point approaches the x=y line (shown in red), the particles approach each other and the classical potential is at a maximum. In three worlds the particles overcome the potential barrier and continue. But one world (the last to approach) is reflected back. Although initially the two particles in each world approach each other with the same velocity, the particles in the last world to approach the barrier have lost energy to their counterparts in other worlds due to the quantum mechanical (interworld) interaction. When viewed in the projection we see two particles rebound instead of passing each other.
Compute Environment: Ubuntu 18.04 (Deprecated)

From Poirier's Bohmian Mechanics without Wavefunctions, to Hall's Many Interacting Worlds -- in More Than One Dimension

Ref:

  1. Quantum Mechanics Without Wavefunctions Jeremy Schiff and Bill Poirier J. Chem. Phys. 136, 031102 (2012)

  2. Quantum Phenomena Modeled by Interactions between Many Classical Worlds Michael J. W. Hall, Dirk-André Deckert and Howard M. Wiseman, PHYSICAL REVIEW X 4, 041013 (23 October 2014)

  3. Verlet integration (Wikipedia)



  4. Explicit, Time Reversible, Adaptive Step Size Control Ernst Hairer and Gustaf Söderlind SIAM Journal on Scientific Computing. 2005, vol. 26, no. 6, p. 1838-1851

#%typeset_mode True from numpy import array,concatenate,isnan from mpmath import erfinv hbar=var('hbar',latex_name='\\hbar') mu=var('mu',latex_name='\mu')

2-D

Poirier’s many-D generalization, eq. (18), when specialized to 2-D and discretized by replacing derivatives with difference operators is the 2-D analog of Hall’s “toy” expression eq. (24). Let x(n,m)x(n,m) and y(n,m)y(n,m) represent the position of a particle (one particle in 2-D space) or the position of two particles (in 1-D configuration space) at an unspecified time in the world indexed by nn and mm.

vars = ['x','y']; d = len(vars) X = map(function,vars); x,y = X n,m = var('n,m'); ind=[n,m] # position of particle in world (n,m) Enm=map(lambda x:x(*ind),X);Enm
[x(n, m), y(n, m)]

For Hall's finite number of interacting worlds we will use the following difference operators instead of derivatives.

def Dminus(x,i):return(x-x.subs(ind[i]==ind[i]-1)) def Dplus(x,i):return(x.subs(ind[i]==ind[i]+1)-x)

The Jacobian replacing derivatives with differences.

Jnm = matrix(map(lambda e:[Dminus(e,i) for i in range(d)],Enm));Jnm
[-x(n - 1, m) + x(n, m) -x(n, m - 1) + x(n, m)] [-y(n - 1, m) + y(n, m) -y(n, m - 1) + y(n, m)]
# inverse Jacobian Knm = matrix(map(lambda x: map(lambda y:y.factor(),x),Jnm^(-1)))
# determinant Dnm=Jnm.determinant(); Dnm
-(x(n, m - 1) - x(n, m))*(y(n - 1, m) - y(n, m)) + (x(n - 1, m) - x(n, m))*(y(n, m - 1) - y(n, m))
# minors Mnm=matrix(map(lambda x: map(lambda y:(y*Dnm).full_simplify(),x),Knm));Mnm
[-y(n, m - 1) + y(n, m) x(n, m - 1) - x(n, m)] [ y(n - 1, m) - y(n, m) -x(n - 1, m) + x(n, m)]
# inverse is minors divided by determinant def Jdet(n,m): return Dnm.subs({var('n'):n,var('m'):m}) map(lambda x:map(lambda y:y.full_simplify(),x),Mnm/Dnm-Knm)==[[0,0],[0,0]]
True

Quantum Force

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

# introduce symbolic determinant det=function('det') eq18n = [ (hbar^2/4/mu)*sum([ sum([ sum([ sum([ Dplus( Mnm[k,i]/det(n,m)*Mnm[p,j]/det(n,m)*Dplus( Dminus( Mnm[l,j]/det(n,m),l ),k ),p ) for p in range(d)]) for k in range(d)]) for j in range(d)]) for l in range(d)]) for i in range(d)] # Too large to conveniently display # eq18n
# the denominators of both components are the same show(eq18n[0].denominator()) bool(eq18n[0].denominator()==eq18n[1].denominator())
4μdet(n+2,m1)det(n+2,m)det(n+1,m+1)det(n+1,m1)det(n+1,m)3det(n1,m+2)det(n1,m+1)det(n1,m)det(n,m+2)det(n,m+1)3det(n,m1)det(n,m)3\displaystyle 4 \, {{\mu}} {\rm det}\left(n + 2, m - 1\right) {\rm det}\left(n + 2, m\right) {\rm det}\left(n + 1, m + 1\right) {\rm det}\left(n + 1, m - 1\right) {\rm det}\left(n + 1, m\right)^{3} {\rm det}\left(n - 1, m + 2\right) {\rm det}\left(n - 1, m + 1\right) {\rm det}\left(n - 1, m\right) {\rm det}\left(n, m + 2\right) {\rm det}\left(n, m + 1\right)^{3} {\rm det}\left(n, m - 1\right) {\rm det}\left(n, m\right)^{3}
True
# expression components are rational functions eq18numer0=eq18n[0].numerator().subs({det(n+i,m+j):Jdet(n+i,m+j) for i in range(-1,3) for j in range(-1,3)}) eq18denom0=eq18n[0].denominator().subs({det(n+i,m+j):Jdet(n+i,m+j) for i in range(-1,3) for j in range(-1,3)}) eq18numer1=eq18n[1].numerator().subs({det(n+i,m+j):Jdet(n+i,m+j) for i in range(-1,3) for j in range(-1,3)}) eq18denom1=eq18n[1].denominator().subs({det(n+i,m+j):Jdet(n+i,m+j) for i in range(-1,3) for j in range(-1,3)})
# denominator explicitly show(eq18denom0)
4((x(n+2,m1)x(n+1,m1))(y(n+2,m1)y(n+2,m2))(x(n+2,m1)x(n+2,m2))(y(n+2,m1)y(n+1,m1)))((x(n+2,m)x(n+1,m))(y(n+2,m1)y(n+2,m))(x(n+2,m1)x(n+2,m))(y(n+2,m)y(n+1,m)))((x(n+1,m+1)x(n,m+1))(y(n+1,m+1)y(n+1,m))(x(n+1,m+1)x(n+1,m))(y(n+1,m+1)y(n,m+1)))((x(n+1,m1)x(n,m1))(y(n+1,m1)y(n+1,m2))(x(n+1,m1)x(n+1,m2))(y(n+1,m1)y(n,m1)))((x(n+1,m)x(n,m))(y(n+1,m1)y(n+1,m))(x(n+1,m1)x(n+1,m))(y(n+1,m)y(n,m)))3((x(n1,m+2)x(n2,m+2))(y(n1,m+2)y(n1,m+1))(x(n1,m+2)x(n1,m+1))(y(n1,m+2)y(n2,m+2)))((x(n,m+2)x(n,m+1))(y(n1,m+2)y(n,m+2))(x(n1,m+2)x(n,m+2))(y(n,m+2)y(n,m+1)))((x(n1,m+1)x(n2,m+1))(y(n1,m+1)y(n1,m))(x(n1,m+1)x(n1,m))(y(n1,m+1)y(n2,m+1)))((x(n,m+1)x(n,m))(y(n1,m+1)y(n,m+1))(x(n1,m+1)x(n,m+1))(y(n,m+1)y(n,m)))3((x(n1,m)x(n2,m))(y(n1,m1)y(n1,m))(x(n1,m1)x(n1,m))(y(n1,m)y(n2,m)))((x(n,m1)x(n,m2))(y(n1,m1)y(n,m1))(x(n1,m1)x(n,m1))(y(n,m1)y(n,m2)))((x(n,m1)x(n,m))(y(n1,m)y(n,m))(x(n1,m)x(n,m))(y(n,m1)y(n,m)))3μ\displaystyle 4 \, {\left({\left(x\left(n + 2, m - 1\right) - x\left(n + 1, m - 1\right)\right)} {\left(y\left(n + 2, m - 1\right) - y\left(n + 2, m - 2\right)\right)} - {\left(x\left(n + 2, m - 1\right) - x\left(n + 2, m - 2\right)\right)} {\left(y\left(n + 2, m - 1\right) - y\left(n + 1, m - 1\right)\right)}\right)} {\left({\left(x\left(n + 2, m\right) - x\left(n + 1, m\right)\right)} {\left(y\left(n + 2, m - 1\right) - y\left(n + 2, m\right)\right)} - {\left(x\left(n + 2, m - 1\right) - x\left(n + 2, m\right)\right)} {\left(y\left(n + 2, m\right) - y\left(n + 1, m\right)\right)}\right)} {\left({\left(x\left(n + 1, m + 1\right) - x\left(n, m + 1\right)\right)} {\left(y\left(n + 1, m + 1\right) - y\left(n + 1, m\right)\right)} - {\left(x\left(n + 1, m + 1\right) - x\left(n + 1, m\right)\right)} {\left(y\left(n + 1, m + 1\right) - y\left(n, m + 1\right)\right)}\right)} {\left({\left(x\left(n + 1, m - 1\right) - x\left(n, m - 1\right)\right)} {\left(y\left(n + 1, m - 1\right) - y\left(n + 1, m - 2\right)\right)} - {\left(x\left(n + 1, m - 1\right) - x\left(n + 1, m - 2\right)\right)} {\left(y\left(n + 1, m - 1\right) - y\left(n, m - 1\right)\right)}\right)} {\left({\left(x\left(n + 1, m\right) - x\left(n, m\right)\right)} {\left(y\left(n + 1, m - 1\right) - y\left(n + 1, m\right)\right)} - {\left(x\left(n + 1, m - 1\right) - x\left(n + 1, m\right)\right)} {\left(y\left(n + 1, m\right) - y\left(n, m\right)\right)}\right)}^{3} {\left({\left(x\left(n - 1, m + 2\right) - x\left(n - 2, m + 2\right)\right)} {\left(y\left(n - 1, m + 2\right) - y\left(n - 1, m + 1\right)\right)} - {\left(x\left(n - 1, m + 2\right) - x\left(n - 1, m + 1\right)\right)} {\left(y\left(n - 1, m + 2\right) - y\left(n - 2, m + 2\right)\right)}\right)} {\left({\left(x\left(n, m + 2\right) - x\left(n, m + 1\right)\right)} {\left(y\left(n - 1, m + 2\right) - y\left(n, m + 2\right)\right)} - {\left(x\left(n - 1, m + 2\right) - x\left(n, m + 2\right)\right)} {\left(y\left(n, m + 2\right) - y\left(n, m + 1\right)\right)}\right)} {\left({\left(x\left(n - 1, m + 1\right) - x\left(n - 2, m + 1\right)\right)} {\left(y\left(n - 1, m + 1\right) - y\left(n - 1, m\right)\right)} - {\left(x\left(n - 1, m + 1\right) - x\left(n - 1, m\right)\right)} {\left(y\left(n - 1, m + 1\right) - y\left(n - 2, m + 1\right)\right)}\right)} {\left({\left(x\left(n, m + 1\right) - x\left(n, m\right)\right)} {\left(y\left(n - 1, m + 1\right) - y\left(n, m + 1\right)\right)} - {\left(x\left(n - 1, m + 1\right) - x\left(n, m + 1\right)\right)} {\left(y\left(n, m + 1\right) - y\left(n, m\right)\right)}\right)}^{3} {\left({\left(x\left(n - 1, m\right) - x\left(n - 2, m\right)\right)} {\left(y\left(n - 1, m - 1\right) - y\left(n - 1, m\right)\right)} - {\left(x\left(n - 1, m - 1\right) - x\left(n - 1, m\right)\right)} {\left(y\left(n - 1, m\right) - y\left(n - 2, m\right)\right)}\right)} {\left({\left(x\left(n, m - 1\right) - x\left(n, m - 2\right)\right)} {\left(y\left(n - 1, m - 1\right) - y\left(n, m - 1\right)\right)} - {\left(x\left(n - 1, m - 1\right) - x\left(n, m - 1\right)\right)} {\left(y\left(n, m - 1\right) - y\left(n, m - 2\right)\right)}\right)} {\left({\left(x\left(n, m - 1\right) - x\left(n, m\right)\right)} {\left(y\left(n - 1, m\right) - y\left(n, m\right)\right)} - {\left(x\left(n - 1, m\right) - x\left(n, m\right)\right)} {\left(y\left(n, m - 1\right) - y\left(n, m\right)\right)}\right)}^{3} {{\mu}}

Create a numeric function that depends on just the nearest neighbors. This requires a function of 2x5x5 = 50 scalar quantities.

xyk=[X[i](n+j,m+k) for j in [-2,-1,0,1,2] for k in [-2,-1,0,1,2] for i in range(d)] xyv=[var(vars[i]+str(2+j)+str(2+k)) for j in [-2,-1,0,1,2] for k in [-2,-1,0,1,2] for i in range(d)] show(zip(xyk,xyv))
[(x(n2,m2)\displaystyle x\left(n - 2, m - 2\right), x00\displaystyle x_{00}), (y(n2,m2)\displaystyle y\left(n - 2, m - 2\right), y00\displaystyle y_{00}), (x(n2,m1)\displaystyle x\left(n - 2, m - 1\right), x01\displaystyle x_{01}), (y(n2,m1)\displaystyle y\left(n - 2, m - 1\right), y01\displaystyle y_{01}), (x(n2,m)\displaystyle x\left(n - 2, m\right), x02\displaystyle x_{02}), (y(n2,m)\displaystyle y\left(n - 2, m\right), y02\displaystyle y_{02}), (x(n2,m+1)\displaystyle x\left(n - 2, m + 1\right), x03\displaystyle x_{03}), (y(n2,m+1)\displaystyle y\left(n - 2, m + 1\right), y03\displaystyle y_{03}), (x(n2,m+2)\displaystyle x\left(n - 2, m + 2\right), x04\displaystyle x_{04}), (y(n2,m+2)\displaystyle y\left(n - 2, m + 2\right), y04\displaystyle y_{04}), (x(n1,m2)\displaystyle x\left(n - 1, m - 2\right), x10\displaystyle x_{10}), (y(n1,m2)\displaystyle y\left(n - 1, m - 2\right), y10\displaystyle y_{10}), (x(n1,m1)\displaystyle x\left(n - 1, m - 1\right), x11\displaystyle x_{11}), (y(n1,m1)\displaystyle y\left(n - 1, m - 1\right), y11\displaystyle y_{11}), (x(n1,m)\displaystyle x\left(n - 1, m\right), x12\displaystyle x_{12}), (y(n1,m)\displaystyle y\left(n - 1, m\right), y12\displaystyle y_{12}), (x(n1,m+1)\displaystyle x\left(n - 1, m + 1\right), x13\displaystyle x_{13}), (y(n1,m+1)\displaystyle y\left(n - 1, m + 1\right), y13\displaystyle y_{13}), (x(n1,m+2)\displaystyle x\left(n - 1, m + 2\right), x14\displaystyle x_{14}), (y(n1,m+2)\displaystyle y\left(n - 1, m + 2\right), y14\displaystyle y_{14}), (x(n,m2)\displaystyle x\left(n, m - 2\right), x20\displaystyle x_{20}), (y(n,m2)\displaystyle y\left(n, m - 2\right), y20\displaystyle y_{20}), (x(n,m1)\displaystyle x\left(n, m - 1\right), x21\displaystyle x_{21}), (y(n,m1)\displaystyle y\left(n, m - 1\right), y21\displaystyle y_{21}), (x(n,m)\displaystyle x\left(n, m\right), x22\displaystyle x_{22}), (y(n,m)\displaystyle y\left(n, m\right), y22\displaystyle y_{22}), (x(n,m+1)\displaystyle x\left(n, m + 1\right), x23\displaystyle x_{23}), (y(n,m+1)\displaystyle y\left(n, m + 1\right), y23\displaystyle y_{23}), (x(n,m+2)\displaystyle x\left(n, m + 2\right), x24\displaystyle x_{24}), (y(n,m+2)\displaystyle y\left(n, m + 2\right), y24\displaystyle y_{24}), (x(n+1,m2)\displaystyle x\left(n + 1, m - 2\right), x30\displaystyle x_{30}), (y(n+1,m2)\displaystyle y\left(n + 1, m - 2\right), y30\displaystyle y_{30}), (x(n+1,m1)\displaystyle x\left(n + 1, m - 1\right), x31\displaystyle x_{31}), (y(n+1,m1)\displaystyle y\left(n + 1, m - 1\right), y31\displaystyle y_{31}), (x(n+1,m)\displaystyle x\left(n + 1, m\right), x32\displaystyle x_{32}), (y(n+1,m)\displaystyle y\left(n + 1, m\right), y32\displaystyle y_{32}), (x(n+1,m+1)\displaystyle x\left(n + 1, m + 1\right), x33\displaystyle x_{33}), (y(n+1,m+1)\displaystyle y\left(n + 1, m + 1\right), y33\displaystyle y_{33}), (x(n+1,m+2)\displaystyle x\left(n + 1, m + 2\right), x34\displaystyle x_{34}), (y(n+1,m+2)\displaystyle y\left(n + 1, m + 2\right), y34\displaystyle y_{34}), (x(n+2,m2)\displaystyle x\left(n + 2, m - 2\right), x40\displaystyle x_{40}), (y(n+2,m2)\displaystyle y\left(n + 2, m - 2\right), y40\displaystyle y_{40}), (x(n+2,m1)\displaystyle x\left(n + 2, m - 1\right), x41\displaystyle x_{41}), (y(n+2,m1)\displaystyle y\left(n + 2, m - 1\right), y41\displaystyle y_{41}), (x(n+2,m)\displaystyle x\left(n + 2, m\right), x42\displaystyle x_{42}), (y(n+2,m)\displaystyle y\left(n + 2, m\right), y42\displaystyle y_{42}), (x(n+2,m+1)\displaystyle x\left(n + 2, m + 1\right), x43\displaystyle x_{43}), (y(n+2,m+1)\displaystyle y\left(n + 2, m + 1\right), y43\displaystyle y_{43}), (x(n+2,m+2)\displaystyle x\left(n + 2, m + 2\right), x44\displaystyle x_{44}), (y(n+2,m+2)\displaystyle y\left(n + 2, m + 2\right), y44\displaystyle y_{44})]
# suitable values constants = {hbar:0.2,mu:1.0} # compile numeric functions ff0=fast_callable(eq18numer0.subs(dict(zip(xyk,xyv))).subs(constants), vars=xyv, domain=RDF) ff1=fast_callable(eq18numer1.subs(dict(zip(xyk,xyv))).subs(constants), vars=xyv, domain=RDF) ffd=fast_callable(eq18denom0.subs(dict(zip(xyk,xyv))).subs(constants), vars=xyv, domain=RDF) def ff(p): den=ffd(*p); return array([ff0(*p)/den,ff1(*p)/den])

Quantum Potential###

Why does Poirier eq. (20) include a second term? If we include it then total energy T+V+QT+V+Q is not conserved!

eq20nm = sum([ sum([ sum([ (hbar^2/4/mu)*(Knm[k,j]*Dplus( Dplus( Knm[l,j],l ),k )+1/2*Dplus(Knm[l,j],l)*Dplus(Knm[k,j],k)) for k in range(d)]) for j in range(d)]) for l in range(d)]) #too complex to display conveniently #eq20nm

Expand terms using symbolic determinants.

eq20nm1 = sum([ sum([ sum([ (hbar^2/4/mu)*(Mnm[k,j]/det(n,m)*Dplus( Dplus( Mnm[l,j]/det(n,m),l ),k )) for k in range(d)]) for j in range(d)]) for l in range(d)]) show(eq20nm1.subs({hbar:1,mu:1}).expand().combine().collect(det(n,m)).collect(det(n+2,m)).collect(det(n+1,m+1)).collect(det(n+1,m)).collect(det(n,m+1)) )
x(n1,m)22x(n1,m)x(n,m1)+x(n,m1)2+y(n1,m)22y(n1,m)y(n,m1)+y(n,m1)24det(n,m)2+x(n+2,m1)x(n,m1)x(n+2,m)x(n,m1)x(n+2,m1)x(n,m)+x(n+2,m)x(n,m)+y(n+2,m1)y(n,m1)y(n+2,m)y(n,m1)y(n+2,m1)y(n,m)+y(n+2,m)y(n,m)4det(n+2,m)det(n,m)+x(n+1,m+1)x(n1,m)x(n+1,m)x(n1,m)+x(n+1,m+1)x(n,m1)x(n,m+1)x(n,m1)2x(n+1,m+1)x(n,m)+x(n+1,m)x(n,m)+x(n,m+1)x(n,m)+y(n+1,m+1)y(n1,m)y(n+1,m)y(n1,m)+y(n+1,m+1)y(n,m1)y(n,m+1)y(n,m1)2y(n+1,m+1)y(n,m)+y(n+1,m)y(n,m)+y(n,m+1)y(n,m)4det(n+1,m+1)det(n,m)+x(n+1,m1)x(n1,m)x(n+1,m)x(n1,m)2x(n+1,m1)x(n,m1)+x(n+1,m)x(n,m1)+x(n+1,m1)x(n,m)+x(n,m1)x(n,m)x(n,m)2+y(n+1,m1)y(n1,m)y(n+1,m)y(n1,m)2y(n+1,m1)y(n,m1)+y(n+1,m)y(n,m1)+y(n+1,m1)y(n,m)+y(n,m1)y(n,m)y(n,m)24det(n+1,m)det(n,m)+x(n1,m+2)x(n1,m)x(n1,m)x(n,m+2)x(n1,m+2)x(n,m)+x(n,m+2)x(n,m)+y(n1,m+2)y(n1,m)y(n1,m)y(n,m+2)y(n1,m+2)y(n,m)+y(n,m+2)y(n,m)4det(n,m+2)det(n,m)2x(n1,m+1)x(n1,m)x(n1,m)x(n,m+1)x(n1,m+1)x(n,m1)+x(n,m+1)x(n,m1)x(n1,m+1)x(n,m)x(n1,m)x(n,m)+x(n,m)2+2y(n1,m+1)y(n1,m)y(n1,m)y(n,m+1)y(n1,m+1)y(n,m1)+y(n,m+1)y(n,m1)y(n1,m+1)y(n,m)y(n1,m)y(n,m)+y(n,m)24det(n,m+1)det(n,m)\displaystyle \frac{x\left(n - 1, m\right)^{2} - 2 \, x\left(n - 1, m\right) x\left(n, m - 1\right) + x\left(n, m - 1\right)^{2} + y\left(n - 1, m\right)^{2} - 2 \, y\left(n - 1, m\right) y\left(n, m - 1\right) + y\left(n, m - 1\right)^{2}}{4 \, {\rm det}\left(n, m\right)^{2}} + \frac{x\left(n + 2, m - 1\right) x\left(n, m - 1\right) - x\left(n + 2, m\right) x\left(n, m - 1\right) - x\left(n + 2, m - 1\right) x\left(n, m\right) + x\left(n + 2, m\right) x\left(n, m\right) + y\left(n + 2, m - 1\right) y\left(n, m - 1\right) - y\left(n + 2, m\right) y\left(n, m - 1\right) - y\left(n + 2, m - 1\right) y\left(n, m\right) + y\left(n + 2, m\right) y\left(n, m\right)}{4 \, {\rm det}\left(n + 2, m\right) {\rm det}\left(n, m\right)} + \frac{x\left(n + 1, m + 1\right) x\left(n - 1, m\right) - x\left(n + 1, m\right) x\left(n - 1, m\right) + x\left(n + 1, m + 1\right) x\left(n, m - 1\right) - x\left(n, m + 1\right) x\left(n, m - 1\right) - 2 \, x\left(n + 1, m + 1\right) x\left(n, m\right) + x\left(n + 1, m\right) x\left(n, m\right) + x\left(n, m + 1\right) x\left(n, m\right) + y\left(n + 1, m + 1\right) y\left(n - 1, m\right) - y\left(n + 1, m\right) y\left(n - 1, m\right) + y\left(n + 1, m + 1\right) y\left(n, m - 1\right) - y\left(n, m + 1\right) y\left(n, m - 1\right) - 2 \, y\left(n + 1, m + 1\right) y\left(n, m\right) + y\left(n + 1, m\right) y\left(n, m\right) + y\left(n, m + 1\right) y\left(n, m\right)}{4 \, {\rm det}\left(n + 1, m + 1\right) {\rm det}\left(n, m\right)} + \frac{x\left(n + 1, m - 1\right) x\left(n - 1, m\right) - x\left(n + 1, m\right) x\left(n - 1, m\right) - 2 \, x\left(n + 1, m - 1\right) x\left(n, m - 1\right) + x\left(n + 1, m\right) x\left(n, m - 1\right) + x\left(n + 1, m - 1\right) x\left(n, m\right) + x\left(n, m - 1\right) x\left(n, m\right) - x\left(n, m\right)^{2} + y\left(n + 1, m - 1\right) y\left(n - 1, m\right) - y\left(n + 1, m\right) y\left(n - 1, m\right) - 2 \, y\left(n + 1, m - 1\right) y\left(n, m - 1\right) + y\left(n + 1, m\right) y\left(n, m - 1\right) + y\left(n + 1, m - 1\right) y\left(n, m\right) + y\left(n, m - 1\right) y\left(n, m\right) - y\left(n, m\right)^{2}}{4 \, {\rm det}\left(n + 1, m\right) {\rm det}\left(n, m\right)} + \frac{x\left(n - 1, m + 2\right) x\left(n - 1, m\right) - x\left(n - 1, m\right) x\left(n, m + 2\right) - x\left(n - 1, m + 2\right) x\left(n, m\right) + x\left(n, m + 2\right) x\left(n, m\right) + y\left(n - 1, m + 2\right) y\left(n - 1, m\right) - y\left(n - 1, m\right) y\left(n, m + 2\right) - y\left(n - 1, m + 2\right) y\left(n, m\right) + y\left(n, m + 2\right) y\left(n, m\right)}{4 \, {\rm det}\left(n, m + 2\right) {\rm det}\left(n, m\right)} - \frac{2 \, x\left(n - 1, m + 1\right) x\left(n - 1, m\right) - x\left(n - 1, m\right) x\left(n, m + 1\right) - x\left(n - 1, m + 1\right) x\left(n, m - 1\right) + x\left(n, m + 1\right) x\left(n, m - 1\right) - x\left(n - 1, m + 1\right) x\left(n, m\right) - x\left(n - 1, m\right) x\left(n, m\right) + x\left(n, m\right)^{2} + 2 \, y\left(n - 1, m + 1\right) y\left(n - 1, m\right) - y\left(n - 1, m\right) y\left(n, m + 1\right) - y\left(n - 1, m + 1\right) y\left(n, m - 1\right) + y\left(n, m + 1\right) y\left(n, m - 1\right) - y\left(n - 1, m + 1\right) y\left(n, m\right) - y\left(n - 1, m\right) y\left(n, m\right) + y\left(n, m\right)^{2}}{4 \, {\rm det}\left(n, m + 1\right) {\rm det}\left(n, m\right)}
eq20nm2 = sum([ sum([ sum([ (hbar^2/4/mu)*(1/2*Dplus(Mnm[l,j]/det(n,m),l)*Dplus(Mnm[k,j]/det(n,m),k)) for k in range(d)]) for j in range(d)]) for l in range(d)])
show(eq20nm2.subs({hbar:1,mu:1}))
18(x(n+1,m1)x(n+1,m)det(n+1,m)x(n,m1)x(n,m)det(n,m))2+18(y(n+1,m1)y(n+1,m)det(n+1,m)y(n,m1)y(n,m)det(n,m))214(x(n+1,m1)x(n+1,m)det(n+1,m)x(n,m1)x(n,m)det(n,m))(x(n1,m+1)x(n,m+1)det(n,m+1)x(n1,m)x(n,m)det(n,m))+18(x(n1,m+1)x(n,m+1)det(n,m+1)x(n1,m)x(n,m)det(n,m))214(y(n+1,m1)y(n+1,m)det(n+1,m)y(n,m1)y(n,m)det(n,m))(y(n1,m+1)y(n,m+1)det(n,m+1)y(n1,m)y(n,m)det(n,m))+18(y(n1,m+1)y(n,m+1)det(n,m+1)y(n1,m)y(n,m)det(n,m))2\displaystyle \frac{1}{8} \, {\left(\frac{x\left(n + 1, m - 1\right) - x\left(n + 1, m\right)}{{\rm det}\left(n + 1, m\right)} - \frac{x\left(n, m - 1\right) - x\left(n, m\right)}{{\rm det}\left(n, m\right)}\right)}^{2} + \frac{1}{8} \, {\left(\frac{y\left(n + 1, m - 1\right) - y\left(n + 1, m\right)}{{\rm det}\left(n + 1, m\right)} - \frac{y\left(n, m - 1\right) - y\left(n, m\right)}{{\rm det}\left(n, m\right)}\right)}^{2} - \frac{1}{4} \, {\left(\frac{x\left(n + 1, m - 1\right) - x\left(n + 1, m\right)}{{\rm det}\left(n + 1, m\right)} - \frac{x\left(n, m - 1\right) - x\left(n, m\right)}{{\rm det}\left(n, m\right)}\right)} {\left(\frac{x\left(n - 1, m + 1\right) - x\left(n, m + 1\right)}{{\rm det}\left(n, m + 1\right)} - \frac{x\left(n - 1, m\right) - x\left(n, m\right)}{{\rm det}\left(n, m\right)}\right)} + \frac{1}{8} \, {\left(\frac{x\left(n - 1, m + 1\right) - x\left(n, m + 1\right)}{{\rm det}\left(n, m + 1\right)} - \frac{x\left(n - 1, m\right) - x\left(n, m\right)}{{\rm det}\left(n, m\right)}\right)}^{2} - \frac{1}{4} \, {\left(\frac{y\left(n + 1, m - 1\right) - y\left(n + 1, m\right)}{{\rm det}\left(n + 1, m\right)} - \frac{y\left(n, m - 1\right) - y\left(n, m\right)}{{\rm det}\left(n, m\right)}\right)} {\left(\frac{y\left(n - 1, m + 1\right) - y\left(n, m + 1\right)}{{\rm det}\left(n, m + 1\right)} - \frac{y\left(n - 1, m\right) - y\left(n, m\right)}{{\rm det}\left(n, m\right)}\right)} + \frac{1}{8} \, {\left(\frac{y\left(n - 1, m + 1\right) - y\left(n, m + 1\right)}{{\rm det}\left(n, m + 1\right)} - \frac{y\left(n - 1, m\right) - y\left(n, m\right)}{{\rm det}\left(n, m\right)}\right)}^{2}

Manually reorganizing terms, we can put it into another form.

eq20nm4=((x(n+1,m-1)-x(n+1,m))^2+(y(n+1,m-1)-y(n+1,m))^2)/8/det(n+1,m)^2 + ((x(n-1,m+1)-x(n,m+1))^2+(y(n-1,m+1)-y(n,m+1))^2)/8/det(n,m+1)^2 + ((x(n-1,m)-x(n,m-1))^2+(y(n-1,m)-y(n,m-1))^2)/8/det(n,m)^2 - ((x(n-1,m+1)-x(n,m+1))*(x(n-1,m)-x(n,m-1)) + (y(n-1,m+1)-y(n,m+1))*(y(n-1,m)-y(n,m-1)))/4/det(n,m+1)/det(n,m) + ((x(n+1,m-1)-x(n+1,m))*(x(n-1,m)-x(n,m-1)) + (y(n+1,m-1)-y(n+1,m))*(y(n-1,m)-y(n,m-1)))/4/det(n+1,m)/det(n,m) - ((x(n+1,m-1)-x(n+1,m))*(x(n-1,m+1)-x(n,m+1)) + (y(n+1,m-1)-y(n+1,m))*(y(n-1,m+1)-y(n,m+1)))/4/det(n+1,m)/det(n,m+1) show(eq20nm4)
(x(n+1,m1)x(n+1,m))2+(y(n+1,m1)y(n+1,m))28det(n+1,m)2+(x(n1,m+1)x(n,m+1))2+(y(n1,m+1)y(n,m+1))28det(n,m+1)2(x(n+1,m1)x(n+1,m))(x(n1,m+1)x(n,m+1))+(y(n+1,m1)y(n+1,m))(y(n1,m+1)y(n,m+1))4det(n+1,m)det(n,m+1)+(x(n1,m)x(n,m1))2+(y(n1,m)y(n,m1))28det(n,m)2+(x(n+1,m1)x(n+1,m))(x(n1,m)x(n,m1))+(y(n+1,m1)y(n+1,m))(y(n1,m)y(n,m1))4det(n+1,m)det(n,m)(x(n1,m+1)x(n,m+1))(x(n1,m)x(n,m1))+(y(n1,m+1)y(n,m+1))(y(n1,m)y(n,m1))4det(n,m+1)det(n,m)\displaystyle \frac{{\left(x\left(n + 1, m - 1\right) - x\left(n + 1, m\right)\right)}^{2} + {\left(y\left(n + 1, m - 1\right) - y\left(n + 1, m\right)\right)}^{2}}{8 \, {\rm det}\left(n + 1, m\right)^{2}} + \frac{{\left(x\left(n - 1, m + 1\right) - x\left(n, m + 1\right)\right)}^{2} + {\left(y\left(n - 1, m + 1\right) - y\left(n, m + 1\right)\right)}^{2}}{8 \, {\rm det}\left(n, m + 1\right)^{2}} - \frac{{\left(x\left(n + 1, m - 1\right) - x\left(n + 1, m\right)\right)} {\left(x\left(n - 1, m + 1\right) - x\left(n, m + 1\right)\right)} + {\left(y\left(n + 1, m - 1\right) - y\left(n + 1, m\right)\right)} {\left(y\left(n - 1, m + 1\right) - y\left(n, m + 1\right)\right)}}{4 \, {\rm det}\left(n + 1, m\right) {\rm det}\left(n, m + 1\right)} + \frac{{\left(x\left(n - 1, m\right) - x\left(n, m - 1\right)\right)}^{2} + {\left(y\left(n - 1, m\right) - y\left(n, m - 1\right)\right)}^{2}}{8 \, {\rm det}\left(n, m\right)^{2}} + \frac{{\left(x\left(n + 1, m - 1\right) - x\left(n + 1, m\right)\right)} {\left(x\left(n - 1, m\right) - x\left(n, m - 1\right)\right)} + {\left(y\left(n + 1, m - 1\right) - y\left(n + 1, m\right)\right)} {\left(y\left(n - 1, m\right) - y\left(n, m - 1\right)\right)}}{4 \, {\rm det}\left(n + 1, m\right) {\rm det}\left(n, m\right)} - \frac{{\left(x\left(n - 1, m + 1\right) - x\left(n, m + 1\right)\right)} {\left(x\left(n - 1, m\right) - x\left(n, m - 1\right)\right)} + {\left(y\left(n - 1, m + 1\right) - y\left(n, m + 1\right)\right)} {\left(y\left(n - 1, m\right) - y\left(n, m - 1\right)\right)}}{4 \, {\rm det}\left(n, m + 1\right) {\rm det}\left(n, m\right)}
bool((eq20nm2.subs({hbar:1,mu:1})-eq20nm4).full_simplify()==0)
True

But for the computation we don't really need the determinants. (Only use the second term.)

eq20nm0 = sum([ sum([ sum([ (hbar^2/4/mu)*(1/2*Dplus(Knm[l,j],l)*Dplus(Knm[k,j],k)) for k in range(d)]) for j in range(d)]) for l in range(d)]) show(eq20nm0.collect(hbar).collect(mu))
((x(n+1,m1)x(n+1,m)x(n+1,m)y(n+1,m1)x(n,m)y(n+1,m1)x(n+1,m1)y(n+1,m)+x(n,m)y(n+1,m)+x(n+1,m1)y(n,m)x(n+1,m)y(n,m)x(n,m1)x(n,m)x(n,m1)y(n1,m)x(n,m)y(n1,m)x(n1,m)y(n,m1)+x(n,m)y(n,m1)+x(n1,m)y(n,m)x(n,m1)y(n,m))2+2(x(n+1,m1)x(n+1,m)x(n+1,m)y(n+1,m1)x(n,m)y(n+1,m1)x(n+1,m1)y(n+1,m)+x(n,m)y(n+1,m)+x(n+1,m1)y(n,m)x(n+1,m)y(n,m)x(n,m1)x(n,m)x(n,m1)y(n1,m)x(n,m)y(n1,m)x(n1,m)y(n,m1)+x(n,m)y(n,m1)+x(n1,m)y(n,m)x(n,m1)y(n,m))(x(n1,m+1)x(n,m+1)x(n,m+1)y(n1,m+1)x(n,m)y(n1,m+1)x(n1,m+1)y(n,m+1)+x(n,m)y(n,m+1)+x(n1,m+1)y(n,m)x(n,m+1)y(n,m)+x(n1,m)x(n,m)x(n,m1)y(n1,m)x(n,m)y(n1,m)x(n1,m)y(n,m1)+x(n,m)y(n,m1)+x(n1,m)y(n,m)x(n,m1)y(n,m))+(x(n1,m+1)x(n,m+1)x(n,m+1)y(n1,m+1)x(n,m)y(n1,m+1)x(n1,m+1)y(n,m+1)+x(n,m)y(n,m+1)+x(n1,m+1)y(n,m)x(n,m+1)y(n,m)+x(n1,m)x(n,m)x(n,m1)y(n1,m)x(n,m)y(n1,m)x(n1,m)y(n,m1)+x(n,m)y(n,m1)+x(n1,m)y(n,m)x(n,m1)y(n,m))2+(y(n+1,m1)y(n+1,m)x(n+1,m)y(n+1,m1)x(n,m)y(n+1,m1)x(n+1,m1)y(n+1,m)+x(n,m)y(n+1,m)+x(n+1,m1)y(n,m)x(n+1,m)y(n,m)y(n,m1)y(n,m)x(n,m1)y(n1,m)x(n,m)y(n1,m)x(n1,m)y(n,m1)+x(n,m)y(n,m1)+x(n1,m)y(n,m)x(n,m1)y(n,m))2+2(y(n+1,m1)y(n+1,m)x(n+1,m)y(n+1,m1)x(n,m)y(n+1,m1)x(n+1,m1)y(n+1,m)+x(n,m)y(n+1,m)+x(n+1,m1)y(n,m)x(n+1,m)y(n,m)y(n,m1)y(n,m)x(n,m1)y(n1,m)x(n,m)y(n1,m)x(n1,m)y(n,m1)+x(n,m)y(n,m1)+x(n1,m)y(n,m)x(n,m1)y(n,m))(y(n1,m+1)y(n,m+1)x(n,m+1)y(n1,m+1)x(n,m)y(n1,m+1)x(n1,m+1)y(n,m+1)+x(n,m)y(n,m+1)+x(n1,m+1)y(n,m)x(n,m+1)y(n,m)+y(n1,m)y(n,m)x(n,m1)y(n1,m)x(n,m)y(n1,m)x(n1,m)y(n,m1)+x(n,m)y(n,m1)+x(n1,m)y(n,m)x(n,m1)y(n,m))+(y(n1,m+1)y(n,m+1)x(n,m+1)y(n1,m+1)x(n,m)y(n1,m+1)x(n1,m+1)y(n,m+1)+x(n,m)y(n,m+1)+x(n1,m+1)y(n,m)x(n,m+1)y(n,m)+y(n1,m)y(n,m)x(n,m1)y(n1,m)x(n,m)y(n1,m)x(n1,m)y(n,m1)+x(n,m)y(n,m1)+x(n1,m)y(n,m)x(n,m1)y(n,m))2)28μ\displaystyle \frac{{\left({\left(\frac{x\left(n + 1, m - 1\right) - x\left(n + 1, m\right)}{x\left(n + 1, m\right) y\left(n + 1, m - 1\right) - x\left(n, m\right) y\left(n + 1, m - 1\right) - x\left(n + 1, m - 1\right) y\left(n + 1, m\right) + x\left(n, m\right) y\left(n + 1, m\right) + x\left(n + 1, m - 1\right) y\left(n, m\right) - x\left(n + 1, m\right) y\left(n, m\right)} - \frac{x\left(n, m - 1\right) - x\left(n, m\right)}{x\left(n, m - 1\right) y\left(n - 1, m\right) - x\left(n, m\right) y\left(n - 1, m\right) - x\left(n - 1, m\right) y\left(n, m - 1\right) + x\left(n, m\right) y\left(n, m - 1\right) + x\left(n - 1, m\right) y\left(n, m\right) - x\left(n, m - 1\right) y\left(n, m\right)}\right)}^{2} + 2 \, {\left(\frac{x\left(n + 1, m - 1\right) - x\left(n + 1, m\right)}{x\left(n + 1, m\right) y\left(n + 1, m - 1\right) - x\left(n, m\right) y\left(n + 1, m - 1\right) - x\left(n + 1, m - 1\right) y\left(n + 1, m\right) + x\left(n, m\right) y\left(n + 1, m\right) + x\left(n + 1, m - 1\right) y\left(n, m\right) - x\left(n + 1, m\right) y\left(n, m\right)} - \frac{x\left(n, m - 1\right) - x\left(n, m\right)}{x\left(n, m - 1\right) y\left(n - 1, m\right) - x\left(n, m\right) y\left(n - 1, m\right) - x\left(n - 1, m\right) y\left(n, m - 1\right) + x\left(n, m\right) y\left(n, m - 1\right) + x\left(n - 1, m\right) y\left(n, m\right) - x\left(n, m - 1\right) y\left(n, m\right)}\right)} {\left(\frac{x\left(n - 1, m + 1\right) - x\left(n, m + 1\right)}{x\left(n, m + 1\right) y\left(n - 1, m + 1\right) - x\left(n, m\right) y\left(n - 1, m + 1\right) - x\left(n - 1, m + 1\right) y\left(n, m + 1\right) + x\left(n, m\right) y\left(n, m + 1\right) + x\left(n - 1, m + 1\right) y\left(n, m\right) - x\left(n, m + 1\right) y\left(n, m\right)} + \frac{x\left(n - 1, m\right) - x\left(n, m\right)}{x\left(n, m - 1\right) y\left(n - 1, m\right) - x\left(n, m\right) y\left(n - 1, m\right) - x\left(n - 1, m\right) y\left(n, m - 1\right) + x\left(n, m\right) y\left(n, m - 1\right) + x\left(n - 1, m\right) y\left(n, m\right) - x\left(n, m - 1\right) y\left(n, m\right)}\right)} + {\left(\frac{x\left(n - 1, m + 1\right) - x\left(n, m + 1\right)}{x\left(n, m + 1\right) y\left(n - 1, m + 1\right) - x\left(n, m\right) y\left(n - 1, m + 1\right) - x\left(n - 1, m + 1\right) y\left(n, m + 1\right) + x\left(n, m\right) y\left(n, m + 1\right) + x\left(n - 1, m + 1\right) y\left(n, m\right) - x\left(n, m + 1\right) y\left(n, m\right)} + \frac{x\left(n - 1, m\right) - x\left(n, m\right)}{x\left(n, m - 1\right) y\left(n - 1, m\right) - x\left(n, m\right) y\left(n - 1, m\right) - x\left(n - 1, m\right) y\left(n, m - 1\right) + x\left(n, m\right) y\left(n, m - 1\right) + x\left(n - 1, m\right) y\left(n, m\right) - x\left(n, m - 1\right) y\left(n, m\right)}\right)}^{2} + {\left(\frac{y\left(n + 1, m - 1\right) - y\left(n + 1, m\right)}{x\left(n + 1, m\right) y\left(n + 1, m - 1\right) - x\left(n, m\right) y\left(n + 1, m - 1\right) - x\left(n + 1, m - 1\right) y\left(n + 1, m\right) + x\left(n, m\right) y\left(n + 1, m\right) + x\left(n + 1, m - 1\right) y\left(n, m\right) - x\left(n + 1, m\right) y\left(n, m\right)} - \frac{y\left(n, m - 1\right) - y\left(n, m\right)}{x\left(n, m - 1\right) y\left(n - 1, m\right) - x\left(n, m\right) y\left(n - 1, m\right) - x\left(n - 1, m\right) y\left(n, m - 1\right) + x\left(n, m\right) y\left(n, m - 1\right) + x\left(n - 1, m\right) y\left(n, m\right) - x\left(n, m - 1\right) y\left(n, m\right)}\right)}^{2} + 2 \, {\left(\frac{y\left(n + 1, m - 1\right) - y\left(n + 1, m\right)}{x\left(n + 1, m\right) y\left(n + 1, m - 1\right) - x\left(n, m\right) y\left(n + 1, m - 1\right) - x\left(n + 1, m - 1\right) y\left(n + 1, m\right) + x\left(n, m\right) y\left(n + 1, m\right) + x\left(n + 1, m - 1\right) y\left(n, m\right) - x\left(n + 1, m\right) y\left(n, m\right)} - \frac{y\left(n, m - 1\right) - y\left(n, m\right)}{x\left(n, m - 1\right) y\left(n - 1, m\right) - x\left(n, m\right) y\left(n - 1, m\right) - x\left(n - 1, m\right) y\left(n, m - 1\right) + x\left(n, m\right) y\left(n, m - 1\right) + x\left(n - 1, m\right) y\left(n, m\right) - x\left(n, m - 1\right) y\left(n, m\right)}\right)} {\left(\frac{y\left(n - 1, m + 1\right) - y\left(n, m + 1\right)}{x\left(n, m + 1\right) y\left(n - 1, m + 1\right) - x\left(n, m\right) y\left(n - 1, m + 1\right) - x\left(n - 1, m + 1\right) y\left(n, m + 1\right) + x\left(n, m\right) y\left(n, m + 1\right) + x\left(n - 1, m + 1\right) y\left(n, m\right) - x\left(n, m + 1\right) y\left(n, m\right)} + \frac{y\left(n - 1, m\right) - y\left(n, m\right)}{x\left(n, m - 1\right) y\left(n - 1, m\right) - x\left(n, m\right) y\left(n - 1, m\right) - x\left(n - 1, m\right) y\left(n, m - 1\right) + x\left(n, m\right) y\left(n, m - 1\right) + x\left(n - 1, m\right) y\left(n, m\right) - x\left(n, m - 1\right) y\left(n, m\right)}\right)} + {\left(\frac{y\left(n - 1, m + 1\right) - y\left(n, m + 1\right)}{x\left(n, m + 1\right) y\left(n - 1, m + 1\right) - x\left(n, m\right) y\left(n - 1, m + 1\right) - x\left(n - 1, m + 1\right) y\left(n, m + 1\right) + x\left(n, m\right) y\left(n, m + 1\right) + x\left(n - 1, m + 1\right) y\left(n, m\right) - x\left(n, m + 1\right) y\left(n, m\right)} + \frac{y\left(n - 1, m\right) - y\left(n, m\right)}{x\left(n, m - 1\right) y\left(n - 1, m\right) - x\left(n, m\right) y\left(n - 1, m\right) - x\left(n - 1, m\right) y\left(n, m - 1\right) + x\left(n, m\right) y\left(n, m - 1\right) + x\left(n - 1, m\right) y\left(n, m\right) - x\left(n, m - 1\right) y\left(n, m\right)}\right)}^{2}\right)} {{\hbar}}^{2}}{8 \, {{\mu}}}

Numeric form.

eq20nm3=eq20nm0.subs(dict(zip(xyk,xyv))) ru=fast_callable(eq20nm3.subs({hbar:0.2,mu:1.0}), vars=xyv, domain=RDF)

Classical force

F0(x)= 9*sech(5*x)^2*tanh(5*x)
plot(F0,(-3,3),axes_labels=['distance','force'],axes_labels_size=1)

Classical potential

P0(x)=9/(10*cosh(5*x)^2)
plot(P0,(-3,3),axes_labels=['distance','energy'],axes_labels_size=1,legend_label='potential')

Numeric form.

F1=fast_callable(F0(x11), vars=[x11], domain=RDF) def FF(xy): FF1=F1(xy[0]-xy[1]); return array([FF1,-FF1]) def PO(xy): return P0(xy[0]-xy[1])

Initial Conditions

Consider two particles in a one dimensional space over 4 worlds. The initial configuration space is shown below.

N=2; M=2 r0=array([[[(i-0.5)/N,(j-0.5)/M] for j in range(-M/2+1,M/2+1)] for i in range(-N/2+1,N/2+1)]) p0=array(map(lambda y:map(lambda x:array([1.0*RDF(erfinv(2.0*x[0])-2.0+x[1]/2),RDF(erfinv(2.0*x[1])+2.0+x[0]/2)]),y),r0)) list_plot([r0[i,j] for i in range(N) for j in range(M)])+list_plot([p0[i,j] for i in range(N) for j in range(M)],color='red',figsize=[5,5])

The particles approach each other with the same initial velocity in each world.

v0=array([[[RDF(1.0),RDF(-1.0)] for j in range(M)] for i in range(N)])

In order to compute the quantum force on the particles in world n,mn,m we introduce ficticious worlds at "approximate infinity" to represent open boundry conditions.

inf = RDF(10.0^12); inf2 = RDF(10.0^6); def rf1(XY,n,m): global N,M p=[[XY[n+i,m+j][0],XY[n+i,m+j][1]] if n+i>=0 and m+j>=0 and n+i<N and m+j<M else [ -inf if n+i<-1 else -inf2 if n+i<=-1 else inf if n+i>N else inf2 if n+i>=N else XY[n+i,m][0], -inf if m+j<-1 else -inf2 if m+j<=-1 else inf if m+j>M else inf2 if m+j>=M else XY[n,m+j][1]] for i in [-2,-1,0,1,2] for j in [-2,-1,0,1,2]] p=reduce(lambda x,y:x+y,p) r=ff(p) if isnan(r[0]): print "r[0] isnan" if isnan(r[1]): print "r[1] isnan" return r

Total classical potential

def PU(x): return sum([PO(x[i,j]) for i in range(N) for j in range(M)])

Total quantum Potential

def ru1(XY,n,m): global N,M nm=n+m+1 p=[[XY[n+i,m+j][0],XY[n+i,m+j][1]] if n+i>=0 and m+j>=0 and n+i<N and m+j<M else [ -inf if n+i<-1 else -inf2 if n+i<=-1 else inf if n+i>N else inf2 if n+i>=N else XY[n+i,m][0], -inf if m+j<-1 else -inf2 if m+j<=-1 else inf if m+j>M else inf2 if m+j>=M else XY[n,m+j][1]] for i in [-2,-1,0,1,2] for j in [-2,-1,0,1,2]] p=reduce(lambda x,y:x+y,p) return ru(*p) def QU(x): return sum([ru1(x,i,j) for i in range(N) for j in range(M)])

Total kinetic energy

def K(x): return RDF(0.5)*sum([sum(x[i,j]^2) for i in range(N) for j in range(M)])

Simulation

Run the integration using Störmer–Verlet discretization with a continuous integrating stepsize controller.

%time t_end = 6.0; step = .04 # initial step size ds = RDF(0.001) # controller alpha = RDF(30.0); beta = RDF(0.0) # step density rho = RDF(1.0); dt=ds/rho def G(p,q): a = sum([sum(p[i,j]*q[i,j]) for j in range(M) for i in range(N)]) b = sum([sum(q[i,j]*q[i,j]) for j in range(M) for i in range(N)]) eps = RDF(10.0^-10) return -alpha*a/max(b,eps)+beta*(1.0-ds/dt) #acceleration = (quantum force + classical force)/mass mm = RDF(1.0) def A(x): return RDF(0.5)*array([[(rf1(x,i,j)+FF(XY[i,j]))/mm for j in range(M)] for i in range(N)]) # initial values XY = p0.copy(); VV = v0.copy(); AX = A(XY) rho = rho - RDF(0.5)*G(AX,VV)*dt t = 0; tt = [t] XV = [ [ [ [XY[i,j][0],XY[i,j][1]] for i in range(N)] for j in range(M)] ] KK=K(VV); PP=PU(XY); QQ=QU(XY) T1=KK+PP+QQ rhos = [rho]; TK={t:KK}; TP={t:PP};TQ={t:QQ}; T=[T1] while t<t_end: t1 = t+step while t<t1: rho = rho + G(AX,VV)*dt if isnan(rho) or rho>100 or rho<0.01: print "Step control failed at %s. rho=%s"%(t,rho) t=t_end; break dt = ds/rho; t += dt # note 0.5 in A above XY += VV*dt + AX*dt^2 VV += AX*dt; AX = A(XY); VV += AX*dt # Check total energy conservation KK=K(VV); PP=PU(XY); QQ=QU(XY) T2=KK+PP+QQ if abs(T1-T2)>0.1: print "Integration failed at %s. Try a shorter time step: %s."%(t,ds/2) #break T1=T2 # record XV = XV + [ [ [ [XY[i,j][0],XY[i,j][1]] for i in range(N)] for j in range(M)] ]; tt += [t] rhos.append(rho); TK[t]=KK; TP[t]=PP; TQ[t]=QQ; T.append(T1)
CPU time: 33.48 s, Wall time: 41.55 s
line2d(zip(tt,rhos),axes_labels=['time','rho'],legend_label='step density') (line2d([[t,TK[t]] for t in tt],axes_labels=['time','energy'],legend_label='kinetic',color='green')+ line2d([[t,TQ[t]] for t in tt],axes_labels=['time','energy'],legend_label='quantum potential',color='red')+ line2d([[t,TP[t]] for t in tt],axes_labels=['time','energy'],legend_label='classical potential',color='blue')) line2d(zip(tt,T),axes_labels=['time','energy'],legend_label='total',color='black')
line2d([[t,TQ[t]] for t in tt],axes_labels=['time','energy'],legend_label='quantum potential',color='red')

You might want something like this: https://addons.mozilla.org/En-US/firefox/addon/toggle-animated-gifs/ to manage the animation.

Start/stop GIF animations through a keyboard shortcut or by clicking them. You can also restart animations from the beginning, or disable animations by default. %md Each point in the following diagram represents the position of two particles. When the point approaches the x=y line

Each moving point in the following diagram represents the position of two particles in a single world (four worlds are shown). As the point approaches the x=y line (shown in red), the particles approach each other and the classical potential is at a maximum. In three worlds the particles overcome the potential barrier and continue. But one world (the last to approach) is reflected back. Although initially the two particles in each world approach each other with the same velocity, the particles in the last world to approach the barrier have lost energy to their counterparts in other worlds due to the quantum mechanical (interworld) interaction. When viewed in the projection we see two particles rebound instead of passing each other.

animate([text("time:"+("%4.1f"%x[0]).rjust(6), (6,3),color='black') + text("kinetic:"+("%.4e"%TK[x[0]]).rjust(10), (6,2.5),color='black') + text("classic:"+("%.4e"%TP[x[0]]).rjust(10), (6,2),color='black') + text("quantum:"+("%.4e"%TQ[x[0]]).rjust(10), (6,1.5),color='black') + line([[-6,-6],[6,6]],color='red',legend_label='potential barrier') + list_plot([x[1][i][j] for j in range(N) for i in range(M)],xmin=-6,xmax=6,ymin=-6,ymax=6) for x in zip(tt,XV)]).show(gif=true,delay=20)
show(sum([line([[tt[i],XV[i][j][k][0]] for i in range(0,len(tt))],color = 'red') for j in range(M) for k in range(N)])+ sum([line([[tt[i],XV[i][j][k][1]] for i in range(0,len(tt))],color = 'blue') for j in range(M) for k in range(N)]) ,figsize=[10,7],axes_labels=['time','location'],axes_labels_size=1)
animate([text("time:"+("%4.1f"%x[0]).rjust(5), (3,3),color='black') + list_plot([[x[1][i][j][0],PO(x[1][i][j])] for j in range(N) for i in range(M)],xmin=-6,xmax=6,ymin=0,ymax=3,color='red')+ list_plot([[x[1][i][j][1],PO(x[1][i][j])] for j in range(N) for i in range(M)],xmin=-6,xmax=6,ymin=0,ymax=3,color='blue') for x in zip(tt,XV)]).show(gif=true,delay=20)