Bohmian particles at a rectangular step

Julian Grond and David Kammerlander

October 2004

Load Packages

<< VQM`;

Off[General::spell1];

Parameters and Definitions

Parameters for the numerical computation

Make the numerical calculations on a sufficiently large bounded domain. The wavefunction is set equal to zero at the border.

numleft = -100.;      (* left border of the domain *)
numright = 100.;      (* right border of the domain *)
dx = Pi/20.;           (* space step *)
dt = 0.05;          (* time step *)

Definition of the wave packet

p0 = 2;		(* average momentum of incoming wave packet *)
x0 = 0;		(* average position at t=0 *)
t0 = -25;	(* initial time *)

This prepares a table of numerical values:

psi0 = N[Table[ QFreeGaussian1D[x,t0,x0,p0,0.01,1], {x,numleft,numright,dx}]];

We can visualize this list of complex numbers using the command QListArgColorPlot provided by the VQM` packages

QListArgColorPlot[psi0,  PlotRange ->  All, 
          QHorizontalRange -> {{-100, 100}, {-90, 30}}];

Definition of the potential

V0 = 1.6;		(* height of the potential step *)

potfunc[x_] = If[x>0, V0, 0];

This prepares a table of numerical values:

pot = N[Table[ potfunc[x], {x,numleft,numright,dx}]];

Initialization

This gives the initial wavefunction, the potential operator and the Hamilton operator.

psi = QNewFunction[Re[psi0],Im[psi0]];

V   = QNewFunction[Re[pot], Im[pot]];

H   = QSchroedinger1D[V, 1, dx];

Velocity Field

We will numerically compute the velocity field, given by the guiding equation: v =(ℏ/m) Im Grad[ψ]/ψ, with ℏ=m=1. A list will be constructed with each line the discretized (dx)  velocity field at a certain time (step dT = 10 dt, because of  QTimeEvolution[H,psi,dt,0,10]). We use an accuracy of 10 digits after the point and average out neighbouring values of the list. When a value in the list is too small, no calculation is done and set to zero.

a = 10 ;

psicom = QGetArray[psi][[1]] + 
          I QGetArray[psi] [[2]];

psicomAv = 
        SetAccuracy[(1/2) (psicom + RotateLeft[psicom]), a];

v0 = 
        Table[If[Abs[psicomAv[[n]]] > 0.01, 
            ComplexExpand[
              Im[((psicomAv[[n + 1]] - psicomAv[[n]])/((1/2) dx (psicomAv[[n]] + psicomAv[[n + 1]])))]], 0], {n, Length[psicom] - 1}];

v1 = 
        Table[Do[QTimeEvolution[H, psi, dt, 0, 10]; 
            psicomAv = 
              (1/2) SetAccuracy[QGetArray[
                        psi][[1]] + I QGetArray[psi][[2]] + 
                    RotateLeft[QGetArray[
                            psi][[1]] + I QGetArray[
                              psi][[2]]], a]; Return[
              Table[If[
                  Abs[psicomAv[[n]]] > 
                    0.01, ComplexExpand[
                    Im[((psicomAv[[n + 1]] - psicomAv[[n]])/((1/2) dx (psicomAv[[n]] + psicomAv[[n + 1]])))]], 0], {n, 1, Length[psicom] - 1}]], {q, 
              1}], {nn, 1, (60/(10 dt))}];

vv = Join[{v0}, v1];

Plot of the velocity field in the region where the wave packet is not too small:

Table[
        ListPlot[vv[[n]],  
          PlotLabel -> "= t" ((n - 1) 10 dt - 25), 
          AxesLabel  ->  {x, v}, PlotRange -> {-5, 5}, 
          Ticks -> {{ {200, SetAccuracy[200 dx - 100, 1]}, {400, 
                  SetAccuracy[400 dx - 100, 1]}, {600, 
                  SetAccuracy[600 dx - 100, 1]}, {800, 
                  SetAccuracy[800 dx - 100, 1]}, {1000, 
                  SetAccuracy[1000 dx - 100, 1]}, {1200, 
                  SetAccuracy[1200 dx - 100, 1]}}, Automatic}], {n, 
          60/(10 dt)}];

