Quantum Oscillations

AP1601

Quantum mechanics describes equations of motion of atomic and subatomic particles. At the small size of atoms, electrons and ions behave both as individual particles and as waves. Like marbles and balls, we can count atoms and electrons one-by-one. Like waves crashing onto the beach, electrons can constructively and destructively interfere with each other forming "crests and valleys" or interference patterns of particle current.

Although it is not possible to learn and to understand quantum mechanics in a few lectures, it is very easy to experiment with quantum mechanics using Mathematica. In order to proceed, you will have to accept a few (unusual) rules about the quantum nature of particles. This is not surprising. We have absolutely no first-hand experience with atomic particles. But, if you accept these few rules about quantum mechanics, you can use Mathematica to begin a very useful study of quantum behavior brought to you by computer simulation.

I have included below a short introduction to quantum mechanics adapted from an article written by the great Prof. Gerald Feinberg of Columbia's Physics Department. Following this introduction, I will review complex arithmetic, present the rules for the interpretation of quantum wave mechanics, and, finally, introduce a method for simulating quantum fluctuations with Mathematica.

A Brief Introduction to Quantum Mechanics

(Adapted from Quantum Mechanics by G. Feinberg.)

Quantum mechanics is the fundamental theory used by 20th- century physicists to describe atomic and subatomic phenomena. It has been successful in tying together a wide range of observations into a coherent picture of the universe.
While quantum mechanics uses some of the concepts of Newtonian mechanics, the previous description of physical phenomena, it differs fundamentally from the Newtonian description. For example, in Newtonian physics, quantities were believed to be continuously variable, able to take on any value in some range. An example is angular momentum, which, for a particle revolving in a circular orbit about some center of attraction, is proportional to the speed multiplied by the distance from the center. Because that distance could have any value in Newtonian mechanics, so could the angular momentum. In quantum mechanics, on the other hand, angular momentum is always restricted to certain discrete values, whose ratios are simple rational numbers.
An even more fundamental difference between quantum mechanics and previous physical theories is that probability enters in a basic way into how quantum mechanics describes the world. This is made evident by the ways that quantum mechanics and Newtonian mechanics deal with predictions of the future. For something described by Newtonian mechanics, such as the motion of planets about our Sun, it is possible, if sufficiently accurate measurements are made at one time, to predict the future behavior of the system to arbitrarily great accuracy. For systems described by quantum mechanics, even one as simple as an atom with a single electron, precise prediction of future behavior is usually impossible. Instead, only predictions of the probability of various behaviors can be made.

Quanta

Quantum mechanics was developed over a period of 30 years, during which it was successively applied to several physical phenomena. 
The first use of quantum ideas was made in the analysis of how electromagnetic radiation is produced. This was done by the German physicist Max Planck in 1900. Planck was trying to account for the distribution, among different frequencies, of the radiation emitted by a hot object, such as the surface of the Sun. He found that to obtain results in agreement with observation, he had to assume that the radiation was not emitted continuously, as was previously believed. Instead, it was emitted in discrete amounts, which he called quanta. For these quanta, there was always a relation between the frequency f, and the amount of energy emitted E, of the form E = hf. Here h is a universal constant introduced by Planck, and now named after him. Planck's constant has the units of energy multiplied by time, known as action. Its numerical value is approximately 6.63 (10^(-34)) joule-seconds. The specific result of Planck's analysis was a formula expressing the amount of energy radiated at any frequency as a function of the temperature of the emitting object. This relation, the blackbody distribution, agrees accurately with observation.

The development of actual quantum mechanics--the mathematical theory--took place in the years 1924-27. Initially there were two seemingly different approaches: matrix mechanics, invented by Werner Heisenberg, and wave mechanics, invented by Erwin Schrodinger. However, it was soon shown that these were distinct aspects of a single theory, which came to be known as quantum mechanics. This unified version was invented by Paul Dirac. In matrix mechanics, physical quantities such as the position of a particle, are represented not by numbers, but by mathematical quantities known as matrices. Matrix mechanics is most useful in dealing with situations in which there is a small number of relevant energy levels, such as a particle with definite angular momentum in a magnetic field.

