SCHRODINGER EQUATION FOR BOUND STATES IN ONE DIMENSION

William M. MacDonald@2001

Mathematica for Scientists and Engineers

Initialization

Clear all previously defined user symbols.

Clear["Global`*"]

Turn off warnings about user symbols that resemble, but are not the same as, each other or system symbols.

Off[General::spell, General::spell1];

Use the Mathematica Browser to go to Section 3.10.3 of The Mathematica Book to see keyboard combinations for Greek letters and mathematical or physical symbols. Or, use File=>Palettes in the menu to open the CompleteCharacters palette.

Introduction

We want to explore the bound states of the one dimensional Schrodinger equation.  We use the standard notation of Mathematica for functions (f[x]) because it removes the ambiguity associated with parentheses in standard mathematical notation.

\(TraditionalForm\`\(-\(ℏ\^2\/\(2\ m\)\)\) d\^2 u[x]\/dx\^2 + V[x] u[x] == 
          E\ u[x]\)

We choose a natural scaling parameter a that has the dimension of length as the unit of length.  This is equivalent to using the dimensionaless length variable

s == s/a .

We choose as the unit of energy the quantity

((ℏ^2)/(2 m (a^2))).

It then has the numerical value unity in the equations that we write.  
	In treating one-dimensional wells we choose the ``width'' of the well as the unit of length.  In these units the Schrodinger equation for the wave function u becomes

u''[s]+(e-v[s])u[s]==0

where e is the dimensional energy and v[s] is the potential in these units.
	It is useful to solve this equation as a pair of coupled first order equations because we will need the derivative of the wave function u'[x] as well as the solution u[x] itself.  This is done by the usual method of replacing a second order differential equation by two coupled first order differential equations.

u'[s]== w[s]
w'[s]== -(e-v[s])u[s]

For bound states the energy will be negative and so we set

e == -(b^2)

giving the equation for u

u''[s]== ((b^2) + v[s]) u[s]

and the equivalent coupled equations

u'[s] == w[s];
w'[
              s] == ((b^2) + v[s]) u[s]

We will start solutions from some point distant from the potential to the left and integrate to some point distant to the right.   At points distant from the potential, where it may be approximated as zero, the solution of the equation for u is

u[s] = 
      D (e^(-b s)) + 
        C (e^(b s))

To the distant left, where v is zero, the solution must decay as s approaches negative infinity, that is,

u[s] = C (e^(b s)) .

To the right of the potential, the solution must decay as s approaches positive infinity, that is

u[s] = 
      D (e^(-b s))

Let us choose symmetrical points about the center of the potential as points for matching the interior solution to the exponentials outside the well.

Infinite Square Well

It is always useful to have some idea of the answer to a problem before one begins.The energy levels and wave functions for an infinite well are very easy to find because the wave function must vanish at the well edges.  We shall always choose the length scale so that the well boundaries occur at s= -1/2 and s=1/2. 
	With an infinite well we choose the zero of energy as the bottom of the well.  The potential is constant between the walls and therefore equal to zero. 
	The Schrodinger equation in dimensionless units becomes

u''[x] == - e u[x]

where e is equal to the kinetic energy in dimensionless units.

Wave functions

The wave functions for the square well must vanish at the edges. The solutions of this equation that vanish at s=±1/2 are

Clear[u];
u[s_, n_] = Cos[(2  n + 1) π s];

The norm of the wave function is the integral of the square.

Clear[norm]; 
    norm[n_] := 
      Simplify[(∫_(-1/2)^(1/2))u[x, 
                  n]^2 \[DifferentialD]x, 
        n ∈ Integers]

Normalize the functions to unity.

u[x_, n_] = u[x, n]/norm[n]^(1/2)

Use = (Set) to reassign values to a symbol. Using := (SetDelayed) will put the kernel into an infinite loop.

Plot the wave functions.

Plot[
        Evaluate[Table[u[s, n], {n, 0, 3}]], {s, -1/2, 
          1/2}];

Energy values

The Schrodinger equation for this infinite square well is

\(\(SchEq = ∂\_{s, 2}u[s, n] == \(-e\)\ u[s, n];\)\)

Solve for the e values, which are equal to energies relative to the bottom of the well.

Levels[n_] = 
      e /. Solve[SchEq, 
            e][[1]]

The lowest four levels are then

Levels[Range[0,3]]//N

Plot these levels by using this form for the potential.

v[x_ /; -1/2 ≤ s ≤ 1/2] := 
        0;
v[x_ /; s < -1/2 || s > 1/2] := 500;

Plot[
        Evaluate[Prepend[u[s, Range[0, 3]], v[s]]], {s, -0.502, 0.502}, 
        Axes -> None, 
        PlotStyle -> 
          Table[{Thickness[0.01], Hue[1 - s, 1, 1]}, {s, 0.2, 1, 
              0.2}]];

Finite Square Well

Consider the case of a finite square well.

v (s)=- v0       -1/2≤s≤1/2
            \(\(=\)\(\ \)\(0\ \ \ \ \ \ \ \ \ \ \ \ s < \(-\(\(/\)\(2\)\)\)\  || s > 1/2\)\)

The wave function must satisfy the equations

u''[x] == ((b^2) - v0) u[x]         - 
        1/2 ≤ s ≤ 1/2
\(\(u''\)[x] == 
        b\^2\ u[x]\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ s < \(-\(\(/\)\(2\)\)\)\  || s > 1/2\)

Since the energy e will be negative for bound states, we shall use the parameter b defined in the introduction
	e == -(b^2).

Clearing symbols

Clear["Global`*"]

Setting up the well

We shall choose a specific finite value for the depth of the well of v[x]

Clear[v];
v[
          s_ /; -0.5 ≤ s ≤ 
              0.5] := -v0;
v[s_ /; s < 0.5 || s > 0.5] := 0;

Here is the plot for a well of depth 40..

PotentialPlot = 
        Plot[v[s] /. v0 ->  40, {s, -2, 2}, 
          PlotRange -> {{-2, 2}, {-40, 10}}, 
          Frame ->  False, FrameLabel -> {"s", "v(s)"}, 
          PlotLabel -> "Plot of Potential", 
          PlotStyle -> Hue[0.7, 1, 0.8]];

Setting up the wave function

Outside the well and to the left, the wave function decays exponentially as s assumes large negative values.

uL[s_] := (Subscript[A, 1]) (E^(b s));

Inside the well the wave function is a linear combination of the two solutions to the Scrodinger equation.

uC[s_] := 
        (Subscript[A, 2]) Sin[(Sqrt[v0 - (b^2)]) s] + 
          (Subscript[B, 2]) Cos[(Sqrt[v0 - (b^2)]) s];

To the right of the well, the wave function decays exponentially as x increases.

uR[s_] := (Subscript[A, 3]) (E^(-b s))

Energy Levels

Eigenvalue equation

The boundary conditions on the wave function are the following:
  1. continuity of the wave function at s=±1/2
  2. continuity of the derivative of the wave function at s=±1/2
	The overall normalization of the wave function is determined by one parameter, (Subscript[A, 1]), so the task is to solve for (Subscript[A, 2]), (Subscript[B, 2]), and (Subscript[C, 3]) in terms of (Subscript[A, 1]).  There are four boundary conditions, so these three parameters are overdetermined.  This means that solutions can exist only for specific values of the parameter b, and these values correspond to the energies of the bound states.  These are the energy levels of the quantum well.
	We first use the boundary conditions to eliminate (Subscript[A, 1]) and (Subscript[A, 3]), obtaining two homogeneous equations for (Subscript[A, 2]) and (Subscript[B, 2]). The condition for solutions will be that the determinant of these equations vanish.  
	The constants (Subscript[A, 1]) and (Subscript[A, 3]) are most easily eliminated by recognizing that the continuity of a wave function u[x] and its derivative u'[x] implies that the ratio u'[x]/u[x] is continuous.  This is called the ``log derivative'' because it is just the derivative of the logarithm of u[x].  Any multiplication factor of an expression drops out of the log derivative, so both (Subscript[A, 1]) and (Subscript[A, 3]) cancel out of the log derivative of uL[x] and uR[x] respectively. 
	The boundary conditions for the square well can be applied at its boundaries because it vanishes identically outside.  The expressions containing (Subscript[A, 2]) and (Subscript[B, 2]) that result from matching the log derivatives at the left and the right are easily found.

eq[1] = (Numerator[uL'[s]/uL[s] - uC'[s]/uC[s] /. 
                s -> -1/2 // Together] // Simplify)

eq[2] = (Numerator[uR'[s]/uR[s] - uC'[s]/uC[s] /. 
                s -> 1/2 // Together] // Simplify)

We now use the symbol Coefficient to pick out the coefficients of (Subscript[A, 2])and (Subscript[B, 2])on the left side of each equation.  The result is a matrix, call it m.

Clear[m]; 
    m[b_] = {Coefficient[eq[1], {(Subscript[A, 2]), (Subscript[B, 2])}], 
        Coefficient[eq[2], {(Subscript[A, 2]), (Subscript[B, 2])}]}

The condition for solutions of eq[1] and eq[2] is that

Det[m[b]]==0,

that is, the determinant of the matrix must vanish.  This determinant is a function of the parameter b, and the bound states occur where it vanishes.  
	Let us define

EigenEq[b_] = Det[m[b]];

Eigenstates for v0=40

We now set v0=40

v0 = 40;

Plot[EigenEq[b],{b,0,Sqrt[v0]}];

There appear to be three values of b where the eignenvalue equation is satisfied. Go to the menu and choose Input:Get Graphics Coordinates to learn how to read them from this graph and save them in the following cell.

bvalues = {{0.0920886, 1.24769}, {4.24512, 1.24769}, {5.82127, 
            0.735651}, {6.32163, 0.735651}};

There are therefore four values of b, corresponding to four bound energy levels.

Use FindRoot to obtain values of b that satisfy EigenEqn[b].

?FindRoot

EigenValues=Map[FindRoot[EigenEq[b],{b,#[[1]]}]&,bvalues]

Is there an eigenvalue near 6.32456?

Table[{b, Re[EigenEq[b]], Im[EigenEq[b]]}, {b, 6.3, (Sqrt[v0]) + 0.1, 
            0.01}] // N // TableForm

Although it appears that the function EigenEq[b] vanishes beyond b=6 , it is only the real part that vanishes because EigenEq[b] becomes purely imaginary when the square of b equals or exceeds v0. So we drop the last value.

EigenValues = Drop[EigenValues, -1]

The energies are

BindingEnergies = -(b^2) /. EigenValues

Compare these values with those for the infinite square well by computing their distance from the bottom of the v0=40 well.  Sort from lowest to highest.

Sort[BindingEnergies + v0]

Here are the levels for the infinite square well.

Wave Functions for Bound States

Relate A2,B2, and A3 to A1

We now find the wave functions for each of the bound states by returning to the boundary conditions and solving for A2,B2, and A3 in terms of A1.  The value of A1 will then be determined by normalizing the wave function, that is, requiring that the integral of the square of the wave function over all values of s be unity.

Clear[A2B2Sol];
A2B2Sol = 
        Solve[{(∂_s) uL[s] == (∂_s) uC[s], 
              uL[s] == uC[s]} /. s -> -0.5, {(Subscript[A, 2]), 
            (Subscript[B, 2])}];

{(Subscript[A, 2]), 
          (Subscript[B, 2])} = {(Subscript[A, 2]), 
            (Subscript[B, 2])} /. A2B2Sol[[1]];

A3Sol = 
        Solve[uC[s] == uR[s] /. s -> 0.5, 
          (Subscript[A, 3])];

(Subscript[A, 3]) = 
        (Subscript[A, 3]) /. A3Sol[[1]];

Wave function for states

The unnormalized wave functions in the well are

uC[s_, 1] = 
        Chop[uC[s] /. EigenValues[[1]]];
uC[s_, 2] = 
        Chop[uC[s] /. EigenValues[[2]]];
uC[s_, 3] = 
        Chop[uC[s] /. EigenValues[[3]]];

Table[uC[s, n], {n, 3}]

The wave functions to the left of the well are

uL[s_,1]=Chop[uL[s]/.EigenValues[[1]]];
uL[s_,2]=Chop[uL[s]/.EigenValues[[2]]];
uL[s_,3]=Chop[uL[s]/.EigenValues[[3]]];

Table[uL[s, n], {n, 3}]

The wave functions to the right are decaying.

uR[s_,1]=Chop[uR[s]/.EigenValues[[1]]];
uR[s_,2]=Chop[uR[s]/.EigenValues[[2]]];
uR[s_,3]=Chop[uR[s]/.EigenValues[[3]]];

Table[uR[s, n], {n, 3}]

Now normalize the wave functions.

Clear[norm]; 
    norm[n_] := norm[
          n] = (∫_(-∞)^(-0.5))uL[s, 
                    n]^2 \[DifferentialD]s + (∫_(-0.5)^0.5) 
                    uC[s, n]^2 \[DifferentialD]s + (∫_0.5^∞) uR[s, n]^2 \[DifferentialD]s

{norm[1], norm[2], norm[3]}

A1Values = (Solve[norm[#1] == 1, 
                (Subscript[A, 1])][[2]] &) /@ {1, 
          2, 3}

Now put this altogether for the lowest state, that has the largest binding energy.

BindingEnergies[[3]]

u[s_ /; s < -(1/2), 3] = 
      uL[s, 3] /. A1Values[[3]]; 
u[s_ /; -(1/2) ≤ s ≤ (1/2), 3] = 
      uC[s, 3] /. A1Values[[3]]; 
u[s_ /; s > (1/2), 3] = 
      uR[s, 3] /. A1Values[[3]];

Plot[u[s, 3], {s, -2, 2}, Frame -> True, 
        FrameLabel -> {"s", None}, 
        PlotLabel -> "E = -34.3711"];

Here is the next state.

BindingEnergies[[2]]

u[s_ /; s < -(1/2), 2] = 
      uL[s, 2] /. A1Values[[2]]; 
u[s_ /; -(1/2) ≤ s ≤ (1/2), 2] = 
      uC[s, 2] /. A1Values[[2]]; 
u[s_ /; s > (1/2), 2] = 
      uR[s, 2] /. A1Values[[2]];

Plot[u[s, 2], {s, -2, 2}, Frame -> True, 
        FrameLabel -> {"s", None}, 
        PlotLabel -> "E = -18.4968"];

This level has the smallest binding energy

BindingEnergies[[1]]

u[s_ /; s < -(1/2), 1] = 
        uL[s, 1] /. A1Values[[1]];
u[s_ /; -(1/2) ≤ s ≤ (1/2), 1] = 
        uC[s, 1] /. A1Values[[1]];
u[s_ /; s > (1/2), 1] = 
        uR[s, 1] /. A1Values[[1]];

Plot[u[s, 1], {s, -2, 2}, Frame -> True, 
        FrameLabel -> {"s", None}, 
        PlotLabel -> "E = -0.0160771"];

Wood-Saxon Well

Clear["Global`*"];

Setting up the well

In this section we explore the energy levels of a rounded well, but we choose parameters that cause it to approximate the finite square well studied above.  We use the form of the well that best fits the energy levels of nuclei, the Woods-Saxon potential. We choose the width of the well to be 1, so the well edges are approximately at -0.5 and +0.5.

v[s_]= - v0/(1 + Exp[(Abs[s]-0.5)/0.05]);

v0=40;

Plot[v[s], {s, -2. , 2}, 
        PlotRange -> {{-2, 2}, {-40, 10}}, 
        AxesLabel -> {"s", "Wood-Saxon Well"}];

Solving the Schrodinger equation

The Schrodinger equation can be written as two coupled first order differential equations and we must solve these.

u'[s]== w[s]
w'[s]== -(e-v[s])u[s]

We shall match the wave function to an exponential function at s=-1, so w[-1] = b u[-1]

eqs[s_, b_]={u'[s]==w[s],w'[s]==((b^2) + v[s]) u[s],u[-1] == 0.1,w[-1] == b u[-1]};

Clear[soln]; 
soln[
        b_] := NDSolve[
          eqs[s, b], {u, w}, {s, -1, 
            1}][[1]]

PlotSol[b_] := 
      Plot[Evaluate[u[s] /. soln[b]], {s, -1, 1}, 
        AxesLabel -> {"s", ""}, 
        PlotLabel -> "b = " <> ToString[b]]

Search for energy levels by plotting solutions and looking for alternations in the sign of the asymptotic wave function (to the right).

Do[PlotSol[b],{b,0.1,6.5,1}]

There are three levels with b values in the intervals {0.1, 1.1}, {4.1, 5.1}, and {5.1, 6.1}.

Searching for the energy levels

We note that the growing exponential has zero amplitude when e = -(b^2) corresponds to a bound state.  So we define a function that must vanish when the coefficient of Exp[b s} is zero asymptotically at positive s values. The condition is that

u'[s] - b u[s] == w[x] - b u[x] == 0

Clear[BC];
BC[b_]:=
((b u[1] + w[1])/.Flatten[soln[b]])

Now plot this function from zero binding (b=0) to ((b^2)=100).

Plot[BC[b], {b, 0. , (Sqrt[v0])}, 
        PlotRange -> {Automatic, Automatic}];

Plot[BC[b], {b, 0, 0.2}];

bValues = {0.0618, 4.09, 5.78};

?FindRoot

EigenValues = 
      Map[FindRoot[BC[b] == 0, {b, {#, # + 0.05}}] &, 
        bValues]

Energies = -(b^2) /. EigenValues

There is a small effect of increasing the radius for the lowest state. For, the least bound state there may be some problem with the integration routine in the region where the wave function is a superposition of a growing and a decaying exponential, although the results agree to the significant figures given.

Compare the binding energies for the Woods-Saxon well with those for the finite square well.

Plots of the bound state wave functions

\(uSol[s_, n_, s0_] := 
      u[s] /. 
\(NDSolve[{\(u''\)[s] == \ \((b\^2 + v[s])\)\ u[s], 
                u[\(-s0\)] == 0.1, \(u'\)[\(-s0\)] == 0.1\ b} /. 
              EigenValues[[n]], 
            u, {s, \(-s0\), 
              s0}]\)[[1]]\)

PlotWaveFct[n_, s0_] := 
        Plot[Evaluate[uSol[s, n, s0]], {s, -s0, s0}, 
          PlotLabel ->  
            StringJoin["Binding Energy = ", 
              ToString[
                Energies[[
                  n]]]]];

Do[PlotWaveFct[n, 2], {n, 1, 3}]