Particle in a moving wall

Quantum Mechanics
 K. Konishi, G. Paffuti

Problem

Packages and utilities

Off[General::"spell",General::"spell1",InterpolatingFunction::"dmval"];
Off[ClearAll::"clloc"];
ClearAll["Global`*"];
Needs["PlotLegends`"];

Short review of the theory

The system is studied in  the notebook Problems[*] of this chapter, the reader can find in that file some more details on the solution.

A free particle can move in a box [(Subscript[x, L]), (Subscript[x, R])], both ends can in principle move with velocities (Subscript[v, L]) and (Subscript[v, R])respectively.  L = (Subscript[x, R])-(Subscript[x, L]) is the length of the box. The Schrödinger equation is simple but has moving boundary conditions

i ℏ(∂/(∂t))ψ =- ((ℏ^2)/(2m)) ((∂^2)/(∂(x^2)))ψ  ;      ψ[(Subscript[x, L])] =  ψ[(Subscript[x, R])] = 0.

In [*] it has been shown that with the obvious change of variables

ξ = (1/L)(x-(Subscript[x, L]));

eq.() can be transformed in a usual fixed boundary problem:

i ℏ((∂ϕ)/(∂t))=- ((ℏ^2)/(2m (L^2))) (((∂^2)ϕ)/(∂(ξ^2))) +i ℏ (( ((OverDot[L])/L)ξ+((Subscript[v, L])/L))((∂ϕ)/(∂ξ)) + (1/2)((OverDot[L])/L)ϕ) .

and the original wave function is recovered as

ψ[t,x]= (1/(Sqrt[L]))ϕ[t, ξ[t,x]].

The first technique that can be applied to the problem is just a strightforward solution with NDSolve of equation ().

A second possibility is to expand the wave function in terms of "adiabatic modes" i.e. in terms of eigenfunctions of H if boundaries would be fixed, a sort of instantaneous picture of the system. The modes are usual trigonometric functions for a free particle, bounded between (Subscript[x, L])and (Subscript[x, R])

(Subscript[u, k])[x;L] = (Sqrt[(2/L)])Sin[((π  k (x-xL))/L)];   (Subscript[E, k]) = ((ℏ^2)/(2m))(π^2)((k^2)/(L^2)).

and the expansion reads

ψ[x,t] = (∑_k)(Subscript[a, k])[t] (Subscript[u, k])[x;L]

The coefficients (Subscript[a, k])satisfy an equation which is obtained substituting () in ():

((Subscript[da, k])/dt)+ (∑_s) (Subscript[F, ks])(Subscript[a, s])= -i (Subscript[E, k])(Subscript[a, k]) ;

where

(Subscript[F, ks])= (∫_0^L)(Subscript[u, k])(∂_t)(Subscript[u, s]) =  ((2k s)/((k^2)-(s^2)))(((OverDot[L])/L) (((-1))^(k-s))  + (vL/L)(-1+ (((-1))^(k-s))))(k≠ s) ; (Subscript[F, kk])=0.

A second possibility is then to consider a truncated version of () and solve a system of ordinary coupled differential equations in time.

The two methods are not so different because Mathematica  internal procedure for NDSolve is often based on a basis expansion, the basis () can be particularly well suited for this particular problem, especially for adiabatic motions, when the basis become exact.

In the folowing we apply both methods and wil compute several physical quantities, as energy and momentum variations, transition for periodic excitation etc.Theoretical formulas for various quantities will be given when the coputations will be performed.

From now on we will use ℏ = m = 1, and the initial length of the box will be taken 1, L[0] = 1. The ξ variable will be called x for simplicity.

NDSolve

General solution for constant velocities

We have just to write in Mathematica   code equation ()