Max Planck

The German physicist Max Karl Ernst Ludwig Planck, b. Apr. 23, 1858, d. Oct. 3, 1947, developed the concept of the quantum, or fundamental increment of energy--basic to quantum mechanics, and a cornerstone of modern physics. After receiving his Ph.D. from the University of Munich in 1879, Planck taught at the University of Kiel (1885-89) and the University of Berlin (1889-1926). His appointment at the latter institution included the directorship of the Institute of Theoretical Physics that was newly founded for him.
Planck began studying blackbody radiation in 1897 and discovered that at long wavelengths it did not obey the distribution laws given by Wilhelm Wien. This discovery led him to announce (1900) his revolutionary idea that an oscillator could emit energy only in discrete quanta, contrary to classical physical theory. The quantum theory--which gained Planck the Nobel Prize for physics in 1918--was used by Albert Einstein to explain (1905) the photoelectric effect and by Niels Bohr to propose (1913) a model of the atom with quantized electronic states; the theory was later developed into quantum mechanics. Planck mastered every aspect of physics from thermodynamics and electrodynamics to relativity and also wrote extensively on the philosophy of science.

The Wave Function

Wave mechanics is more useful than matrix mechanics in a situation where the number of energy levels is infinite, as with an electron in an atom. It is based on the idea originally suggested by Louis de Broglie, that particles such as electrons have waves associated with them. The wavelength, γ, of the wave is related to the mass, m, and speed, v, of the particle by the relation γ = h/mv. This implies that for electrons moving at 10 percent the speed of light, such as those produced by some television tubes, the wavelengths are about (10^(-10)) meters, or about the distance between atoms in a crystalline solid. This prediction of deBroglie was verified by Clinton Davisson and by George Thomson, who were able to pass the electron waves through metallic crystals, and so produce diffraction patterns similar to those produced by X rays.

In 1925, Erwin Schrodinger developed an equation, now bearing his name, which describes how a wave associated with an electron or other subatomic particle varies in space and time as the particle moves under the influence of various forces. This equation has many types of solutions, and Schrodinger imposed the condition that for a particle bound in an atom, the solution should be mathematically well defined everywhere. When applied to the case of an electron in a hydrogen atom, Schrodinger's equation immediately gave the correct energy levels previously calculated by Bohr. However, the equation could also be applied to more complicated atoms, and even to particles not bound in atoms at all. It was soon found that in every case, Schrodinger's equation gave a correct description of a particle's behavior, provided that the particle was not moving at a speed near that of light.
In spite of this success, Schrodinger initially misunderstood the meaning of the quantum wave function and believed that the intensity of the wave at a point in space represented the "amount" of the electron that was present at that point. That is, the electron was spread out, not concentrated at a point. However, this interpretation was wrong, because even if a particle was originally concentrated on a small region, in most cases it would soon spread over an increasingly larger region, in contradiction to what was the observed behavior of particles.

The correct interpretation of the waves was proposed by Max Born. While studying how quantum mechanics describes collisions between particles, he realized that the intensity of the deBroglie-Schrodinger wave was a measure of the probability of finding the particle at each point in space. In other words, a measurement would always find a whole particle, rather than a fraction of one, but in regions where the wave intensity was low, the particle would rarely be found, whereas in regions of high intensity, the particle would often be found.

Erwin Schrodinger

The Austrian theoretical physicist Erwin Schrodinger, b. Aug. 12, 1887, d. Jan. 4, 1961, published (1926) four papers that laid the foundation of the wave-mechanics approach to quantum theory and set forth his now-famous wave equation. Schrodinger earned a doctorate at the University of Vienna in 1910. He succeeded (1927) Max Planck in the chair of theoretical physics at the University of Berlin but left Germany in 1933 because of Nazi threats--the same year he shared the Nobel Prize for physics with Paul Dirac for his contributions to atomic theory. In 1939 he joined the newly formed Institute for Advanced Studies in Dublin. There he continued his studies of the application and statistical interpretation of wave mechanics, the mathematical character of the new statistics, and the relationship of these statistics to statistical thermodynamics.

