Quantum Modes

AP1601

Introduction

One of the well-known consequences of quantum mechanics is the "energy level." Within potential wells (or fixed domains), oscillations of the wave-function can be expressed as a collection of discrete "modes" of the system. These are called Eigensystems, and Mathematica has built-in commands to find these normal modes and frequencies.

For our quantum ball confined to a box, Schrodinger's equations can be approximated on a grid. The matrix form of these coupled equations are:
	I (∂_t) ψ =  M . ψ
where we have again assigned ℏ = m = 1 for simplicty.This equation appears remarkably simple, but, remember that Schrodinger's wavefunction, ψ,  is complex. 
The time dependence of an eigenmode of fixed frequency (i.e. energy) varies as
	ψ(x, t) ∝ Exp[- I ω t] g(x),
where g(x) is the time-independent spatial struncture of the normal mode.The quantum fluctuations of the wavefunction are faster at higher energy and slower at lower energy. Substituting this expression into our matrix version of Schrodinger's equation gives
	ω ψ =  M . ψ
The eigenvalues of Schrodinger's equation are proportional to the frequency of the oscillations of the wavefunction and to the energy of the eigenmode.

Our Quantum Model (on a grid)

Last week, we introduced Schrodinger's equation on a grid. Computers can not represent a continuous range of numbers (i.e. x) unless they deal with a symbolic or algebraic solution to the quantum problem. The approach we have taken is to replace the x-axis with a grid with the wavefunction represented by complex numbers at each grid point. Spatial derivatives, (∂_x) (∂_x) ψ, are replaced with differences of ψ evaluated on neighboring grid points. (This creates an equation having the same form as the equation for the acceleration of balls on a string.) The equations for the wavefunction at the grid points are:
	I (∂_t) ψ[1][t]  = - (1/2) (- 2 ψ[1][t] +  ψ[2][t]) + V[2] ψ[2][t]
	I (∂_t) ψ[i][t]  = - (1/2) ( ψ[i-1][t] - 2 ψ[i][t] +  ψ[i+1][t]) +V[i] ψ[i][t]
	I (∂_t) ψ[N][t]  = - (1/2) ( ψ[N-1][t] - 2 ψ[N][t] ) + V[N] ψ[N][t]
Notice the equations at the two ends of the grid, i = 1 and i = N, differ from the equations for quantum oscillations located at points separated from the ends.
Important: Remember, when we solve Schrodinger's equation on a grid, this is only an approximation. In reality, space is continuous. The grid approximation is valid only when the quantum oscillations have a wavelength much longer than the spacing between grids. As the energy and momentum of the quantum ball increases, the wavelength of its wavefunction decreases. Therefore, at some level as the energy increases, our computer solution becomes invalid. The finer grid spacing we can use (limited by the speed of our computer and the amount of available computer memory), the more accurate our computer-generated wavefunctions become!

The matrix form of these equations are:

nGrid = 50;

mMatrix = Table[Which[j == i, 1+ v[i], j == i+1, -0.5, j == i-1, -0.5, True, 0], {i,1,nGrid}, {j, 1, nGrid}];

Shallow[mMatrix,{6,10}]

This is a "matrix representation" of Schrodinger's equation expressed on a grid. If the potential is zero everywhere, this matrix has identical eigenmodes as the matrix used to represent a string. 
As mentioned above, the eigenvalues must be interpreted differently for quantum oscillations. The eigenvalues are proportional to the energy, e = ω; whereas, for the string modes, the eigenvalues where proportional to (ω^2).

Energy Levels

Since we are solving quantum mechanics on a grid, the number of normal modes (and, therefore, the number of energy eigenvalues) are are equal to the number of grid points.

When the potential vanishes everywhere within the quantum box,

?Eigensystem

