Particle on a ring

Quantum Mechanics
K. Konishi, G. Paffuti

Study the motion of a particle along a ring with non periodic boundary conditions varying with time.

Particle on a ring

Packages

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

Short review of the theory

Block[{pt,delta=1.2`,R=5},pt={-R,0,0};Show[Graphics3D[{EdgeForm[],RGBColor[0,0,0],Specularity[RGBColor[0,0,0],2.8`],{RGBColor[0,0,1],Polygon[{pt+{-delta,0,-delta},pt+{delta,0,-delta},pt+{delta,0,delta},pt+{-delta,0,delta}}]},RGBColor[1,0,0],Thickness[0.01`],Line[Table[{R Cos[((2 π k)/100)],R Sin[((2 π k)/100)],0},{k,0,100}]]},SphericalRegion->False,ViewPoint->2 {2,3,3},Boxed->False]]]

Let φ the angle along the circle. The system can be described by:

H = ((ℏ^2)/(2m (R^2)))((((1/i)(∂/(∂φ))))^2);            ψ[2π] = (e^(-i 2π α ))ψ[0].

In the following we use units ℏ = 1, m=1, R=1, φ → x. Eigenstates and eigenvalues of H  are given by

(Subscript[u, n])[x;α] = (1/(Sqrt[2π ]))Exp[i n x - i α x];        (Subscript[E, n])= (1/2)(((n-α))^2);

In the following we will study numerically the time evolution of a state |n> when α varies with time. For simplicity we limit ourselves to a linear time dependence, starting from α = 0,

α = 0 + v*t;               0 ≤ t ≤ Tmax

The general form of a state at time t is

Ψ[x,t] = (∑_n)(Subscript[b, n])[t](Subscript[u, n])[x;α[t]]

where (see Problems)

((Subscript[db, k])/dt) + (∑_s)(Subscript[F, ks])(Subscript[b, s]) = - i (Subscript[E, k])[α](Subscript[b, k]);
(Subscript[F, ks])=(1/(k-s))((dα)/dt) (1-(Subscript[δ, ks])) + (Subscript[δ, ks])(-i π ((dα)/dt)).

In the adiabatic approximation the time evolution of a state |n> involves only the diagonal term of the above matrix (the state remain frozen with the Hamiltonian) and the solution of () reduces to

(Subscript[b, k])[t] = Exp[-i (∫_0^t)(Subscript[E, k])[α]𝕕t] Exp[i π α[t]] ≡ Exp[-i ((Subscript[δ, k]) - γ)].

The first factor is the usual dynamical phase, the second one is a Berry phase.

In the following we will provide two different methods of numerical solution. The first one use directly NDSolve, the second one use the mode expansion.

NDSolve

The boundary condition () should be written as ψ[2π,t] = Exp[-i α 2π] ψ[0,t], for all times, and in this form is not accepted by NDSolve. We can avoid the problem by a change of variable

ψ[x,t] = Exp[-i x α] f[x,t],

f[x,t] is now periodic. It satisfy a Schrödinger equation obtained by substitution in the Schrödinger equation for ψ:

i (∂_t) f + ((dα)/dt) x f = - (1/2)((∂_(x,x))f -2 i α f - (α^2) f)

We will start from a given state, nStart, and follow the evolution with α = v t, up to Tmax =  1/|v|. The final α is changed by 1. The initial condition for f[x,t] will be of the form:

fIn[x,t] = (Subscript[u, nStart])[x,t]*Exp[i x α[0]] = (Subscript[u, nStart])[x,t].

We take dα/dt = v  =< 0 , in such a way that for adiabatic evolution n → n+1 is expected.

For linear change in time the total phase acquired by the state n in adiabatic approximation is

(Subscript[δph, n])[t] = (1/2)(∫_0^t) (((n- v t))^2) - π v t = (1/2) ((n^2) t-n (t^2) v+(((t^3) (v^2))/3)) - π v t =
                  - (1/(6v))(((( n - t v))^3) - (n^3)) - π v t .

Parameters

nStart = 1;
v = {-5.,-0.5,-0.1,-0.05,-0.01,-0.005};
Tmax = Abs[1/v];

δph[t_]:= (1/2) ((nStart^2) t-nStart (t^2) v+(((t^3) (v^2))/3)) - π v t;

Initial function and adiabatic state

fIn[x_]:= Exp[I nStart x];

ψAd[x_,t_,v_]:=  Exp[I (nStart - v*t) x]

Evolution

mySol computes a vector of solutions, one for each velocity. MaxSteps → 30000 has been selected to handle the low velocity case. The adiabatic limit (v -> 0) requires more computational effort.

mySol=Table[
αdot = v[[k]];
eq =  {I (∂_t) f[x,t] +αdot x*f[x,t]+ 1/2((∂_(x,x))f[x,t]-
2I αdot*t*(∂_x) f[x,t] -((αdot)^2) (t^2) f[x,t])== 0,
 f[x,0]== fIn[x],
f[2π,t]== f[0,t]};
 (f/.First[NDSolve[eq,f,{x,0,2π},{t,0,Tmax[[k]]},
MaxSteps-> 30000,
Method->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", 
"MinPoints"-> 30,DifferenceOrder->"Pseudospectral" 
}}]]),
{k,1,Length[v]}];

Graphics

We report the real part of the function after a  change of 1 in α. The function needed to come back to original wave function is

g[x_,t_,αdot_]:= Exp[- I αdot *t *x]

In the graph we show:

A quantitative check for adiabatic behavior can be the value of the function at the origin. In the normalization used (without 1/(Sqrt[2π]) )

(Subscript[ψ, ad])[x,t]= (e^(i(n-vt)x)) (e^(-i (Subscript[δ, T])[t]));  ℝe[(Subscript[ψ, ad])[0,t]] = Cos[(Subscript[δ, T])[t]];

The phase δ is the sum of dynamic phase and Berry phase, as defined above.

Show[GraphicsGrid[Partition[Table[
string1 = "v="<>ToString[v[[k]]];
string2 = " ψ[0]="<>ToString[
NumberForm[Re[mySol[[k]][0,Tmax[[k]]]*g[0,Tmax[[k]],v[[k]]]],3]];
string3 = " ad: "<>ToString[NumberForm[Cos[δph[Tmax][[k]]],3]];
thisLabel = string1<>string2<>string3;
Plot[{
Re[fIn[y]],
Re[ψAd[y,Tmax[[k]],v[[k]]]*Exp[-I δph[Tmax[[k]]][[k]]]],
Re[mySol[[k]][y,Tmax[[k]]]*g[y,Tmax[[k]],v[[k]]]]},{y,0,2 π},
PlotLabel-> thisLabel,DisplayFunction-> Identity,
PlotStyle-> {{RGBColor[0,0,1],Dashing[{0.05,0.05}]},
{RGBColor[1,0,0],Dashing[{0.02,0.02}]},RGBColor[0,0,0]}],
{k,1,Length[v]}],2]],ImageSize-> 72*8]

The qualitative picture is quite clear: for high velocity we are in an almost sudden change of the Hamiltonian and the wave function tend to be frozen, for low velocities we approach the adiabatic evolution.

This animation is to give a qualitative idea of the change of the wave function. Let us note that for small dα/dt, the effective wave number has a small change and the overall change of ψ is relatively small. For low v the wave number double (in our case) which means in the temporal domain a much faster wave. We give a parameter for the plot, frozenWave, if True the adiabatic phase is subtracted and ψ more or less is a steady wave, if frozenWave = False, we plot the whole time development, i.e. we have traveling waves and the reader can verify that they seem to go much faster on the screen.

selectedWave = 4;
numFrames = 100;
frozenWave = True;

evolutionPlot = Block[{j,dt,par},
(If[frozenWave,par =1;,par=0;];
j = Min[selectedWave,Length[v]];
dt = Tmax[[selectedWave]]/numFrames;
Table[Plot[
Re[mySol[[j]][y,dt*k]*
g[y,dt*k,v[[j]]]
Exp[I* par* δph[dt*k][[j]]]],{y,0,2 π},
PlotRange-> {{0,2π},{-1.5,1.5}}],{k,0,numFrames}])];

ListAnimate[evolutionPlot,AnimationRunning->False]

Mode expansion

We expand the wave function in a series of adiabatic eigenstates () and solve the system of ordinary differential equations (). The number of modes is nModes, ranging from -nP to +nP.  To enumerate these modes we will use positive integers, as usual in Mathematica , the state corresponding to index value k is  k - nP -1, i.e. the first is -nP and the last is +nP. startMode, which will put 1 in the following, will denote our starting state, the corresponding index is startMode + nP+1.

startMode = 1;
nP = 5;nModes = 2 nP + 1;

The assign the initial value for the coefficients  and adiabatic phases that will be used as a check. v is the rate of change for α.

b0 = Table[0,{nModes}]; b0[[nP+1+startMode]]=1.;
δ[t_,v_]:=Table[ (1/(6v))( (((k-nP-1))^3) - ((((k-nP-1)-v t))^3)) - π v*t,{k,1,nModes}];

The adiabatic base is

ψAd[x_,t_,v_]:= Table[Exp[I ((k-nP-1) - v*t)x],{k,1,nModes}];

Solution

v = {-5.,-0.5,-0.1,-0.05,-0.01,-0.005};

ene=(1/2)Table[(((k-nP-1 - α))^2),{k,1,nModes}];

The solution of equations  () can be obtained using a vector function X[t] (the coefficients) for NDSolve. The function functMatF which just do matrix multiplication is written in such a way that NDSolve take its argument only if X[t] is a vector, as it must be. The matrix (Subscript[F, ks]) for equation () is computed inside this function.

mySolModes=Table[
ω = v[[s]];
matF= Normal[SparseArray[{{i_,i_} ->   -I π  ω,{i_,j_}/;j≠ i->  (ω/(i-j))},{nModes,nModes}]];

functMatF[c_?VectorQ,t_]:= matF.c +I  c*(ene/.α->ω t);
X/.First[NDSolve[{X'[t] + functMatF[X[t],t]== 0,X[0]==b0},X,{t,0,1/Abs[ω]}]],
{s,1,Length[v]}];

Module[{jS,ω,b,Tplot},
GraphicsGrid[Partition[Table[
jS=startMode+nP+1;
ω = v[[j]];
b = mySolModes[[j]];
Tplot = 1/Abs[ω];
string1="v="<>ToString[ω];
string2=" ψ[0]="<>ToString[NumberForm[
Re[ψAd[0,ω,Tplot].b[Tplot]],3]];
string3=" ad: "<>ToString[NumberForm[Cos[δ[Tplot,ω][[jS]]],3]];
thisLabel=string1<>string2<>string3;
Plot[{
Re[ψAd[y,ω,0][[jS]]],
Re[ψAd[y,ω,Tplot][[jS]]* Exp[-I δ[Tplot,v[[j]]][[jS]]] ],
Re[ψAd[y,ω,Tplot].b[Tplot]]},{y,0,2 π},
PlotLabel-> thisLabel,DisplayFunction-> Identity,
PlotStyle-> {{RGBColor[0,0,1],Dashing[{0.05,0.05}]},
{RGBColor[1,0,0],Dashing[{0.02,0.02}]},RGBColor[0,0,0]}],
{j,1,Length[v]}],2],ImageSize-> 72*8]]

We can follow the history in time of the  occupation probabilities of the modes, |(Subscript[b, k])((|)^2) .

data = Table[Table[Abs[mySolModes[[k]][1/Abs[v[[k]]]][[s]]]^2,{s,1,nModes}],{k,1,Length[v]}];

GraphicsGrid[Partition[Table[
Show[
ss=Table[If[data[[k]][[jj]] > 0.1,jj-nP-1,""],{jj,1,nModes}];
Table[
Plot[Abs[mySolModes[[k]][t][[j]]]^2,{t,0,1/Abs[v[[k]]]},
PlotRange-> {All,{0,1.05}},PlotStyle-> Hue[0.1*j]],{j,1,nModes}],
Graphics[Table[
Text[ss[[k1]],{1./Abs[v[[k]]],data[[k]][[k1]]}],{k1,1,nModes}]]
,
AxesLabel-> {"t","(Subscript[c, k]^2)"},PlotLabel->"v = "<> ToString[v[[k]]],
DisplayFunction-> Identity],{k,1,Length[v]}],2],ImageSize-> 72*8]

Evolution of momentum

From mode expansion it is possible to compute the evolution of every operator. Let us consider in particular the momentum operator P. As states () are eigenstates of P we have:

<Ψ[t]|P|Ψ[t]> = (∑_k)|(Subscript[b, k])[t]((|)^2)(k-α[t])

In the problems it has been shown that averaging over fast oscillations

((d<P>)/dt) = -(( dα)/dt) ⇒ <P[t]> = <P[0]> - (( dα)/dt) t.

In this section these equations will be verified on a simple state of the form |1> + 1.2 |3>.  Let us proceed in exactly the same way as in the previous section, the only difference being the initial state. As usual we do not care about overall normalization.

b0 = Table[0,{nModes}]; 
startMode=1;
b0[[nP+1+startMode]]=1.; b0[[nP+1+startMode + 1]]=1.2;

v = {-5.,-0.5,-0.1,-0.05,-0.01,-0.005};
ene=(1/2)Table[(((k-nP-1 - α))^2),{k,1,nModes}];
pModes =Table[k-nP-1 - α,{k,1,nModes}];

We find coefficients b[t] with the new initial condition.

mySolModes=Table[
ω = v[[s]];
matF= Normal[SparseArray[{{i_,i_} ->   -I π  ω,{i_,j_}/;j≠ i->  (ω/(i-j))},{nModes,nModes}]];

functMatF[c_?VectorQ,t_]:= matF.c +I  c*(ene/.α->ω t);
X/.First[NDSolve[{X'[t] + functMatF[X[t],t]== 0,X[0]==b0},X,{t,0,1/Abs[ω]}]],
{s,1,Length[v]}];

Graphics

We plot <P> vs t for different values of dα/dt. The initial and final momentum are obtained averaging over modes the vector pModes, which depends on α, at initial time α=0 and at final time, in our case, α = 1/Sign[v], i.e. ±1 depending on sign of v.

p0 =Table[ Re[
Conjugate[mySolModes[[k]][0]].(mySolModes[[k]][0]*pModes/.α->0)],{k,1,Length[v]}]

p1 =Table[ T1= 1/Abs[v[[k]]];
 Re[
Conjugate[mySolModes[[k]][T1]].(mySolModes[[k]][T1]*pModes/.α-> Sign[v[[k]]])],{k,1,Length[v]}]

The dashed line in the plot is the adiabatic prediction (or what we expect after average on high frequency modes).

Show[GraphicsGrid[Partition[Table[b=mySolModes[[k]];ω=v[[k]];Plot[Re[Conjugate[b[t]].(b[t] pModes/.α->ω t)],{t,0,(1/(Abs[v[[k]]]))},AxesLabel->{"t","P"},PlotLabel->"v = "<>ToString[ω],DisplayFunction->Identity,Epilog->{RGBColor[1,0,0],Dashing[{0.02,0.02}],Line[{{0,p0[[k]]},{(1/(Abs[ω])),p1[[k]]}}]}],{k,1,Length[v]}],2]],ImageSize->72 6]

As a last example of calculation let us verify explicitly the equations of motion for P. As an operator (see problems)

((d P)/dt) = -((dα)/dt)2πδ(x).

The mean value is simply obtained:

(∑_k) ((Subscript[b, k]^('* ))(Subscript[b, k])+(Subscript[b, k]^( * )) (Subscript[b, k]^(')))(k-α) - ((dα)/dt)(∑_k) |(Subscript[b, k])((|)^2) = - ((dα)/dt)|(∑_k) (Subscript[b, k])((|)^2).

We plot the two sides of the equation:  they are indistinguishable. The graphs also give an idea on the force acting on the particle. The rapid oscillations as dα/dt → 0 suggest a quite singular adiabatic limit.

GraphicsGrid[Partition[Table[b=mySolModes[[k]];ω=v[[k]];Plot[{Re[Conjugate[b'[t]].(b[t] pModes/.α->ω t)]+Re[Conjugate[b[t]].(b'[t] pModes/.α->ω t)]-ω Re[Conjugate[b[t]].b[t]],-ω ((Abs[Total[b[t]]])^2)},{t,0,(1/(Abs[ω]))},PlotStyle->{{},{RGBColor[1,0,0],Dashing[{0.02}]}},AxesLabel->{"t","dP/dt"},PlotLabel->"v = "<>ToString[ω],PlotRange->All,DisplayFunction->Identity],{k,1,Length[v]}],2],ImageSize->72 7]