The Schrodinger Equation

For our Mathematica simulations, we will use the symbol, ψ, to represent Edwin Schrodinger's wave function. ψ(x, t) is a complex function of space and time. Since ψ is complex, it's value at any point has both real and imaginary parts. We will solve Schrodinger's equation to see how the complex function ψ(x, t) changes in time and space.
The interpretation of the wave function is what gives ψ its usefulness. As Schrodinger wrote it, ψ is a function of space, or x. In this coordinate system, the product of ψ with its complex conjugate is proportional to a particle's probability distribution function at the location x. This probability must be normalized such that if we sum the probabilities at all possible locations, this sum must equal to 1, or 100%.

Momentum and Energy Representations

In order to determine the probability distribution functions for momentum or energy, we need to transform ψ into a different coordinate system (sometimes called a different "representation"). For momentum, the appropriate coordinate system is the wavenumber of the spatial oscillations of the ψ (i.e. the inverse of the wavelength). The shorter the wavelength of the oscillations of ψ, the higher is the particle's momentum. For energy, the appropriate coordinate system (or "representation") is the frequency of the temporal oscillations. The more quickly ψ oscillates in time, the higher the energy of the particle. These transformations can be performed by Fourier (or Laplace) transforms. If we were to Fourier transform the spatial coordinate of ψ(x), the result would be a new complex wave function, ψ(k), where k = 2π/λ. The probability of finding a particle of a particular momentum, p = ℏk, is Conjugate[ψ(k)] × ψ(k). (ℏ = h/2π.) If we were to Fourier transform the time variation of  ψ(t), the result would be a new complex wavefunction, ψ(ω). The probability of measuring a particle at a particular energy, e = ℏω, is proportional to Conjugate[ψ(ω)] × ψ(ω).

We can summarize the rules used to interpret Schrodinger's wavefunction:
1. The wavefunction, ψ(x,, t), is a complex function which can be used to calculate the probability of the measurements of particles. 
2. The probability of finding a particle at a location, x, at the instant of time, t, is proportional to ((Abs[ψ])^2) = Conjugate[ψ(x)] × ψ(x).
3. The probability of finding a particle having a particular momentum at an instant of time, t, is proportional to the square of the spatial Fourier transform of ψ. We write this probability as ((Abs[ ψ(k, t)])^2)= Conjugate[ψ(k)] × ψ(k).
4. The probability of measuring a particle having a particular energy at a particular location, x, is proportional to the square of the temporal Fourier transform of ψ. We write this probability as ((Abs[ ψ(x, ω)])^2)= Conjugate[ψ(ω)] × ψ(ω).

A very small ball in a quantum box

For simplicity, we will consider only a particle of mass, m, constrained to move only in one dimension, x, and contained within a box of length, l.

If our box were the size of a billiards table, the ball's motion could be accurately described without the need to use quantum mechanics. However, if the box were the size of an atom and if the ball were an electron, then quantum mechanics would be required to describe its motion. 
In the following, we will learn (1) how to compute the probabilities for measurements of the particle's position, x, momentum, p, and total energy, e, and (2) how to compute the evolution of these probabilities in time.

Schrodinger's equation on a grid

In general, the simulation of quantum behavior on a computer is difficult. Nevertheless, the same techniques developed for simulations of the oscillations of a string can be used for a quantum ball in a box.
Our first step is to define a grid on which to represent the wavefunction, ψ. For real particles, the coordinate system used to represent quantum oscillations is continuous. (The same can be said for oscillations on a string.) However, we can represent quantum oscillations approximately whenever the wavelength of the oscillations is long compared to the spacing between grid points.

