Chem 3111                                                                                                               Eric R. Bittner
Univ. of Houston                                                                                                  

Thanks to Hannes Jónsson and Tom Marchioro, U. Washington

Making Plots in Mathematica

In this exercise we will learn to make plots using Mathematica.  We will examine the quantum states of a particle in a one and two dimensional box.  If you haven't had quantum mechanics, don't worry too much.  We're primarily interested in making pictures of the wavefunctions rather than learning to solve the equations.  


A. Particle in a 1-D box:

The time-independent Schroedinger equation in one dimension is
                -h^2/(8 mπ^2)(d^2 ψ(x))/dx^2 + V(x) ψ(x) = E ψ(x)           
where ψ(x) is the  wavefunction, h is Planck's constant, m is the mass of the particle and E is  the total energy.  For a one-dimensional box, where the potential energy is  V(x) = 0 within the region 0 < x < L and infinity outside the region, the  wavefunction must vanish at the edges of the box, ψ(0) = 0 and ψ(L)  = 0, giving two boundary conditions which the wavefunction must also satisfy.   The solutions represent stationary states (analogous to the 'standing waves'  of a string) for the particle.  Solutions to both the differential equation  and the two boundary conditions can only be found for certain values of the  energy, E. The normalized wavefunctions are
                  ψ(x) = 2/L^(1/2)sin(nπx/L)
where n is an integer, n=1, 2, 3, ... The probability of  finding the particle in different regions within the box is related to ψ^2(x), in such a way that the total probability of finding the particle  within the region a < x < b is
                  Pab = ∫ _ a^b ψ^2 d x

We will use Mathematica to plot the first few stationary state wavefunctions and probability  distributions, and calculate the probability of finding the particle in  various regions within the box.  First, define a function corresponding to  the wavefunction. It is good practice to clear the symbols that will be used for the function and the variables

Clear[psi, x, n] psi[x_, n_] := Sqrt[2/L] * Sin[n * Pi * x/L]

The factor 2/L^(1/2) ensures that the total probability of finding the particle within the box  is unity. This can be verified by integration, for example for the n=1 state

Integrate[psi[x, 1]^2, {x, 0, L}]

Q:   Do the integration for states, for example with quantum number n=2, n=3  and n=10.  Is the result unity in all cases?


In order to plot the wavefunctions, it is necessary to assign numerical values to the parameters.  We can choose units such that the mass is unity, Plancks constant is unity, and the box length is unity

L = 1 ; m = 1 ; h = 1 ;

Then, a plot of the ground state wavefunction and the probability distribution is easily obtained

Plot[{psi[x, 1], psi[x, 1]^2}, {x, 0, 1}]

It is nice to add color here to distinguish between the two curves.  This is done by specifying PlotStyle. Also, it can be nice to add a frame by setting the option Frame->True

Plot[{psi[x, 1], psi[x, 1]^2}, {x, 0, 1}, Frame -> True, PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 0, 1]}}]

Q:   Which curve is the wavefunction and which curve is the probability  distribution?

It is clear from this plot that when the particle is in the ground state it is most likely found near the center of the box. The probability of finding the particle in the region from x = 0.4 to x = 0.6, an interval that spans one fifth of the box length, can be found by integration of the probability density

Integrate[psi[x, 1]^2, {x, 0.4, 0.6}]

Q:  Calculate the probability of finding the particle in the region from x  = 0.0 to x = 0.2 (a region of the same length as the one above but located at  the edge of the box).  How much smaller is the probability of finding the  particle in the region [0.0, 0.2] as compared with the region [0.4, 0.6]?  (copy and paste the integration command above).

Q:  Consider a classical particle, i.e. a particle that is described by  Newton's law, F=ma.  What would the motion of the classical particle in the  box be like?  How would the probability of finding the classical particle in  the region  [0.0, 0.2] compare with the probability of finding it in the  region [0.4, 0.6]?  (just give a qualitative answer).

Q:   Make a plot of the wavefunction and the probability density for the  first excited state.  Calculate the probability of finding the particle in  the region [0.0, 0.2] and in the region [0.4, 0.6] and compare (copy and  paste the integration command from above).

Q:   Make a plot of the wavefunction and the probability density for the  n=10 state.  Calculate the probability of finding the particle in the region  [0.0, 0.2] and in the region [0.4, 0.6] and compare.  How many nodes does the  wavefunction have?

The energy of the particle in a box increases as n^2 .   It is given by
                                     E _ n = (h^2 n^2)/(8 mL^2)
The higher the energy is, the less important the wave properties of the  particle become because the wavelength becomes shorter, and the classical  description of the probability distribution becomes better.  