Motion

The motion of a particle is given by (dQ/dt) = v(Q) (we use also Q as the initial coordinate of a particle). Runge - Kutta 4 is used.

ClearAll[s1, s2, s3, s4, Q]; dT = 10 dt;

s1[Q_, t_] := 
        dT vv[[
              Floor[t/dT]]][[
            Floor[((Q+100 )/dx)]]];
s2[Q_, t_] := 
        dT vv[[
              Floor[ (t + dT/2)/
                  dT]]][[
            Floor[(((Q + s1[Q, t]/2 + 100) )/dx)]]] ;
s3[Q_, t_] := 
      dT  vv[[
            Floor[(t + dT/2)/
                dT]]][[
          Floor[(((Q + s2[Q, t]/2 + 100) )/dx)]]]  ; 
    s4[Q_, t_] := 
      dT vv[[Floor[(t + dT)/
                dT]]][[
          Floor[(((Q + s3[Q, t] + 100) )/dx)]]];

Plot of the position of a particle with Gaussian distributed initial values versus time. 40/d+1 gives the number of particles to be plotted.

d = 1; 
    QMotion = 
      Table[Join[
          Q = ((m + 50)/Abs)[QFreeGaussian1D[m, t0, x0, p0, 0.01, 1]]/
                28 - 50; {Q}, 
          Table[Do[
              Return[Q = 
                  Q + (1/6) (dT vv[[
                                Floor[
                                  n ]]][[
                              Floor[((Q+100 )/dx)]]] + 2 dT vv[[
                                Floor[ (n  + 
                                      (1/2))]]][[Floor[(Q + (1/2) s1[Q, n dT] + 100) /
                                  dx]]] + 
                          2 dT vv[[
                                Floor[(n  + 
                                      (1/2))]]][[
                              Floor[(((Q + (1/2) s2[Q, n dT] + 100) )/dx)]]] + 
                          dT vv[[
                                Floor[(n  + 
                                      1)]]][[
                              Floor[(((Q + s3[Q, n dT] + 100) )/dx)]]])], {1}], {n, 
              60/(10 dt)}]], {m, -70, -30, d}];

Table[
        ListPlot[QMotion[[n]],  
          AxesLabel  ->  {prop t,  "x"},  
          PlotRange -> {-80, 50}, 
          PlotLabel ->  "=Q,Particle motion" 
              Evaluate[((n - 70 + 50 - 1)/Abs)[
                          QFreeGaussian1D[n - 70 - 1, t0, x0, p0, 0.01, 1]]/
                    28 - 50]], {n, 40/d + 1}];

Movies of Bohmian Particles

The vertical position of a particle corresponds to the Gaussian initial probability density and the particle density to the probability density at the time.

StepGraph := 
      Line[{{-100, 0}, {0, 0}, {0, V0/4 }, {100, V0/4}}]

\(\(\(Do[
      Plot[potfunc[x]/4, {x, \(-100\), 100}, 
        AxesLabel -> \ {x, "V"}, 
        Ticks -> \ {Automatic, {{0.4, V0}}}, 
        PlotRange -> \ {{\(-100\), 100}, {0, 1}}, 
        Background -> 
          GrayLevel[
            0.1], 
Epilog -> 
{Table[{PointSize[0.02], Hue[m\ \((d + 1)\)/150], 
                Point[{\(QMotion[\([m]\)]\)[[
                      n]], 
                    10  Abs[QFreeGaussian1D[\(\(m - 70 + 50 - 1\)\/Abs[QFreeGaussian1D[m - 70 - 1, t0, x0, p0, 0.01, 1]]\)/28 - 50, t0, x0, p0, 0.01, 1]]\^2
}]}, {m, 1, 
                40/d + 1}]}
, 
        PlotLabel\  -> \ "Motion of the Bohmian Particles - Step 
 
"
], {n, 1, 60/\((10\ dt)\)}]\)\(\ \)\)\)