If ψ(x, t) were a function of continuous variables x and t (i.e. the usual case), then Schrodinger's equation would be 
	I ℏ (∂_t) ψ = - ((ℏ^2)/2m) (∂_x) (∂_x) ψ + V(x) ψ
In this equation, V(x) is the potential energy of the quantum ball at x.
What does this equation mean? 
Schrodinger's equation defines the relationship between the spatial and temporal oscillations of the wavefunction. Recall our quantum rules: the spatial wavelength is inversely proportional to the momentum and the temporal frequency is proportional to the energy. If we were to substitute for the wavefunction, 
	ψ(x, t) ∝ Exp[- I ω t] Exp[ I k x],
then Schrodinger's equation can be rewritten as
	ℏω = (((ℏk))^2)/2m + V
This is a statement which defines the total energy, e = ℏω, in terms of the momentum, p = ℏk, and the potential energy, V.

Schrodinger's equation on a grid requires that we replace the spatial derivatives, (∂_x) (∂_x) ψ, 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 ith grid points are:
	I ℏ (∂_t) ψ[1][t]  = - ((ℏ^2)/2m) (- 2 ψ[1][t] +  ψ[2][t]) + V[2] ψ[2][t]
	I ℏ (∂_t) ψ[i][t]  = - ((ℏ^2)/2m) ( ψ[i-1][t] - 2 ψ[i][t] +  ψ[i+1][t]) +V[i] ψ[i][t]
	I ℏ (∂_t) ψ[N][t]  = - ((ℏ^2)/2m) ( ψ[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 away from the ends.
Compare these equations to those we used to simulate oscillations on a string. There are two significant differences and several similarities. The differences are: (1) the quantum wave function is complex, and (2) the quantum equations are first order in time (whereas the equations for a string are second-order).

Finite Difference Approximations

For this topic, we will not use NDSolve[...]. Instead, we will integrate the wave-dynamics using various approximations to the time derivative, (∂_t) ψ. We will examine three methods. For simplicity, we will assign ℏ -> 1 and m -> 1, and we will also set our length and time scales such that δt = δx = 1.

Explicit

The most straight-forward finite-difference approximation to Schrodinger's equation 
	I (∂_t) ψ = - (∂_x) (∂_x) ψ + V(x) ψ
is
	I (ψ(t + 1, x) - ψ(t, x)) ≈ - (ψ(t, x-1) - 2 ψ(t, x) + ψ(t, x+1)) + V(x) ψ(t, x)
We can find ψ(t+1, x) from ψ(t, x):
	ψ(t + 1, x) ≈ ψ(t, x) + I (ψ(t, x-1) - 2 ψ(t, x) + ψ(t, x+1)) - I V(x) ψ(t, x)

Implicit

If the right-hand-side were evaluated at the advanced time, then Schrodinger's equation would be approximately:
	I (ψ(t + 1, x) - ψ(t, x)) ≈ - (ψ(t+1, x-1) - 2 ψ(t+1, x) + ψ(t+1, x+1)) + V(x) ψ(t+1, x)
We must invert a matrix to find ψ(t+1, x) from ψ(t, x):
	ψ(t + 1, x)  - I (ψ(t+1, x-1) - 2 ψ(t+1, x) + ψ(t+1, x+1)) + I V(x) ψ(t+1, x)  ≈  ψ(t, x)

Centrally-Balanced (or Crank-Nicolson)

Finally, we can take the average of both times to find the right-hand-side:
	I (ψ(t + 1, x) - ψ(t, x)) ≈ [- (∂_x) (∂_x) + V(x)] (1/2)(ψ(t+1, x) + ψ(t, x))
Arranging terms, we find an approximate relationship:
	ψ(t + 1, x)  - I (ψ(t+1, x-1) - 2 ψ(t+1, x) + ψ(t+1, x+1))/2 + I V(x) ψ(t+1, x)/2  ≈  
		ψ(t, x) + I (ψ(t, x-1) - 2 ψ(t, x) + ψ(t, x+1))/2 - I V(x) ψ(t, x)/2

We will find that only the last method will be satisfactory. The explicit method is unstable. The implicit method is inaccurate. And, the Crank-Nicolson method is both numerically efficient and accurate.

Initial Condition

We choose an initial condition with a Gaussian spatial extent and a complex variation corresponding to a traveling momentum.

form[k_,x0_,s_][x_] := Exp[I k x] Exp[-((x-x0)/(2 s))^2];

The integral of the norm must be normalized to unity.

myConjugate[x_] := x /. Complex[re_,im_] -> Complex[re,-im]

myConjugate[form[k,x0,s0][x]]

Integrate[form[k,x0,s0][x]  myConjugate[form[k,x0,s0][x]], {x, -Infinity, Infinity}]

initial[k_,x0_,s_][x_] := Exp[I k x] Exp[-((x-x0)/(2 s))^2]/Sqrt[ Sqrt[2 π s^2]];

Notice: In our units, the spatial grid size are integers.

ψ0=Table[initial[0.6,0.0,8.0][x],{x,-50,50}];
ListPlot[Abs[ψ0],PlotRange->All]

We need to check that the initial condition integrates to a unity norm (or 100% probability of being located within the "box"):

Plus @@ (ψ0 Conjugate[ψ0])

Plotting the WaveFunction

p0=0.6;
x0=0.0;
s0=8.0;
pC=Plot[{Re[initial[p0,x0,s0][x]],Im[initial[p0,x0,s0][x]]},{x,-50,50},PlotStyle->{Blue,Red}];
pA=Plot[{Abs[initial[p0,x0,s0][x]],-Abs[initial[p0,x0,s0][x]]},{x,-50,50}];
Show[pA,pC,PlotLabel->"Initial Condition"]

Since we will be working with ψ(x, t) on a grid, we define a useful plotting utility:

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]]