Q:  How is this reflected in the probability calculations above? (How does  the probability finding a highly excited particle in the interval [0.0,0.2]  and compare with the interval [0.4,0.6], and how does this compare with the  classical prediction).


B. Particle in a 2-D box:

The Schroedinger equation in two dimensions is
                        -h^2/(8 mπ^2)(∂^2 ψ(x, y)/∂ x^2  + ∂^2 ψ(x, y)/∂ y^2) + V(x,y) ψ(x,y) = E ψ(x,y)
Since the wavefunction depends on  two variables now, x and y, the second derivatives are partial derivatives  which are writen with ∂ rather than d.

For a particle in a two-dimensional box, the potential is taken to be 0 in  the region  where 0 < x < a  and  0 < y < b  and infinite everywhere else.  The potential can be written as a sum of a 1D "box potential" in the x  direction and a 1D "box potential" in the y direction.  When a  two-dimensional potential function can be written as a sum of two potential  functions each depending on only one of the Cartesian coordinates, then the  Hamiltonian separates and the wavefunction simply becomes the product of the  two one-dimensional wavefunctions.  Therefore, in the case of a  two-dimensional particle in a box, the wavefunctions are
                              ψ(x,y) = 2/a^(1/2)sin(n _ xπx/a)2/b^(1/2)sin(n _ yπx/b)
where n _ x and n _ y are quantum numbers, n _ x=1,2,..., and n _ y=1,2,...  The energy of the particle is      
                                 E _ (n _ x n _ y) = h^2/(8 m) (n _ x^2/a^2+n _ y^2/b^2)

For certain box shapes, where the lengths a and b are in just the right  ratio, the same value of energy can be obtained from two combination of  quantum numbers n _ x and n _ y, that is, two distinct quantum states have the same energy.  In such a  case it is said that the energy level is degenerate.

It is instructive to plot the wavefunctions for the two dimensional particle in a box. First, it is good to make sure the variables involved do not already have values assigned to them.

Clear[nx, ny, x, y, a, b, h, m] ;

Then, define the wavefunctions (a product of two one-dimensional particle-in-a-box wavefunctions)

psi2d[x_, nx_, y_, ny_] := (2/a)^(1/2) * Sin[nx * Pi * x/a] * (2/b)^(1/2) * Sin[ny * Pi * y/b] ;

and the energy levels

en[nx_, ny_] := (h^2/(8 * m)) * (nx^2/a^2 + ny^2/b^2)

The wavefunctions can be plotted using the built-in Mathematica three-dimensional plotting routine. Remember, values need to be assigned to all parameters in the expression for the wavefunction before a plot can be made.  For example, choosing the box to be rectangular with one side twice as long as the other,

a = 1 ; b = 2 ; h = 1 ; m = 1 ;

The plot the lowest energy state, the ground state, can be obtained by

Plot3D[psi2d[x, 1, y, 1], {x, 0, a}, {y, 0, b}] ;

This command could be repeated for each of the different functions, but it  would be tedious typing this in over and over, or else having to change the  values of n _ x and n _ y by hand. Instead, the process can be automated with Mathematica's Table function in the following way

graphs = Table[Plot3D[psi2d[x, nx, y, ny], {x, 0, a}, {y, 0, b}], <br /> {nx, 1, 3}, {ny, 1, 3}] ;

Mathematica now loops through nine combinations of values for n _ x and n _ y. The 3D plots are all stored in the variable 'graphs', so they can be  recalled and displayed in a variety of ways for comparison. For example,

Show[graphs[[3, 2]]] ;

shows the next to the last plot only, while

Show[GraphicsArray[graphs]] ;

puts all nine plots together in an array (you will probably need to enlarge the graph to see it clearly, select the graph and drag out the lower right corner).

Q:  Find out which pair of quantum numbers, n _ x and n _ y, each of the nine graphs corresponds too. Present your results as a 3x3  table of pairs of quantum numbers    (n _ x,n _ y)  (n _ x,n _ y)  (n _ x,n _ y)
                       (n _ x,n _ y)  (n _ x,n _ y)  (n _ x,n _ y)
                       (n _ x,n _ y)  (n _ x,n _ y)  (n _ x,n _ y)

in a text cell.
A similar approach can be taken to study the energy levels and their degeneracy. The command

energy = Table[en[nx, ny], {nx, 1, 5}, {ny, 1, 5}]

gives a two dimensional table (a matrix) of the energy of the first 25 quantum states. Because the ratio of a and b was set to 2 there are degenerate energy levels as can be seen by inspection of the matrix

MatrixForm[energy]

or, better yet, by sorting the energies

Sort[Flatten[energy]]

Q: Examine the above functions for larger values of the quantum numbers n _ x and n _ y. How many of the first 50 energy levels are degenerate?

Converted by Mathematica  (August 29, 2002)