mySolve1[vL_,vR_,Tmax_,ψIn_]:= Module[{L,eq},
L[t_]:= 1 + (vR - vL)*t;
eq = { I (∂_t) f[x,t] + (1/(2 ((L[t])^2)))(∂_(x,x))f[x,t] -I (((vR-vL)/(L[t])) x + (vL/(L[t])))(∂_x) f[x,t]-
 I ((vR-vL)/(2 L[t]))f[x,t]== 0,f[x,0]== ψIn[x],f[0,t]== 0,f[1,t]== 0};
 (f/.First[NDSolve[eq,f,{x,0,1},{t,0,Tmax},
PrecisionGoal-> 5,
Method->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", DifferenceOrder->"Pseudospectral" }}]])
]

We choose as a method of solution the option {"TensorProductGrid", DifferenceOrder -> "Pseudospectral" }, which is based on a sort of Fourier expansion, this method seems to offer better stability than others, but the reader can explore different possibilities.

To present the results we have different possibilities. The first one is to select one or more velocities and show the variation of |ψ| after some time. We first give these kind of graphs. The time T is choosen in such a way that the length of the box is changed of L/2. The left boundary is held fixed. Initial state is choosen as the fundamental level

vRun =Select[ Range[-5,4,1],#≠0&];
timeRun = 0.5/Abs[vRun];

ψIn[x_]:= (Sqrt[2])Sin[π  x]

Almost all the folowing intructions are for beautyfy the graphs, the only real computation is done with the assignement of sol= mySolve1[], which effectively solve differential equation.

GraphicsGrid[Partition[Table[
v = vRun[[k]];
L[t_]:= 1 + v*t;
Tmax = timeRun[[k]];
Lb = L[Tmax];
deltaL = 0.02Max[1,Lb];
string1 = "v="<>ToString[vRun[[k]]];
thisLabel = string1;
thisArrow = If[v > 0,{
RGBColor[0,0,1],Thickness[.015],Arrowheads[0.07],Arrow[{{1+0.05,1},{1+0.3,1}}]},
{RGBColor[0,0,1],Thickness[.012],Arrowheads[0.07],Arrow[{{1-0.02,1},{1-0.15,1}}]}];
sol=mySolve1[0,v,Tmax,ψIn];
Show[
Plot[Abs[ψIn[x]],{x,0,1},
PlotRange-> All,
Axes-> False,PlotStyle-> {RGBColor[1,0,0],Dashing[{0.02,0.02}]},
PlotLabel-> thisLabel,
Epilog->  {Line[{{-deltaL,0},{2,0}}],
{GrayLevel[0.8],Rectangle[{1,0},{1+deltaL,1.8}]},
RGBColor[0,0,0],Rectangle[{-deltaL,0},{0,1.8}],
RGBColor[1,0,0],Rectangle[{Lb,0},{Lb+deltaL,1.8}],
thisArrow}
],
Plot[1/(Sqrt[L[Tmax]])Abs[sol[x/L[Tmax],Tmax]],{x,0,L[Tmax]}]],
{k,1,Length[vRun]}],3],ImageSize-> 72*8]

The reader must pay atention to possibile instabilities for  v< -5 (high velocity compressions). In this case he have to play with some option of NDSolve, as MaxPoints etc. (see Mathematica help).

A second possibility is to construct an interactive procedure, with Manipulate, in such a way that the user can select the  velocity and see an animation of the motion.

functTmax[vL_,vR_]:= Module[{v,tL,tR,tCrash},
v= vR-vL;
tL =If[vL≥  0,10^5,-1/vL];
tR =If[vR≤ 0,10^5,1/vR];
tCrash =If[v≥ 0,10^5,-0.5/v];
If[ vL== vR == 0,1, N[Min[{tL,tR,tCrash}]]]]

We choose to divide the time interval in 30 frames and use ListAnimate to plot the sequence. Snapshot button open a window with controls for each variable.

nFramesNDSolve = 30;

Manipulate[Module[{L,ψIn,sol,t,dt,temp,TmaxPlot1,frame,deltaL,evolutionPlot},
L[t_]:= 1 + (vR - vL)*t;
TmaxPlot1=functTmax[vL,vR];
dt = TmaxPlot1/nFramesNDSolve;
ψIn[x_]:= (Sqrt[2.])Sin[π  x];
sol = mySolve1[vL,vR,TmaxPlot1,ψIn];
deltaL = 0.02;
evolutionPlot=
Table[t = frame*dt;
Plot[Abs[(1/(Sqrt[L[t]]))sol[((x-vL*t)/(L[t])),t]],{x,vL frame*dt,1+vR frame*dt},PlotRange-> {{-1,2},{0,2.5}},
Epilog-> {RGBColor[0,0,0],Rectangle[{-deltaL+vL t,0},{vL t,2}],
RGBColor[1,0,0],Rectangle[{1+vR t,0},{1+vR t+deltaL,2}]}],{frame,0,nFramesNDSolve}];
ListAnimate[evolutionPlot]
],
{{vL,0},-5,5,0.5,Appearance->"Labeled"},
{{vR,-1},-5,5,0.5,Appearance->"Labeled"},
TrackedSymbols:> {vL,vR}]

The reader can notice that for (Subscript[v, L])= (Subscript[v, R]) the graph is not invariant, as naively expected from Galilei invariance but oscillates. The oscillations must be there, is not a bad approximation to the solution, the reader can explain why?

In notebook [*] it has been shown that the change of mean momentum P and of energy, can be expressed by the value of Ψ derivative at the border:

(dE/dt)=  (( (ℏ^2))/(2m)){(Subscript[v, L])((|ψ'[(Subscript[x, L])]|)^2) - (Subscript[v, R])((|ψ'[(Subscript[x, R])]|)^2)  };    (dP/dt) = (( (ℏ^2))/(2m)){((|ψ'[(Subscript[x, L])]|)^2) -((|ψ'[(Subscript[x, R])]|)^2)  }

Let us consider the particular case of fixed left boundary, (Subscript[v, L])= 0. We compute the force up to a time Tmax necessary to change L of 50%. We plot the whole time interval and the transient, which is very widely oscilating for high velocities. To save space we use a Manipulate panel, the reader can check various velocities.

vRun =Union[Select[ Range[-3,3,1],#≠0&],Select[ Range[-0.5,0.5,0.1],#≠0&]];

Manipulate[Module[{plot1,plot2,Tmax,sol},
Tmax = 1/Abs[2 v];
sol = mySolve1[0,v,Tmax,(Sqrt[2])Sin[#1 π]&];
plot1=Plot[(1/(2 (1+v*t)))(Abs[Derivative[1,0][ sol][0,t]]^2- Abs[Derivative[1,0][ sol][1,t]]^2),{t,0,Tmax},PlotRange-> {{0,Tmax},All},AxesLabel-> {"t","F"}];
plot2=Plot[(1/(2 (1+v*t)))(Abs[Derivative[1,0][ sol][0,t]]^2- Abs[Derivative[1,0][ sol][1,t]]^2),{t,0,0.0025},PlotRange-> {{0,0.0025},All},AxesLabel-> {"t","F"}];
GraphicsGrid[{{Show[plot1],Show[plot2]}},ImageSize-> 72 6]],{{v,0.5},vRun},
TrackedSymbols:> {v}]

A check

A check for the numerical computation and/or our analytical result can be done by computing both side of equation () and comparing the results.

For given velocities the energy at a time t is computed by performing the integral of kinetic term (we write this in a symmetric form):

E[t] = (∫_(Subscript[x, L])^(Subscript[x, R]))𝕕 x (1/L)((|(∂_x)ψ[((x-(Subscript[x, L]))/L),t]|)^2) = (1/(L^2))(∫_0^1)dξ ((|(∂_ξ)ψ[ξ,t]|)^2).

Then can use directly the solution of () to compute E[t]. It is easier to compute a list of values and then constuct an interpolating function.

Let us fix a value for the velocities and compute the solution

vL1 = -0.1; vR1=1.; 
Tmax1=0.9functTmax[vL1,vR1];
L1[t_]:= 1 + (vR1 - vL1) t;
sol1= mySolve1[vL1,vR1,Tmax1,(Sqrt[2])Sin[π #1]&];

We make a list of times and compute energy end momentum

getEnergyMomentum[sol1_]:= Module[{interpEne,interpP,nptTime,listT,listEne,listP},
nptTime = 100; 
listT = Range[0,Tmax1,Tmax1/nptTime];
listEne=Table[ (1/(2((L1[t])^2)))NIntegrate[((Abs[(∂_x) sol1[x,t]])^2),{x,0,1},PrecisionGoal-> 3,AccuracyGoal-> 3],{t,listT}];
interpEne = Interpolation[Transpose[{listT,listEne}]];
listP = Table[ (1/(L1[t]))NIntegrate[Im[Conjugate[sol1[x,t]](∂_x) sol1[x,t]],{x,0,1},
PrecisionGoal-> 3,AccuracyGoal-> 3],{t,listT}];
interpP =  Interpolation[Transpose[{listT,listP}]];
{interpEne,interpP}]

{interpEne,interpP} = getEnergyMomentum[sol1];

For the r.h.s of  () we use the identities

(∂_x)ψ[(Subscript[x, L])]= (1/L)ϕ'[0];             (∂_x)ψ[(Subscript[x, R])]= (1/L)ϕ'[1];

dEdt[t_]:= (1/(2((L1[t])^3)))(vL1*(((Abs[(∂_x) sol1[x,t]])^2)/.x-> 0)-vR1*(((Abs[(∂_x) sol1[x,t]])^2)/.x-> 1));
dPdt[t_]:= (1/(2((L1[t])^3)))((((Abs[(∂_x) sol1[x,t]])^2)/.x-> 0)-(((Abs[(∂_x) sol1[x,t]])^2)/.x-> 1))

GraphicsGrid[{{Plot[{interpEne'[t],dEdt[t]},
{t,0,Tmax1},AxesLabel-> {"t","dE/dt"},
PlotStyle-> {Black,{Red,Dashed}}],Plot[{interpP'[t],dPdt[t]},
{t,0,Tmax1},AxesLabel-> {"t","dP/dt"},PlotStyle-> {Black,{Red,Dashed}}]}},
PlotLabel-> "{vL,vR}= "<>ToString[{vL1,vR1}],ImageSize-> 72 8]

The agreement is acceptable, the reader can check  other velocities.

Adiabatic expansion

In previous section we have checked that energy and momentum variations agree with theoretical expectations. In particular we stress that eq.() predicts dE/dt < 0 for an expansion (as an instance vL = 0 and vR>0) while energy increases for a contraction (e.g. vL= 0 and vR < 0).

In the adiabatic regime we expect that transitions between modes are severly suppressed and, see problem [*],

(dE/dt) = - 2 (1/L)(dL/dt) E;     (dP/dt) = -  (1/L)(dL/dt) P .

Let us explore this regime and check the range of validity of these formulas. Let us suppose to start with the fundamental state of the system. In adiabatic regime the energy is given by

adiabticEnergy[t_]:= ((π^2)/(2 ((L1[t])^2)))

Let us evolve the system with vR = 0.1 and vR = 1.0

vL1 =0; vR1=0.1; 
Tmax1=0.9functTmax[vL1,vR1];
L1[t_]:= 1 + (vR1 - vL1) t;
sol1= mySolve1[vL1,vR1,Tmax1,(Sqrt[2])Sin[π #1]&];

{interpEnevR01,interpPvR01} = getEnergyMomentum[sol1];

plot1=Plot[{interpEnevR01[t],adiabticEnergy[t]},{t,0,0.2},PlotStyle-> {Black,{Red,Dashed}},AxesLabel-> {"t","E[t]"}]

vR1=1.0;
sol1= mySolve1[vL1,vR1,Tmax1,(Sqrt[2])Sin[π #1]&];

{interpEnevR1,interpPvR1} = getEnergyMomentum[sol1];

plot2=Plot[{interpEnevR1[t],adiabticEnergy[t]},{t,0,0.2},PlotStyle-> {Black,{Red,Dashed}},AxesLabel-> {"t","E[t]"}]

We see that in both cases the adiabatic regime breaks down, i.e.  transistions to higher modes start to be relevant.

Expansion in adiabatic modes

We start now the study of the problem through expansion in adiabatic modes, as described in eq's ()-(). As first task we write a function which gives basis functions and matrices F. For later use  we compute only the time independent part of these matrices, and call them matF1, matF2:

matF1=  ((2k s)/((k^2)-(s^2))) (((-1))^(k-s))  ;  matF2=((2k s)/((k^2)-(s^2)))(-1+ (((-1))^(k-s))))(k≠ s) .

The basis functions are denoted by a vector of functions uB[x,L,xL], the k-th component,  uB[x,L,xL][[k]] is the k-th basis function. L and xL denote the length of the box and the left boundary. We call the function getBasis and as input it needs the number of modes we want to use. As output we get the basis, the matrices and a list of energies for a unit length box.

getBasis[nModes_]:= Module[{uB,matF1,matF2,ene1},
ene1 = (1/2)(π^2)Table[j^2,{j,1,nModes}];
matF1= Normal[SparseArray[{{i_,j_}/;j≠ i->  (((((-1))^(i-j))2 i*j)/((i^2)-(j^2)))},{nModes,nModes}]];
matF2= Normal[SparseArray[{{i_,j_}/;j≠ i->  (((-1+(((-1))^(i-j)))2 i*j)/((i^2)-(j^2)))},{nModes,nModes}]];
uB[x_,ℓ_,xL_]:= Table[(Sqrt[(2./ℓ)])Sin[((π  k (x-xL))/ℓ)],{k,1,nModes}];
{uB,matF1,matF2,ene1}];

Now we have to solve the system of ordinary differential equations (). At first we consider only the case of constant velocities, as we have done in previous section, and write a function mySolModes to perform this task. As input we use: nModes, vL, vR, Tmax, ψIn, respectively the number of modes, velocities of left and right boundaries, maximum time requested for the solution and finally the initial condition, i.e. the wave function at t=0. Let us describe briefly how the function works.

A typical use of this function will be:

{uB, ene1, a0, ak} = mySolModes[numModes, vL, vR, Tmax, ψIn]

mySolModes[nModes_,vL_,vR_,Tmax_,ψIn_]:=Module[
{uB,matF1,matF2,ene1,v,L,a0,functMatF,X,temp},
{uB,matF1,matF2,ene1} = getBasis[nModes];
v = vR-vL;
L[t_]:= 1 + v*t;
a0=Table[NIntegrate[ψIn[x]*uB[x,1,-1/2][[s]],{x,-1/2,1/2},
PrecisionGoal-> 3,AccuracyGoal-> 7],{s,1,nModes}];
functMatF[a_?VectorQ,t_]:= ((v/(L[t]))*matF1 + (vL/(L[t]))matF2).a +I  a*(ene1/((L[t])^2));
{uB,ene1,a0,X/.First[NDSolve[{X'[t] + functMatF[X[t],t]== 0,X[0]==a0},X,{t,0,Tmax},
AccuracyGoal->4,MaxSteps-> 20000]]}]

Some examples

As first case we reproduce the graphs obtained with NDSolve. We use 15 modes, fixed by the parameters numModes below:

vRun =Select[ Range[-5,4,1],#≠0&];
timeRun = 0.5/Abs[vRun];
ψIn[x_]:= (Sqrt[2.])Sin[π  (x+1/2)];
numModes=15;

The following piece of code for the graphs is almost identical to the one used with NDSolve. The only difference in graphics part is that here we have fixed our attention on "laboratory" coordinates x, so we put the initial  box in [-1/2,1/2].

The essential difference is in the computation is how is written the function. Here is obtained by the scalar product of coefficients (Subscript[a, k]) with basis functions. The general instruction is

ak[t].uB[x, L[t], xL[t]]

Here we use Tmax for the time as we want to plot the final shape of |ψ|. The result looks identical to the previous one.

GraphicsGrid[Partition[Table[
v = vRun[[k]];
L[t_]:= 1 + v*t;
xL[t_]:= -1/2;
xR[t_]:=1/2+v*t;
Tmax = timeRun[[k]];
Lb = L[Tmax];
deltaL = 0.02Max[1,Lb];
string1 = "v="<>ToString[vRun[[k]]];
thisLabel = string1;
thisArrow = If[v > 0,{
RGBColor[0,0,1],Thickness[.015],Arrowheads[0.07],Arrow[{{1/2+0.05,1},{1/2+0.3,1}}]},
{RGBColor[0,0,1],Thickness[.012],Arrowheads[0.07],Arrow[{{1/2-0.02,1},{1/2-0.15,1}}]}];
{uB,ene1,a0,ak}=mySolModes[numModes,0,v,Tmax,ψIn];
Show[
Plot[Abs[ψIn[x]],{x,-1/2,1/2},
PlotRange-> All,
Axes-> False,PlotStyle-> {RGBColor[1,0,0],Dashing[{0.02,0.02}]},
PlotLabel-> thisLabel,
Epilog->  {Line[{{-deltaL-1/2,0},{3/2,0}}],
{GrayLevel[0.8],Rectangle[{1/2,0},{1/2+deltaL,1.8}]},
RGBColor[0,0,0],Rectangle[{-1/2-deltaL,0},{-1/2,1.8}],
RGBColor[1,0,0],Rectangle[{xR[Tmax],0},{xR[Tmax]+deltaL,1.8}],
thisArrow}
],
Plot[Abs[ak[Tmax].uB[x,L[Tmax],xL[Tmax]]],{x,xL[Tmax],xR[Tmax]}]],
{k,1,Length[vRun]}],3],ImageSize-> 72*8]//Timing

As done with NDSolve we can collect various animations using a Manipulate instruction. To have a Tmax for the plot which avoid crash of the walls or to much distance covered by the mouving walls we use the function functTmax which choose the time in such a way that everything is limited inthe graphs and walls do not crash.

functTmax[vL_,vR_]:= Module[{v,tL,tR,tCrash},
v= vR-vL;
tL =If[vL≥  0,10^5,-1/vL];
tR =If[vR≤ 0,10^5,1/vR];
tCrash =If[v≥ 0,10^5,-0.5/v];
If[ vL== vR == 0,1, N[Min[{tL,tR,tCrash}]]]]

We choose 30 frames for the animation

nFramesNDSolve = 30;

Manipulate[Module[{L,ψIn,sol,t,dt,temp,Tmax,frame,deltaL,evolutionPlot},
xL[t_]:= -1/2+vL*t;
xR[t_]:=1/2+vR*t;
L[t_]:= 1 + (vR-vL)*t;
Tmax = 0.9functTmax[vL,vR];
Lb = L[Tmax];
dt = Tmax/nFramesNDSolve;
ψIn[x_]:= (Sqrt[2.])Sin[π  (x+1/2)];
{uB,ene1,a0,ak}=mySolModes[numModes,0,v,Tmax,ψIn];
deltaL = 0.02;
evolutionPlot=
Table[t = frame*dt;
Plot[Abs[ak[t].uB[x,L[t],xL[t]]],{x,xL[t],xR[t]},PlotRange-> {{-1.6,1.6},{0,2.5}},
Epilog-> {RGBColor[0,0,0],Rectangle[{-deltaL+xL[t],0},{xL[t],2}],
RGBColor[1,0,0],Rectangle[{xR[t],0},{xR[t]+deltaL,2}]}],{frame,0,nFramesNDSolve}];
ListAnimate[evolutionPlot]
],
{{vL,0},-5,5,0.5,Appearance->"Labeled"},
{{vR,-1},-5,5,0.5,Appearance->"Labeled"},
TrackedSymbols:> {vL,vR}]

Galilei invariance

We owe to the reader the solution of the puzzle about oscillations of wave function for uniform motion of both walls. The solution is that initial condition was wrong, for testing Galilei invariance. In a uniform motion with velocity V also the intial wave function has to move with that velocity. Then the stationary ψIn used above must be mutiplied by the phase factor Exp[i v x] (in our units m = ℏ = 1). Animation belo shows that in this case everything works, and at the same time this is a check for the program. The instructions below are just a slight modification of previous ones, simply vR is forced to be equal to vL.

Manipulate[Module[{L,ψIn,sol,t,dt,temp,Tmax,frame,deltaL,evolutionPlot,v},
v=vL;
xL[t_]:= -1/2+v*t;
xR[t_]:=1/2+v*t;
L[t_]:= 1 ;
Tmax = 0.9functTmax[v,v];
Lb = L[Tmax];
dt = Tmax/nFramesNDSolve;
ψIn[x_]:= (Sqrt[2.])Sin[π  (x+1/2)]Exp[I v*x];
{uB,ene1,a0,ak}=mySolModes[numModes,v,v,Tmax,ψIn];
deltaL = 0.02;
evolutionPlot=
Table[t = frame*dt;
Plot[Abs[ak[t].uB[x,1,xL[t]]],{x,xL[t],xR[t]},PlotRange-> {{-1.6,1.6},{0,2.5}},
Epilog-> {RGBColor[0,0,0],Rectangle[{-deltaL+xL[t],0},{xL[t],2}],
RGBColor[1,0,0],Rectangle[{xR[t],0},{xR[t]+deltaL,2}]}],{frame,0,nFramesNDSolve}];
ListAnimate[evolutionPlot]
],
{{vL,2},-2,2,0.1,Appearance->"Labeled"},
TrackedSymbols:> {vL}]

Energy and momentum

Let us now check energy and momentum variations, eq's(), using expansion in modes. In notebook [*] we have also written an explicit expression for time derivatives of energy and momentum in terms of modes:

(dE/dt)= - (((π^2)(ℏ^2))/((L^2)m)){((Subscript[v, R])/L)((|(∑_k)(((-1))^k) k (Subscript[a, k])|)^2)  - ((Subscript[v, L])/L)((|(∑_k)k (Subscript[a, k])|)^2)}=(( (ℏ^2))/(2m)){(Subscript[v, L])((|ψ'[(Subscript[x, L])]|)^2) - (Subscript[v, R])((|ψ'[(Subscript[x, R])]|)^2)  }

(dP/dt) =   (((ℏ^2)(π^2))/(m (L^3))) (((| (∑_k) (Subscript[a, k])k |)^2)-((| (∑_k) (Subscript[a, k])k (((-1))^k)|)^2)  ) = (( (ℏ^2))/(2m)){((|ψ'[(Subscript[x, L])]|)^2) -((|ψ'[(Subscript[x, R])]|)^2)  }

We take the same initial conditions taken with NDSolve

numModes = 15;
vL1 = -0.1; vR1=1.; 
Tmax1=  0.9functTmax[vL1,vR1];
L1[t_]:= 1 + (vR1 - vL1) t;
xL1[t_]:= -1/2 + vL1*t;
xR1[t_]:= 1/2 + vR1*t;
ψIn[x_]:= (Sqrt[2.])Sin[π  (x+1/2)];

Compute the solution

{uB,ene1,a0,ak}=mySolModes[numModes,vL1,vR1,Tmax1,ψIn];

This is the ψ function

ψ[x_,t_]:= ak[t].uB[x,L1[t],xL1[t]]

As with NDSolve we compute energy and momentum with an interpolation

nptTime = 150; 
listT = Range[0,Tmax1,Tmax1/nptTime];
listEne=Table[ 1/2*NIntegrate[((Abs[(∂_x) ψ[x,t]])^2),{x,xL1[t],L1[t]+xL1[t]}],{t,listT}];
interpEne = Interpolation[Transpose[{listT,listEne}]];
listP = Table[ NIntegrate[Im[Conjugate[ψ[x,t]](∂_x) ψ[x,t]],{x,xL1[t],L1[t]+xL1[t]},
PrecisionGoal-> 4,AccuracyGoal-> 7],{t,listT}];
interpP =  Interpolation[Transpose[{listT,listP}]];

and compare with theoretical expectations

dEdt[t_]:= (1/2)(vL1*(((Abs[(∂_x) ψ[x,t]])^2)/.x-> xL1[t])-vR1*(((Abs[(∂_x) ψ[x,t]])^2)/.x-> xR1[t]));
dPdt[t_]:= (1/2)((((Abs[(∂_x) ψ[x,t]])^2)/.x-> xL1[t])-(((Abs[(∂_x) ψ[x,t]])^2)/.x-> xR1[t]))

The following graphs show the agreement. Below is reported the previous graph obtained wih NDSolve, they are compatible. The reader can check that using NDSolve with PrecisionGoal → 3 the graph is smoother but do not agree so well with mode expansion, i.e. one needs a relatively high precision in the general pourpose routine to reconstruct the result of a few modes (15 in this example) which are particularly well adapted to this problem.

GraphicsGrid[{{
Plot[{interpEne'[t],dEdt[t]},{t,0,Tmax1},AxesLabel-> {"t","dE/dt"},
PlotStyle-> {Black,{Red,Dashed}}],
Plot[{interpP'[t],dPdt[t]},{t,0,Tmax1},AxesLabel-> {"t","dP/dt"},
PlotStyle-> {Black,{Red,Dashed}}]}},
PlotLabel-> "{vL,vR}= "<>ToString[{vL1,vR1}],ImageSize-> 72 8]

The check of the first expression for dE/dt and dP/dt is immediate (and easier to compute):

kModes = Range[1,numModes];
kModesSign = Table[k (((-1))^k),{k,1,numModes}];

dEdt[t_]:=- ((π^2)/((L1[t])^2))( (vR1/(L1[t]))((Abs[kModesSign.ak[t]])^2) - (vL1/(L1[t]))((Abs[kModes.ak[t]])^2));
dPdt[t_]:= - ((π^2)/((L1[t])^3))( ((Abs[kModesSign.ak[t]])^2) - ((Abs[kModes.ak[t]])^2));

GraphicsGrid[{{
Plot[{interpEne'[t],dEdt[t]},{t,0,Tmax1},AxesLabel-> {"t","dE/dt"},
PlotStyle-> {Black,{Red,Dashed}}],
Plot[{interpP'[t],dPdt[t]},{t,0,Tmax1},AxesLabel-> {"t","dP/dt"},
PlotStyle-> {Black,{Red,Dashed}}]}},
PlotLabel-> "{vL,vR}= "<>ToString[{vL1,vR1}],ImageSize-> 72 8]

Macroscopic motion induced by the particle

As last example of motion let us take  a bit unusual case. In one dimension conservations law completely determine the dynamics. If we imagine that the walls have mass M then imposing anergy conservation and momentum conservation we have two equations for vR and vL and the wall motion is self-consistent with the particle motion. The velocities, and by consequence the positions of walls enter in a very nonlinear way as they are part of the definition of basis function. This approach has some vague similarity to Bohm hidden variables model.

A bit simplified version of the argument consists in identifyng left and right wall, i.e. consider the particle on a circle. The radius of the circle  is L and do not varies, as vR = vL.  We have only one function to determine, the overall velocity V, and we need only one equation, which can be taken that coming from total momentum conservation, i.e.

M (dV/dt)+ (dP/dt) = 0.

P is taken as the mean momentum of the particle inside the box.

We choose the number of modes to be used and get corresponding matrices, basis, etc. We put (arbitrarly) Mwall = 4.

To solve equation () and at the same time compute (Subscript[a, k])[t]we have to modify a bit the function mySolModes. We adopt this strategy:

We pass to the new function, which will be named mySolveRing, the initial mode composition.

In notebook [*] it has been shown that matrix elements of momentum are

⟨s|P|k⟩= (i/L)((2 ks)/((s^2)-(k^2)))((((-1))^(k-s))-1)≡(i/L) (Subscript[matF2, sk])

Then we can compute with

(∑_sk)(Subscript[a, s]^( * )) (Subscript[P, sk])(Subscript[a, k])

the initial momentum of the particle and the momentum at later time, using (Subscript[a, k])[t].

From the computed particle momentum we compute wall velocity V = - P/M, if initial momentum is zero, otherwise V = (- P+(Subscript[P, tot]))/M. This fixes initial velocity of the wall, given the total momentum.

We compute V[t] using equation () with dP/dt taken from ().

We assume zero initial total momentum. The length of the box, as explained, is fixed, let us say L = 1.

These are our general parameters

numModes = 20;
Mwall = 10;

To compute the equation of motion for V we need the sums apperaing in () so we define the two vectors:

kModes = Range[1,numModes];
kModesSign = Table[k (((-1))^k),{k,1,numModes}];

The solution for (Subscript[a, k]) coefficients now needs only one term in the F matrix, as dL/dt = 0. The second term, matF2, has to be multiplied by vL, which now is one of the unknown terms in NDSolve.

Up to version 5.2 NDSolve admitted mixed data structures as unknown functions, in our case the vector of coefficients (Subscript[a, k])[t] and the velocity of the wall, v[t]. The first release of Mathematica 6 seems to have problems in this respect, than we simply collect both (Subscript[a, k])and v in a single big vector, the first numModes components contain the coefficients (Subscript[a, k]) and  the last component is v.

A typical call will be

Y = mySolRing[nModes, a0, Tmax];

The derivative of v is computed with :

(dv/dt) + (1/M) (dP/dt) = 0; (dP/dt) =   (((ℏ^2)(π^2))/(m (L^3))) (((| (∑_k) (Subscript[a, k])k |)^2)-((| (∑_k) (Subscript[a, k])k (((-1))^k)|)^2)  ) ;

mySolRing[nModes_,stateVector_,Tmax_]:= Module[{v,Lfix,xTemp,vTemp,varTemp},
Lfix=1.;
a1 =N[ stateVector/Norm[stateVector]];
{uB,matF1,matF2,ene1} = N[getBasis[nModes]];
p0 =Re[ (I/( Lfix))Conjugate[a1].matF2.a1];
v0 = - p0/Mwall;
(* on the first nModes components we have a matrix multiplication. the last term
is used in dv/dt + p/M ==0 , it computes p/M *)
funcMatF[a_?VectorQ]:= (xTemp = a[[1;;nModes]] ;vTemp= a[[-1]];
varTemp=((Abs[Total[kModes*xTemp]])^2)-((Abs[Total[kModesSign*xTemp]])^2);
Append[vTemp(matF2.xTemp)+I xTemp*ene1,+ π^2/Mwall*varTemp]);
Y/.First[NDSolve[{Y'[t]+funcMatF[Y[t]]==0,Y[0]== Append[a1,-p0/Mwall]},Y,{t,0,Tmax},
AccuracyGoal->4,MaxSteps-> 20000]]
];

We choose an initial state  (it will be normalized in mySolRing ) and the total evolution time

a0 = Table[0,{numModes}]; a0[[1]] = 1; a0[[2]]=0.7+I/2;
TmaxWithWalls = 3;

and compute the evolution

ysol = mySolRing[numModes,a0,TmaxWithWalls];

Let us extract for clarity the coefficients end the velocity :

bk[t_?NumberQ]:= ysol[t][[1;;numModes]];
vWall[t_?NumberQ]:= ysol[t][[-1]];

We changed the usual name for coefficients and named them (Subscript[b, k]). The word (Subscript[a, k])was used before and as been "Tagged" as an Interpolating Function, we have or remove it or change convention, we have choosen the second possibility. We have noticed this because if the reader make changes in the Notebook he must pay attention at naming variables.

Here is a graph showing the variations in velocity induced on the wall :

Plot[vWall[t],{t,0,TmaxWithWalls},AxesLabel-> {"t","vWall"}]

We can compute separate energies, etc.

EneWall[t_]:= (1/2)Mwall ((Re[vWall[t]])^2);
EneParticle[t_]:= Re[Conjugate[bk[t]].(ene1*bk[t])];

GraphicsGrid[{{Plot[EneParticle[t],{t,0,TmaxWithWalls},AxesLabel-> {"t","Epart"}],
Plot[EneWall[t],{t,0,TmaxWithWalls},AxesLabel-> {"t","Ewall"}]}},ImageSize-> 72 8]

The reader can check that the sum of energies is numerically constant.

A problem left o the reader is the sudy of time averaged quantities.

Oscilations and transitions

We can allow more complicated motions to box boundaries, in particular periodic motions can be seen as a mechanical prototype for resonant induced transitions between modes of the box. We consider small oscillations of the total length L[t], around a mean value 1:

L[t] = 1 + A  Sin[Ω t].

We choose two different ways for obatining ()

In the first case only parity preserving transitions are induced, in the second all transitions are allowed.

We write a function  myOscillatingBox which compute time evolution for modes admitting a Boolean variable, leftOsc, distingushing the two cases. leftOsc = False means fixed left boundary. The structure of the function must be obvious at this point of the problem. It needs as input the initial vector of coefficients (Subscript[a, k])the frequency, the amplitude, the time of evolution and the Boolean variable. We use the fact that function getBasis compute matrices (Subscript[F, ks]) without  time variable prefactors, v/L etc. and we can make a call to this function outside myOscillatingBox. These matrices, and the number of modes, will be global variables are defined at first. It would be better to pass these values to mySolOscillating but here we choose to work with global defined variables to simplify the life to the reader. We choose to work with only the lowest 8 modes.

numModes = 8;
{uB,matF1,matF2,ene1} = getBasis[numModes];
Tmax = 9;

myOscillatingBox[stateVector_,Ω_,amp_,leftOsc_]:= Module[
{X,vL,vR,v,a1},
a1=  stateVector/Norm[stateVector];
{uB,matF1,matF2,ene1} = getBasis[numModes];
L[t_]:= (1+amp*Sin[Ω t]);
v[t_]:=  amp*Ω*Cos[Ω t];
If[leftOsc == False,
vL[t_]:= 0;vR[t_]:= amp*Ω*Cos[Ω t];,
vR[t_]:= 0.5amp*Ω*Cos[Ω t];vL[t_]:= -0.5amp*Ω*Cos[Ω t];];
functMatF1[a_?VectorQ,t_]:= ((v[t])/(L[t]))matF1.a  + ((vL[t])/(L[t]))matF2.a;
functMatF2[a_?VectorQ,t_]:=I  a*(ene1/((L[t])^2));
X/.First[NDSolve[{X'[t] +functMatF1[X[t],t]+ functMatF2[X[t],t]== 0,
X[0]==a1},X,{t,0,Tmax},MaxSteps-> 20000
]]];

Consider to start with an oscillation in resonance with the transition 1→2 and start with a system in state 1 (the ground state).

amplitude = 0.05;
a0 = Table[0,{numModes}];
a0[[1]]=1;
Ω = (ene1[[2]]-ene1[[1]]);

aModes= myOscillatingBox[a0,Ω,amplitude,False];

To help the reader in distinguishing between modes we use different colors and add a small legend. For | (Subscript[a, k]) | higher than 0.05 we add the number of the mode to the curve. All this stuff is done in a function showModes[k] which plot modes up to k.

modeColor=Table[RGBColor[0,0,0],{10}];
modeColor[[Range[2,5]]] = {RGBColor[1,0,0],RGBColor[0,1,0],RGBColor[0,0,1],RGBColor[1,0,1]};
modeLegend[k_]:={Table[{Graphics[{modeColor[[j1]],Line[{{0,j1},{j1,j1}}]}],ToString[j1]},{j1,1,k}],
LegendPosition-> {1.05,-0.3},LegendSize-> {0.25,0.6}};

showModes[k_]:= ShowLegend[
data = ((Abs[aModes[Tmax]])^2);
ss=Table[If[data[[jj]] > 0.05,jj,""],{jj,1,numModes}];
Show[Table[Plot[((Abs[aModes[t][[j]]])^2),{t,0,Tmax},PlotStyle->modeColor[[j]]],{j,1,k}],
Graphics[Table[Text[ToString[ss[[k1]]],{1.01Tmax,1.05data[[k1]]}],{k1,1,k}]],
PlotRange-> {{0,Tmax},{0,1.05}},AxesOrigin-> {0,0},
AxesLabel-> {"t","|(Subscript[a, k])((|)^2)"},ImageSize -> 72 7],modeLegend[k]]

showModes[5]

We see a classical flip-flop between the states 1 and 2. The survaving probability for initial state must change as Cos[ℱ ((t])^2) wher ℱ is the "strength" of the oscillatory perturbation term. In Problems [*] it has been shown that for a transition k→s  we have

ℱ = ((k s)/((k^2)-(s^2)))((Ω A )/L)

where A is the amplitude. Te probabilities must vanish for T = π/(2ℱ), in our case

Tminimum = (π/2)/(((k s)/((k^2)-(s^2)))((Ω A )/L))/.{L-> 1,A-> amplitude,k-> 2,s-> 1}

From our calculation :

prob[k_?IntegerQ,t_?NumberQ]:= Abs[aModes[t][[k]]]^2;

FindMinimum[prob[1,t],{t,Tminimum,2,4},AccuracyGoal-> 4]

which is in good agreement.

In the transition 1-> 2 is parity has changed, so it must disappears with a symmetric oscillations of the walls. To produce a symmetric oscillations we set to True the last parameter in mySolOscillating. This is the result:

aModes= myOscillatingBox[a0,Ω,amplitude,True];
showModes[5]

As a second example we consider the parity preserving transition 1-> 3

Ω = (ene1[[3]]-ene1[[1]]);
aModes= myOscillatingBox[a0,Ω,amplitude,True];
showModes[5]

Apart from expected resonance 1-3 we see a substantial contribution from state n=5. In units (π^2)/2(L^2) the energy differences for the former and the second transition are ((3^2) - 1) = 8 and ((5^2)- 1) = 24, then a third armonic effect is in resonance for 1-5 transition. That the effect is a non linear one is apparent if we half the amplitude:

amplitude = 0.025;
Ω = (ene1[[3]]-ene1[[1]]);
aModes= myOscillatingBox[a0,Ω,amplitude,True];
showModes[5]

The same non linear resonance can be directly tested using a "non resonant" frequency, as ((Subscript[E, 3])- (Subscript[E, 1]))/2 here the transition must appear as a double harmonic:

amplitude = 0.08;
Ω = (ene1[[3]]-ene1[[1]])/2;
aModes= myOscillatingBox[a0,Ω,amplitude,False];
showModes[5]

Exercises