{values, vectors} = Eigensystem[Evaluate[mMatrix /. v[_] -> 0// N]];

Mathematica orders the modes from the most complex (or rapidly varying) mode to the mode having the slowest spatial variations:

ListPlot[values,PlotStyle->{PointSize[0.02]},AxesLabel->{"Mode","Energy"}]

Show[%,PlotRange->{{30,50},{0,0.7}}]

The lowest energy (frequency) mode corresponds to the longest wavelength. 
Notice also, the lowest energy possible for a ball is not equal to zero!

values[[nGrid-5]]

If there is a ball in the box, it can not have zero energy.

ListPlot[vectors[[nGrid-5]],PlotRange->{{0,nGrid+1},{-0.4,0.4}},PlotStyle->PointSize[0.02],Frame->True,Ticks->None,FrameTicks->None]

The first six modes are

modes=Table[ListPlot[vectors[[i]],PlotRange->{{0,nGrid+1},{-0.4,0.4}},Joined->True,PlotStyle->{Hue[((30-i)/5)],Thickness[0.01]},Frame->True,Ticks->None,FrameTicks->None,DisplayFunction->Identity],{i,nGrid,nGrid-5,-1}];

Show[modes,DisplayFunction->$DisplayFunction]

The probabilities for finding the quantum ball at any given position within the box are:

allModes=Table[ListPlot[((Abs[vectors[[i]]])^2),PlotRange->{{0,nGrid+1},{0,0.05}},Joined->True,PlotStyle->{Hue[((30-i)/5)],Thickness[0.01]},Frame->True,Ticks->None,FrameTicks->None,DisplayFunction->Identity, FrameLabel -> "Mode = "<>ToString[nGrid-i+1]<>" Energy = "<>ToString[values[[i]]]],{i,nGrid,nGrid-30,-1}];

Manipulate[Show[allModes[[n]]], {n,1,30,1}]

Momentum representation

As explained during last week's class...
The momentum may be positive or negative depending upon whether the ball is moving to the right or to the left. The function momentumψ re-arranges the output from Fourier so that negative momentum is followed by the positive momentum values.

rearrangeFourier[p_List] := Reverse[p[[Range[nGrid/2 + 1, nGrid]]] ~ Join~ p[[Range[1,nGrid/2]]]]

momentumψ[ψ_List] := rearrangeFourier[Fourier[ψ]];

The corresponding values of momentum (or k) are given by the list below:

klist = Table[(2π/nGrid)i // N, {i, -nGrid/2+1,nGrid/2}];

plotMomentum[ψ_List,options___]:=Module[{mlist},mlist=momentumψ[ψ];ListPlot[Table[{klist[[i]],((Abs[mlist[[i]]])^2)},{i,1,nGrid}],Joined->True,Ticks->{{-π,0,π},None},options]]

plotMomentum[vectors[[nGrid-2]],DisplayFunction->Identity,Frame -> True,PlotRange -> All,FrameLabel -> "Mode = "<>ToString[nGrid-(nGrid-2)+1]<>" Energy = "<>ToString[values[[nGrid-2]]]]

As the energy increases, so does the momentum...

kModes=Table[plotMomentum[vectors[[i]],DisplayFunction->Identity,PlotRange -> {0,0.4},Frame -> True,FrameLabel -> "Mode = "<>ToString[nGrid-i+1]<>" Energy = "<>ToString[values[[i]]]],{i,nGrid,nGrid-30,-1}];

Manipulate[Show[kModes[[n]]], {n,1,30,1}]

Combination of Two Modes

Recall that in quantum mechanics we are trying to compute the probabilities for values of our measurements of the particle's position, momentum, and energy. If the ball were arranged so that it would be likely to be located within a particular eigenmode, then the ball is considered to be within a particular energy state.
The time dependence of the wavefunction of a single energy state is 
	ψ(x, t) ∝ Exp[- I ω t] g(x),
where g(x) is the time-independent spatial struncture of the normal mode.

Remember the identities

Exp[-I ω t] // ExpToTrig

((Abs[Exp[- I]])^2)

When a ball is likely  to be within a single eigenmode, the probabilities for making any particular measurement of the ball will not change in time.

What about a ball that is equally likely to be within two eigenmodes?
In this case,

twoψ[t_] := (vectors[[nGrid]] Exp[-I values[[nGrid]] t] + vectors[[nGrid-1]] Exp[-I values[[nGrid-1]] t])/Sqrt[2]

ListPlot[((Abs[twoψ[0.0]])^2),Joined->True]

Animation of Probability

ListAnimate[Table[ListPlot[((Abs[twoψ[t]])^2),PlotRange->{{0,nGrid+1},{0,0.1}},Joined->True,Frame->True,Ticks->None,FrameTicks->None],{t,0,1000,25}], AnimationRunning -> False]

Animation of WaveFunction

plotψ[ψ_List,options___]:=Module[{p1,p2,p3,p4},p1=ListPlot[Re[ψ],PlotStyle->Blue,Joined->True,options];p2=ListPlot[Im[ψ],PlotStyle->Red,Joined->True,options];p3=ListPlot[Abs[ψ],Joined->True,options];p4=ListPlot[-Abs[ψ],Joined->True,options];Show[p1,p2,p3,p4,options]]

ListAnimate[Table[plotψ[twoψ[t], PlotRange-> {-0.3, 0.3}],{t,0,1000,25}], AnimationRunning -> False]

A Potential Barrier

For your second Mathematica project, some of you may examine the influence of nonuniform potentials on quantum oscillations. 
Using Eigensystem, the energy levels are relatively easy to find.
First, we define a potential barrier, and then use Eigensystem to determine the energy levels and the spatial structures of the corresponding wavefunctions.

vBarrier[i_Integer] := If[i > nGrid/2 -3 && i < nGrid/2 +3, 0.2, 0.0]

ListPlot[Table[vBarrier[i],{i,1,nGrid}],PlotStyle->PointSize[0.02]]

When the potential vanishes everywhere within the quantum box,

{values, vectors} = Eigensystem[Evaluate[mMatrix /. v[i_] ->vBarrier[i]// N]];

Mathematica orders the modes from the highest frequency (i.e. highest energy) to the lowest:

ListPlot[values,PlotStyle->{PointSize[0.02]},AxesLabel->{"Mode","Energy"}]

Notice the barrier creates energy levels having nearly identical eigenvalues.

ListPlot[values,PlotStyle->{PointSize[0.02]},PlotRange->{{30,50},{0.,0.5}},AxesLabel->{"Mode","Energy"}]

allModes=Table[ListPlot[((Abs[vectors[[i]]])^2),PlotRange->{{0,nGrid+1},{0,0.1}},Joined->True,PlotStyle->{Hue[((30-i)/5)],Thickness[0.01]},Frame->True,Ticks->None,FrameTicks->None,DisplayFunction->Identity, FrameLabel -> "Mode = "<>ToString[nGrid-i+1]<>" Energy = "<>ToString[values[[i]]]],{i,nGrid,nGrid-30,-1}];

Manipulate[Show[allModes[[n]]], {n,1,30,1}]

A Potential Well

Start with a potential square well with a depth of -0.2.

vBarrier[i_Integer] := If[i > nGrid/2 - 5 && i < nGrid/2 +5, -0.2, 0.0]

ListPlot[Table[vBarrier[i],{i,1,nGrid}],PlotStyle->PointSize[0.02]]

When the potential vanishes everywhere within the quantum box,

{values, vectors} = Eigensystem[Evaluate[mMatrix /. v[i_] -> vBarrier[i]// N]];

Mathematica orders the modes from the highest frequency (i.e. highest energy) to the lowest:

ListPlot[values,PlotStyle->{PointSize[0.02]},AxesLabel->{"Mode","Energy"}]

allModes=Table[ListPlot[((Abs[vectors[[i]]])^2),PlotRange->{{0,nGrid+1},{0,0.2}},Joined->True,PlotStyle->{Hue[((30-i)/5)],Thickness[0.01]},Frame->True,Ticks->None,FrameTicks->None,DisplayFunction->Identity, FrameLabel -> "Mode = "<>ToString[nGrid-i+1]<>" Energy = "<>ToString[values[[i]]]],{i,nGrid,nGrid-30,-1}];

Manipulate[Show[allModes[[n]]], {n,1,30,1}]

Notice the potential well creates two "negative energy" levels. As the potential well becomes deeper (i.e.  more negative), additional modes can become trapped within the potential well.

ListPlot[values,PlotStyle->{PointSize[0.02]},PlotRange->{{30,50},{-0.4,0.4}},AxesLabel->{"Mode","Energy"}]

kModes=Table[plotMomentum[vectors[[i]],DisplayFunction->Identity,PlotRange -> {0,0.4},Frame -> True,FrameLabel -> "Mode = "<>ToString[nGrid-i+1]<>" Energy = "<>ToString[values[[i]]]],{i,nGrid,nGrid-30,-1}];

Manipulate[Show[kModes[[n]]], {n,1,30,1}]

What happens as the potential well depth decreases? Try V = -0.05 and -0.02.

A Parabolic Well (A Quantum Harmonic Oscillator)

The shape of the potential, V(x), has a profound effect on the structure and energy of the quantum eigenstates.
In this section, we examine the quantum equivalent of a parabolic potential well, V(x) = (Subscript[V, 0])(x^2). We will center the potential well at the average grid location in our simulation.

vBarrier[i_Integer] := 0.025 (((i - (nGrid+1)/2))^2)

ListPlot[Table[vBarrier[i],{i,1,nGrid}],PlotStyle->PointSize[0.02]]

When the potential vanishes everywhere within the quantum box,

{values, vectors} = Eigensystem[Evaluate[mMatrix /. v[i_] -> vBarrier[i]// N]];

ListPlot[values,PlotStyle->{PointSize[0.02]},AxesLabel->{"Mode","Energy"}]

Mathematica orders the modes from the highest frequency (i.e. highest energy) to the lowest:

allModes=Table[ListPlot[((Abs[vectors[[i]]])^2),PlotRange->{{0,nGrid+1},{0,0.25}},Joined->True,PlotStyle->{Hue[((30-i)/5)],Thickness[0.01]},Frame->True,Ticks->None,FrameTicks->None,DisplayFunction->Identity, FrameLabel -> "Mode = "<>ToString[nGrid-i+1]<>" Energy = "<>ToString[values[[i]]]],{i,nGrid,nGrid-30,-1}];

Manipulate[Show[allModes[[n]]], {n,1,30,1}]

(These are Hermite polynomials, and Mathematica is able to find these solutions exactly.)

Table[plotψ[vectors[[i]],PlotRange-> {-0.7,0.7}],{i,nGrid,nGrid-5,-1}]

The parabolic potential well confines the quantum ball more or less to the center of the box.
The energy levels are approximately linearly dependent upon the mode index, provided that the grid spacing is sufficient to accurate represent the mode. (Compare this to the uniform box. In this case, the energy levels increase as the square of the mode index.)

ListPlot[values,PlotStyle->{PointSize[0.02]},PlotRange->{{30,50},{0,4}},AxesLabel->{"Mode","Energy"}]

kModes=Table[plotMomentum[vectors[[i]],DisplayFunction->Identity,PlotRange -> {0,0.15},Frame -> True,FrameLabel -> "Mode = "<>ToString[nGrid-i+1]<>" Energy = "<>ToString[values[[i]]]],{i,nGrid,nGrid-30,-1}];

Manipulate[Show[kModes[[n]]], {n,1,30,1}]

Multiple Modes

twoψ[t_] := (vectors[[nGrid-5]] Exp[-I values[[nGrid-5]] t] + vectors[[nGrid-6]] Exp[-I values[[nGrid-6]] t])/Sqrt[2]

ListPlot[((Abs[twoψ[0.0]])^2),Joined->True]

ListAnimate[Table[ListPlot[((Abs[twoψ[t]])^2),PlotRange->{{0,nGrid+1},{0,0.2}},Joined->True,Frame->True,Ticks->None,FrameTicks->None],{t,0,100,5}], AnimationRunning -> False]

ListAnimate[Table[plotψ[twoψ[t], PlotRange-> {-0.5, 0.5}],{t,0,100,5}], AnimationRunning -> False]

Summary

We have used Mathematica's built-in Eigensystem[…] command to analyze the normal modes of the quantum wavefunction. For a "particle in a box," the normal modes are the sine and cosine functions. For parabolic potential wells, they are Hermite polynomials. For potential wells, trapped particle modes appear, provided that the potential well is sufficiently deep.