plotψ[ψ0, Frame-> True]

How big "really"?

Can we put some real units on the dimensions?

Needs["PhysicalConstants`"];

m = ElectronMass

PlanckConstantReduced

An electron moving 100,000 m/s = 100 km/s (about room temperature)... E= ℏω and p = ℏ k.

energy = (1/2)m (1.0*^5 Meter /Second)^2  /.Kilogram  -> Joule Second^2 /Meter^2

%/BoltzmannConstant

frequency = energy/PlanckConstantReduced /(2 π)

momentum = m (1.0*^5 Meter / Second)/.Kilogram  -> Joule Second^2 /Meter^2

wavenumber = momentum/PlanckConstantReduced

wavelength = 2 π / wavenumber

For our initial condition, the wavenumber is p = 0.5.  This sets the length scale. The wavelength is about 2 π × 2 ~ 12 grid points.

Schrodinger's Equation: Explicit Method

The explicit method is the most straight-forward finite-difference approximation to Schrodinger's equation. However, as we will see, this method is unstable and it can not be used. For the explicit method, we find ψ(t+1, x) from ψ(t, x) using the formula below:
	ψ(t + 1, x) ≈ ψ(t, x) + I (ψ(t, x-1) - 2 ψ(t, x) + ψ(t, x+1)) - I V(x) ψ(t, x)

We'll define our grid size and our initial condition:

nGrid = 100;

ψ0 = Table[initial[0.5,nGrid/2,5.0][x], {x,1,nGrid}];
plotψ[ψ0, PlotRange-> All, Frame-> True]

If we let the potential energy be V = 0, then each explicit time-step of the Schrodinger's equation on a grid can be written using ListConvolve[…]. (This is very efficient and  fast.)

Using ListConvolve[…]

ListConvolve[{c,b,a},Table[ψ[i],{i,1,10}], 2, 0] // MatrixForm

We use "0" as padding. This represents the hard "walls" at each side of the computational grid.

Explicit Step

eStep[ψ_List, v_List] := ListConvolve[{I 0.5, 1.0 - I , I 0.5},ψ, 2,0.0] - I ψ v

eStep[Table[ψ[i], {i, 1, 10}], Table[0.0, {i, 1, 10}]] //Chop // MatrixForm

