Excercise sheet 9 -- solutions

Schrödinger equation

1. Try to find the first excited state of the hydrogen atom and the corresponding wavefunction.

2. When solving the eigenvalue E we do not necessarily need to write a loop ourselves. We are only looking for a root of the function Δ(E), so, if we can write Δ as a Mathematica function, we can use the FindRoot command to locate the energies E. Try to do this. You can use the Module command to wrap the commands for calculating Δ into a function. Module command takes two arguments, first a list of local variables (that is, variables that you only use inside the function) and second, the commands separated by semicolons (;). Inside the Module command, you can use the Return -command to return the value for delta. Your function Δ in Mathematica could look something like this:

Delta[e_] := Module[
	{sol0, sol1, …},
	sol0=NDSolve[{eq, …];
	…
	diff = ...
	Return[diff]];
	
You can now plot the function Delta, and use FindRoot to get the energy eigenvalues!

Solutions

1

So, we are trying to solve the equation
	-y'' (r) - (2/r) y(r) = 
        E y(r),
or, more specifically, the first excited state of this.

First, a few words on the excited states. The discrete states of a quantum mechanical system can be characterized by so-called quantum numbers, which then serve to determine the energy, momentum, angular momentum etc. of the system. Suppose a system is completely determined by its energy, or its energy quantum number. This applies well to the states in question here, i.e. the lowest energy states of the hydrogen atom. The energy levels of the system, in this case the electron, are of the form
f(k), where k is the energy quantum number. The range of possible values of k determines the range of the energy values.

Usually, in numerical calculations, we don't know the function f. From the outset of the problem, from the zero of the Coulomb potential, we can however deduce that E < 0, so as k increases, the energy approaches zero from below. The ground state energy was found to be -1, so all (an infinite number of) excited states are in the interval -1 < E < 0. There is, a priori, no way of telling where the 'second' state is. The only way of finding out is to calculate the values of the function Δ(E) in sufficiently small steps, starting from -1 and see where its next zero lies.

What is the significance of the function Δ, defined as
	\(TraditionalForm\`Δ(
          E) = \(y\_0' \((r\_inter)\)\)\/\(\(y\_0\)(r\_inter)\) - \(y\_∞' \((r\_inter)\)\)\/\(\(y\_∞\)(r\_inter)\)\),
where (Subscript[y, 0])(r) is the solution in the interval (Subscript[r, 0]) < r < (Subscript[r, inter]) and (Subscript[y, ∞])(r) gives the solution for (Subscript[r, inter]) < r < (Subscript[r, 1])? Here (Subscript[r, 0]) is the point corresponding to the origin in our numerical calculations, and (Subscript[r, 1]) corresponds to infinity (actual zero and infinity are hard to handle in numerics). The solution of the above equation starts from the formulation of the boundary conditions, in this case they are y((Subscript[r, 0])) = (Subscript[C, 1]) (Subscript[r, 0]) and y((Subscript[r, 1])) = (Subscript[C, 2]) (e^(-(Sqrt[E]) (Subscript[r, 1]))), where (Subscript[C, 1]) and (Subscript[C, 2]) are some constants, which are determined by matching of the partial solutions and the normalization (for the beginning part of the solution, these constants can be set to unity). Corresponding to these boundary conditions, the differential equation has solutions only for a certain set of values of E. A possible solution y of the equation, corresponding to some (guessed) value of E, has to satisfy the conditions that y and y' are continuous, i.e. (Subscript[y, 0])((Subscript[r, inter]))=(Subscript[y, ∞])((Subscript[r, inter])) and (Subscript[y, 0])' ((Subscript[r, inter])) = 
        (Subscript[y, ∞])' ((Subscript[r, inter])). Dividing these equations yields the condition that Δ(E) = 0. What this all means is that if Δ(E) ≠ 0, there are no solutions corresponding to this value of E.

Alright, let's get on with the solution. First, we'll input some definitions. The crucial ones here are the values of r demarcating the different regions. With an inappropriate choice of (Subscript[r, 1]), the function Δ might behave in a way making it difficult to find the zeros. Again, there is no way of telling what is an appropriate choice, but as a rule, you should make (Subscript[r, 1]) large enough and (Subscript[r, inter]) not too large. In the second excercise, a better way of finding the zeros of Δ is presented. Note also that the constants (Subscript[C, 1]) and (Subscript[C, 2])have been set equal to one in the boundary condition functions y0 and y1.

Clear[e, eq];
eq := -y''[r] - (2/r) y[r] == 
          e y[r];
r0 = (10.0^(-8)); r1 = 50.0; rinter = 1;
y0[r_] = r; dy0[r_] = D[y0[r], r];
y1[r_] = Exp[-(Sqrt[-e]) r]; dy1[r_] = D[y1[r], r];

Then, give some off-the-hat values to E and calculate the corresponding value of Δ. You can do this simply by hand, setting one value after the other, or devise some neat automatic method. For example, you could give E the values -0.9, -0.8,  ...  , -0.2 and discover that Δ(E) changes sign in the interval -0.3 < E < -0.2.

In the following, diff1=Δ(e1) and diff2=Δ(e2).

e1 = -0.30; e2 = -0.20;
e = e1;
sol0 = 
        Part[NDSolve[{eq, 
              y[r0] == y0[r0], y'[r0] == dy0[r0]}, 
            y[r], {r, r0, rinter}, AccuracyGoal -> 15], 
          1];
sol1 = 
        Part[NDSolve[{eq, 
              y[r1] == y1[r1], y'[r1] == dy1[r1]}, 
            y[r], {r, rinter, r1}, AccuracyGoal -> 15], 
          1];
u0[r_] = y[r] /. sol0; 
    u1[r_] = y[r] /. sol1;
d0[r_] = D[u0[r], r]/u0[r]; 
    d1[r_] = D[u1[r], r]/u1[r];
Print["diff1 = ", 
      diff1 = d0[rinter] - d1[rinter]]
e = e2;
sol0 = 
        Part[NDSolve[{eq, 
              y[r0] == y0[r0], y'[r0] == dy0[r0]}, 
            y[r], {r, r0, rinter}, AccuracyGoal -> 15], 
          1];
sol1 = 
        Part[NDSolve[{eq, 
              y[r1] == y1[r1], y'[r1] == dy1[r1]}, 
            y[r], {r, rinter, r1}, AccuracyGoal -> 15], 
          1];
u0[r_] = y[r] /. sol0; 
    u1[r_] = y[r] /. sol1;
d0[r_] = D[u0[r], r]/u0[r]; 
    d1[r_] = D[u1[r], r]/u1[r];
Print["diff2 = ", diff2 = d0[rinter] - d1[rinter]]

Now, since diff1 and diff2 have opposite signs, it has to be the case that Δ(E) = 0 for some value of E for which e1 < E < e2. We will search a more accurate approximation of the eigenvalue by dividing the interval [e1, e2] in half and calculate Δ(((e1 + e2)/2)), from the sign of which we can then deduce the location of the eigenvalue with double precision, i.e. it is either in the interval [e1, ((e1 + e2)/2)] or [((e1 + e2)/2), e2]. This process of dividing the interval in half is then repeated until some predefined, desired accuracy is reached.

To decide which of the subintervals to examine next, we calculate the product diff*diff1, where diff stands for Δ(((e1 + e2)/2)). If this product is negative, it means that diff and diff1 have different signs and the eigenvalue must lie between e1 and ((e1 + e2)/2). If it is positive, the eigenvalue is in the other interval.

acc = (10^(-12));
\(While[
      Abs[e1 - e2] > 
        acc, 
e\  = \ \((e1 + e2)\)/
          2; 
sol0 = 
        Part[NDSolve[{eq, 
              y[r0] == y0[r0], \(y'\)[r0] == dy0[r0]}, 
            y[r], {r, r0, rinter}, AccuracyGoal -> 15], 
          1]; 
sol1 = 
        Part[NDSolve[{eq, 
              y[r1] == y1[r1], \(y'\)[r1] == dy1[r1]}, 
            y[r], {r, rinter, r1}, AccuracyGoal -> 15], 
          1]; 
u0[r_] = y[r] /. sol0; 
      u1[r_] = y[r] /. sol1; 
d0[r_] = D[u0[r], r]/u0[r]; 
      d1[r_] = D[u1[r], r]/u1[r]; 
diff = 
        d0[rinter] - d1[rinter]; 
If[diff*diff1\  >= \ 0, 
        e1 = e; diff1 = diff, e2 = e; diff2 = diff;]
]\)
Print["eigenvalue = ", e]
Print["difference of the logarithmic derivatives = ", 
      diff]

This solution can be verified by comparison with the analytical solution, (Subscript[E, n]) = ((-1)/(n^2)), i.e. (Subscript[E, 2]) = -(1/4).

Then we can construct the wave function. A preliminary solution has been calculated above: y(r)={({{(Subscript[u, 0])(r),  when (Subscript[r, 0]) < r < (Subscript[r, inter])}, {(Subscript[u, 1])(r),  when (Subscript[r, inter]) < r < (Subscript[r, 1])}}). This is however not yet satisfactory. The matching is not perfect, because the condition Δ(E) = 0 allows any possible constants to be cancelled out. We need two constants, one for the matching and one for normalization of the wave function.

fig0 = 
        Plot[u0[r], {r, r0, rinter}, 
          DisplayFunction -> Identity];
fig1 = 
        Plot[u1[r], {r, rinter, 15}, 
          DisplayFunction -> Identity];
Show[fig0, fig1, 
        DisplayFunction -> $DisplayFunction];

We'll start the search of the constant from the normalization condition. We'll denote y(r)={({{(Subscript[C, 0]) (Subscript[u, 0])(r),  
                  when (Subscript[r, 0]) ≤ r ≤ (Subscript[r, inter])}, {(Subscript[C, 1]) (Subscript[u, 1])(r),  
                  when (Subscript[r, inter]) < r ≤ (Subscript[r, 1])}}). Then the normalization condition becomes
	\(TraditionalForm\`1 = \(∫\+V\(\(|\)\(ψ(
                r)\)\( | \^2\)\(dV\)\) = \(∫\+Ω 
                dΩ \(\(∫\+0\)\&∞\( r\^2\) 
                    dr\ \(ψ(r)\)\^2\) = 
            4  π \(\(∫\+0\)\&∞ r\^2\ \(\(y(r)\)\^2\/r\^2\) dr\)\)\)\)
	\(TraditionalForm\`1\/\(4  π\) = \(\(∫\+\(r\_0\)\)\&\(r\_1\)\((y(r))\)\^2\ dr = \(\(∫\+\(r\_0\)\)\&\(r\_inter\)\(\((\(C\_0\) \(\(u\_0\)(r)\))\)\^2\) 
                  dr + \(∫\+\(r\_inter\)\)\&\(r\_1\)\(\((\(C\_1\) \(\(u\_1\)(r)\))\)\^2\) 
                  dr = \(C\_0\^2\)(\(∫\+\(r\_0\)\)\&\(r\_inter\)\(\((\(u\_0\)(r))\)\^2\) 
                    dr + \(∫\+\(r\_inter\)\)\&\(r\_1\)\(\((\(C\_1\/C\_0\) \(\(u\_1\)(r)\))\)\^2\) dr)\)\)\).
Then, denote (Subscript[C, 2]) = ((Subscript[C, 1])/(Subscript[C, 0])), and the matching condition (Subscript[C, 0]) (Subscript[u, 0])(
            (Subscript[r, inter])) = (Subscript[C, 1]) (Subscript[u, 1])((Subscript[r, inter])) becomes \(TraditionalForm\`C\_2 = \(\(u\_0\)(r\_inter)\)\/\(\(u\_1\)(r\_inter)\)\).

C2 = u0[rinter]/u1[rinter];
fig0 = 
        Plot[u0[r], {r, r0, rinter}, 
          DisplayFunction -> Identity];
fig1 = 
        Plot[C2 u1[r], {r, rinter, 15}, 
          DisplayFunction -> Identity];
Show[fig0, fig1, 
        DisplayFunction -> $DisplayFunction];

Then we calculate the normalization coefficient. Note that the actual wave function is the function y(r)/r. That is plotted in the following.

nn = 
      NIntegrate[u0[r]^2, {r, r0, rinter}] + 
        NIntegrate[(((C2 u1[r]))^2), {r, rinter, r1}]; 
    C0 = 1/Sqrt[4  π nn];
fig0 = 
        Plot[C0 u0[r]/r, {r, r0, rinter}, 
          DisplayFunction -> Identity];
fig1 = 
        Plot[C0 C2 u1[r]/r, {r, rinter, 15}, 
          DisplayFunction -> Identity];
Show[fig0, fig1, 
        DisplayFunction -> $DisplayFunction];

This can then be compared to the analytical solution. In Wikipedia, there is some visualization of the orbitals, and you can do a comparison if you wish. http://en.wikipedia.org/wiki/Hydrogen_atom#Visualizing_the_hydrogen_electron_orbitals

2

Here is the implementation of the Δ function, illustrating the handyness of the Module command.

Clear[Delta, e];
Delta[e_] := 
      Module[{}, 
Clear[y, 
            r]; 
yhtalo = y''[r] - (-2/r  - e) 
                  y[r] == 0; 
c1 = 
            1; 
c2 = 1; 
yal[r_] = 
            c1*r; 
dyal[r_] = 
            D[yal[r], r]; 
ylo[r_] = 
            c2*E^(-Sqrt[-e]*r); 
dylo[r_] = 
            D[ylo[r], r]; 
tark = 
            10^(-12) ; 
ralku =  
            10^(-8) ; 
rloppu  = -Log[
                      eps]/(Sqrt[-e]) /. {eps -> 
                  1*^-12};  (* <- 
                  this guarantees that the tail of the function is decaying at rloppu *) 
rvali = 
            rloppu/15; 
rata = 
            Part[NDSolve[{yhtalo, 
                  y[ralku] == yal[ralku], y'[ralku] == 
                    dyal[ralku]}, y[r], {r, ralku, rvali}, 
                AccuracyGoal -> 15], 1]; 
alku[r_] = 
            y[r] /. rata; 
logalk[r_] = 
            D[alku[r], r]/alku[r]; 
ratl = 
            Part[NDSolve[{yhtalo, 
                  y[rloppu] == ylo[rloppu], y'[rloppu] == 
                    dylo[rloppu]}, y[r], {r, rloppu, rvali}, 
                AccuracyGoal -> 15], 1]; 
loppu[r_] = 
            y[r] /. ratl; 
loglop[r_] = 
            D[loppu[r], r]/loppu[r]; 
ero = 
            logalk[rvali] - 
              loglop[rvali]; 
ero
] /; 
        NumericQ[e]
 (* 
      This last command tells Mathematica NOT to evaluate the function before it has a numerical argument *) 

    
Plot[Delta[e], {e, -1.5, -0.1}, 
      PlotRange -> {-2, 2}]
FindRoot[
      Delta[e] == 
        0, {e, -1.1, -0.9}]