Sample Explicit Solution

p0 = 0.6;
x0 = nGrid/2;
s0 = 8.0;
ψ0 = (Table[initial[p0,x0,s0][x], {x, 1, nGrid}] // N);
v0 = Table[0.0, {x, 1, nGrid}];
plotψ[ψ0, PlotRange-> All, Frame-> True]

ψ1 = (eStep[ψ0,v0]; // Timing)

ψ10 = Nest[eStep[#,v0]&, ψ0, 10];
plotψ[ψ10, PlotRange-> All, Frame-> True]

sol = NestList[eStep[#,v0]&, ψ0, 100];

ListDensityPlot[Abs[sol],Mesh->False,ColorFunction->Hue,PlotRange->{0,0.3}]

ListPlot3D[Abs[sol],PlotRange->{0,0.3},Mesh->False]

This scheme goes unstable after about 35 iterations with δt = 1.0.

Schrodinger's Equation: Implicit Method

If the finite-difference approximation of the time-derivative (on the right-hand-side) were evaluated at the advanced time, then Schrodinger's equation is approximated by an implicit formula. This method is stable, but it can not be used because it introduces "numerical diffusion". The sum of the probabilities do not remain at 100%. 
For the implicit method, we must invert a matrix to find ψ(t+1, x) from ψ(t, x):
	ψ(t + 1, x)  - I (ψ(t+1, x-1) - 2 ψ(t+1, x) + ψ(t+1, x+1)) + I V(x) ψ(t+1, x)  ≈  ψ(t, x)

?LinearSolve

Tridiagonal

In order to solve this equation efficiently,we need to invert a "tridiagonal matrix". A tridiagonal matrix is a "sparse array", with zeros everywhere except near the diagonal.

?SparseArray

p = Table[a[i], {i, 1, 5}]

sparseA = SparseArray[{{i_,i_} :> 1 + I (1.0 + p[[i]]), {i_,j_}/; Abs[i-j] ==  1 -> - I 0.5}, {5, 5}, 0.0]

MatrixForm[sparseA]

iStep[ψ_List,v_List]:= Module[{sparseA,len, sol},
len = Length[ψ];
sparseA = SparseArray[{{i_,i_} :> 1 + I (1.0 + v[[i]]), {i_,j_}/; Abs[i-j] ==  1 -> - I 0.5}, {len, len}, 0.0];
(* fixed boundaries *)
sol = LinearSolve[sparseA,ψ];
sol]

Sample Implicit Solution

p0 = 0.6;
x0 = nGrid/2;
s0 = 8.0;
ψ0 = (Table[initial[p0,x0,s0][x], {x, 1, nGrid}] // N);
v0 = Table[0.0, {x, 1, nGrid}];
plotψ[ψ0, PlotRange-> All, Frame-> True]

(ψ1 = iStep[ψ0,v0];) // Timing

plotψ[ψ1, PlotRange-> All, Frame-> True]

After 10 steps:

ψ10 = Nest[iStep[#,v0]&, ψ0, 10];
plotψ[ψ10, PlotRange-> All, Frame-> True]

sol = NestList[iStep[#,v0]&, ψ0, 100];

ListDensityPlot[Abs[sol],Mesh->False,ColorFunction->Hue,PlotRange->{0,0.3}]

ListPlot3D[Abs[sol],PlotRange->{0,0.3},Mesh->False]

This scheme is stable, but it does not preserve the norm, or the probability of the quantum particle.

Is Probability Preserved?

norm[ψ_List] := Plus @@ (ψ Conjugate[ψ]);

norms = (norm /@ sol) // Chop;

norms // Shallow

ListPlot[norms]

Schrodinger's Equation: Crank-Nicolson Method

Finally, we can take the average of both times to find the right-hand-side:
	I (ψ(t + 1, x) - ψ(t, x)) ≈ [- (∂_x) (∂_x) + V(x)] (1/2)(ψ(t+1, x) + ψ(t, x))
Arranging terms, we find an approximate relationship:
	ψ(t + 1, x)  - I (ψ(t+1, x-1) - 2 ψ(t+1, x) + ψ(t+1, x+1))/2 + I V(x) ψ(t+1, x)/2  ≈  
		ψ(t, x) + I (ψ(t, x-1) - 2 ψ(t, x) + ψ(t, x+1))/2 - I V(x) ψ(t, x)/2

Notice that this method combines the operations used in both the explicit and implicit methods. We will find that the Crank-Nicolson method is both numerically efficient and accurate.

cnStep[ψ_List, v_List] := Module[{sparseA,r,len},
len = Length[ψ];
sparseA = SparseArray[{{i_,i_} :> 1 + I (1.0 + v[[i]]) 0.5, {i_,j_}/; Abs[i-j] ==  1 -> - I 0.25}, {len, len}, 0.0];
r =  ListConvolve[{I 0.25,1.0 - I 0.5, I 0.25},ψ, 2, 0.0] - I ψ v 0.5;
(* fixed boundaries *)
sol = LinearSolve[sparseA,r];
sol]

Sample Crank-Nicolson Solution

p0 = 0.6;
x0 = nGrid/2;
s0 = 8.0;
ψ0 = (Table[initial[p0,x0,s0][x], {x, 1, nGrid}] // N);
v0 = Table[0.0, {x, 1, nGrid}];
plotψ[ψ0, PlotRange-> All, Frame-> True]

(ψ1 = cnStep[ψ0,v0];) // Timing

After 10 steps:

ψ10 = Nest[cnStep[#,v0]&, ψ0, 10];
plotψ[ψ10, PlotRange-> All, Frame-> True]

sol = NestList[cnStep[#,v0]&, ψ0, 100];

ListDensityPlot[Abs[sol],Mesh->False,ColorFunction->Hue,PlotRange->{0,0.3}]

ListPlot3D[Abs[sol],PlotRange->{0,0.3},Mesh->False]

This scheme is stable, but it does not preserve the norm, or the probability of the quantum particle.

Is Probability Preserved?

norm[ψ_List] := Plus @@ (ψ Conjugate[ψ]);

norms = (norm /@ sol) // Chop;

norms // Shallow

ListPlot[norms-1.]

and this is about as perfect as numerical approximations can get!

Longer Time Solution

(sol = NestList[cnStep[#,v0]&, ψ0, 400]; // Timing)

ListDensityPlot[Abs[sol],Mesh->False,ColorFunction->Hue,PlotRange->{0,0.3}]

ListPlot3D[Abs[sol],PlotRange->{0,0.3},Mesh->False]

Numerical errors accumulate over time, but the accuracy is still very good.

norms=Chop[norm/@sol];
ListPlot[norms-1.]

Other Graphics and Animations

plotψ[sol[[1]], PlotRange-> {-0.25,0.25}]

plotψ[sol[[100]], PlotRange-> {-0.25,0.25}]

ListAnimate[Table[plotψ[sol[[i]],Frame-> True, PlotRange-> {-0.25,0.25}], {i, 1, 400,4}], AnimationRunning -> False]

Summary

We have investigated three methods to solve the time-dependent Schrödinger's equation in one-dimension using finite-difference approximations. We find explicit methods do not work since the solution becomes numerically unstable. We find fully implicit methods do not work because they are inaccurate. Fully implicit methods introduce "numerical diffusion". When we approximate the right-hand-side by the average of the present and forward-time values, then the solution is well-behaved and accurate. This is the method is called the Crank-Nicolson method.

The forms of the solutions for a particle trapped within a box depend on the initial conditions. We looked at one initial condition showing a particle bouncing off the sides of the box. As a consequence of the "wave-particle duality" of quantum mechanics, an interference patter is seen when the particle strikes the wall. The leading-edge of the reflected particle interferes with the trailing edge of the incident particle. The wave-length of the interference pattern depends on the particle's "wavelength" or momentum.