Harmonic and anharmonic motion in the microscopic, quantum world

We will develop some intuition for the description of motion in the microscopic, i.e., quantum domain. We start with the harmonic oscillator, for which we develop solutions to Schroedinger's equation in analytic form, and then we move on to solve the problem numerically for an anharmonic oscillator.

Without further motivation we state that the allowed energies for a quantum particle bound by a potential V ( x ) is given by a differential-equation eigenvalue problem:

> V:=k/2*x^2;

> SE:=-h_^2/(2*m)*diff(psi(x),x$2)+V*psi(x)=lambda*psi(x);

Let us make life simple by adopting Bohr units (h_=m_e=e=1), and by choosing k=1:

> h_:=1: m:=1: k:=1:

> SE;

The physical constraints on this problem are that the wavefunction psi has to vanish at infinity, i.e., for large positive and negative values of x .

We will be interested only in such solutions, i.e., we will discard solutions that do not vanish asymptotically.

On the other hand one can show that the solutions have to be either symmetric, or antisymmetric functions, i.e.,

either: psi(x) = psi(-x), or: psi(x) = - psi(-x)

for the acceptable solutions. These conditions translate into conditions at x = 0 as follows:

either: psi( x ) = const, D(psi)(0)=0 ; or: psi(0)=0, D(psi)(0) = const.

The two undetermined constants fix the amplitude of the solution, and we recognize that this amplitude is arbitrary from a mathematical point of view. This is the consequence of solving a homogeneous problem. Let us start with the symmetric solutions:

> ICs:=psi(0)=1,D(psi)(0)=0;

> sol:=dsolve({SE,ICs},psi(x));

So, it's not as simple as that in Maple... Let us substitute a number for lambda (we are using special units, so the lambda - value had dimensions of energy.

Let us try whether the energy can be zero, i.e., whether we can have a particle sitting at the bottom of the well, which corresponds to having the mass simply at rest in the equilibrium position of the spring:

> sol:=dsolve({subs(lambda=0,SE),ICs},psi(x));

> plot(rhs(sol),x=-5..5,view=[-5..5,0..5]);

> evalf(subs(x=-1,rhs(sol)));

We don't get an acceptable answer for negative x , but presumably we could use the positive- x results, and use the symmetry to extend the solution to negative x .

The reason for not plotting at negative x is that the solution has a sqrt( x ) factor, which is imaginary for negative x values.

The main point is that the solution is not acceptable for large x . It blows up. The question is: Can we find energy values such that the solution is normalizable?

Strategy: try lambda=1/4, and lambda=1:

> sol:=dsolve({subs(lambda=1/4,SE),ICs},psi(x));

> plot(rhs(sol),x=-5..5,view=[-5..5,0..5]);

> sol:=dsolve({subs(lambda=3/4,SE),ICs},psi(x));

> plot((rhs(sol)),x=-5..5,view=[-5..5,-5..5],numpoints=500);

Somewhere in-between there should be a solution that approaches the x -axis asymptotically.

> sol:=dsolve({subs(lambda=1/2,SE),ICs},psi(x));

> plot(rhs(sol),x=-5..5,view=[-5..5,0..2]);

To complete the interpretation of the Schroedinger wave equation result we compute the probability density for this state. The square of the wave function is providing the probability for finding a particle at location x . For this interpretation to hold we need to normalize the solution of the homogeneous Schroedinger eigenvalue problem in the following way:

> int(rhs(sol)^2,x=-infinity..infinity);

> WF:=1/sqrt(%)*rhs(sol);

> int(WF^2,x=-infinity..infinity);

Now we graph the probability in such a way as to show that there is a finite probability to find the quantum particle outside the classically allowed range. For this we plot the density superimposed on the potential while using the energy value as a baseline:

> plot([V,1/2,WF^2+1/2],x=-2..2,color=[red,blue,green]);

We see that the classical turning point for a particle of energy 1/2 is at x = 1, and x = -1. For | x | > 1 a classical particle is not to be found with this energy value, since the kinetic energy is positive (by definition), and the potential energy exceeds the value of 1/2 in this range. Yet the quantum particle (which does not follow any trajectory) has a finite probability to be found in the classically forbidden regime. This regime is called the quantum tunneling region. It is a feature of Schroedinger solutions of all finite potentials that quantum particles penetrate beyond the classical turning points for the given amount of energy.

Exercise 1:

Search for the next two acceptable valuea of lambda. Hint: use the asymptotic behaviour of the solution to determine whether an energy (lambda) value can be found such that the solution merges with the x -axis asymptotically. Tabulate the lambda-values together with the one found above.

Graph the probability distributions for the properly normalized wavefunctions.

We proceed with the determination of antisymmetric solutions:

> ICa:=psi(0)=0,D(psi)(0)=1;

> sol:=dsolve({subs(lambda=1/2,SE),ICa},psi(x));

> plot(rhs(sol),x=-5..5,view=[-5..5,0..2]);

> sol:=dsolve({subs(lambda=5/2,SE),ICa},psi(x));

> plot(rhs(sol),x=-5..5,view=[-5..5,-2..2]);

> sol:=dsolve({subs(lambda=3/2,SE),ICa},psi(x));

> plot(rhs(sol),x=-5..5,view=[-5..5,-2..2]);

We carry out the normalization of the state, and display the probability distribution:

> int(rhs(sol)^2,x=-infinity..infinity);

> WF:=1/sqrt(%)*rhs(sol);

> int(WF^2,x=-infinity..infinity);

Now we graph the probability in such a way as to show that there is a finite probability to find the quantum particle outside the classically allowed range. For this we plot the density superimposed on the potential while using the energy value as a baseline:

> plot([V,3/2,WF^2+3/2],x=-2.5..2.5,color=[red,blue,green]);

We observe again the fact that the quantum particle penetrates into the classically forbidden region. We also observe that the particle has zero probability to be found at the origin ( x =0) when it is in an antisymmetric state. Given that there are no continuous trajectories for quantum particles this is not a problem, however: the interpretation of the quantum probability distribution is the following:

If repeated measurements are made on an ensemble of particles prepared in a given eigenstate (the ensemble can be simultaneous, or there could be one particle at a time with repeated measurements), then each time the probability is recorded (a histogram is made of the position for the ensemble), the distribution will be an approximation to the above curve. There is no meaning to to question as to how an individual quantum particle moved in such a way as to cross x =0.

Exercise 2:

Find two more anti-symmetric solutions, and tabulate the eigenvalues.

Combine these eigenvalues with the ones found for the symmetric case, and observe how they behave.

Graph the probability distributions using the properly normalized wavefunction.

Exercise 3:

Change the value of the spring constant k . Pick k =2 and k =4 respectively, and determine the lowest 6 eigenvalues in each case. Draw your conclusions about the energy spectrum.

Exercise 4:

Keep the value of the spring constant k =1, but change the mass to m =2 (in Bohr units). Determine the 6 lowest energy values, and draw your conclusions about the energy spectrum.

Let us go back to the general case, and show what the exact eigenfunctions look like:

> V:=1/2*m*omega^2*x^2;

> SE:=-h_^2/(2*m)*diff(psi(x),x$2)+V*psi(x)=lambda*psi(x);

> alpha:=sqrt(omega);

> with(orthopoly);

> phi:=n->exp(-1/2*omega*x^2)*H(n,alpha*x);

> subs(psi(x)=phi(0),SE);

> simplify(%);

> solve(%,lambda);

Clearly, the ground-state energy equals 1/2 omega. Let's continue:

> subs(psi(x)=phi(1),SE);

> simplify(%);

> solve(%,lambda);

For the first excited state we find that the energy equals 3/2 omega. Substitute different values of n and verify the eigenvalue spectrum.

The only remaining point concerns the normalization of the eigenfunctions. The interpretation of the functions phi( n , x ) is that their square provides the probability to find the particle at location x . Therefore the sum of phi( x , n )^2 over all locations x , i.e., the integral is chosen to equal one. The functions as defined do not meet that requirement, and one has to find an omega- and n -dependent factor to achieve this goal. For example:

> assume(omega>0);

> int(phi(5)^2,x=-infinity..infinity);

The normalization factor can be found to be:

>

For other potentials (e.g., an x^4 anharmonic oscillator) the eigenvalues and eigenfunctions cannot be found in terms of simple functions. Nevertheless, we can apply numerical techniques to find approximate eigenvalues and eigenfunctions using dsolve[numeric]. We need to apply an iteration in the energy eigenvalue:

> restart;

> V:=x^4;

> SE:=-1/2*diff(psi(x),x$2)+V*psi(x)=E_t*psi(x);

> unprotect(Psi);

> ICs:=psi(0)=1,D(psi)(0)=0;

> E_t:=0.668;

> sol:=dsolve({SE,ICs},psi(x),numeric,output=listprocedure):

> Psi:=subs(sol,psi(x)):

> plot('Psi(x)',x=0..4,-0.5..1.5);

The energy trial value was determined by trial and error. An energy value (and the corresponding eigenfunction) is considered acceptable when the function is close to the x -axis for some range. This range can be extended by tuning the energy value to higher accuracy.

Exercise 5:

Find approximate eigenenergies for two more eigenvalues (above the ground state). Observe the nodal structure of the wavefunctions and compare with the exact harmonic oscillator result above. Calculate also several anti-symmetric eigenfunctions by changing the initial conditions for the integration appropriately.

Compare the spacing between these eigenvalues with the spacing property found for the harmonic oscillator.

Exercise 6:

Solve numerically for the lowest six eigenvalues and eigenfunctions of the potential V ( x )=| x |. Discuss the properties of the spacing of the eigenvalues as compared to the harmonic oscillator.

Note that a systematic method to improve the trial values for the eigenenergies involves a bisection: one can systematically cut the energy range defined by two brackets (with a solution going to positive infinity at large x for one bracketing value, and to negative infinity for the other) in half, and determine which half contains the asymptotically acceptable solution to reduce the bracket systematically by factors of two.

>

Here we do the exercise for the linear potential, and write a bisection algorithm.

> restart;

> V:=x;

> SE:=-1/2*diff(psi(x),x$2)+V*psi(x)=E_t*psi(x);

> unprotect(Psi);

> ICs:=psi(0)=1,D(psi)(0)=0;

By trial and error we find that the lowest eigenvalue is bracketed by 0.8 and 0.9:

> E_t:=0.825;

> sol:=dsolve({SE,ICs},psi(x),numeric,output=listprocedure):

> Psi:=subs(sol,psi(x)):

> plot('Psi(x)',x=0..4,-0.5..1.5);

We choose a sufficiently large value of x for which we want the solution to vanish:

> X:=4;

> SEroot:=proc(lambda) local SE,sol,Psi; global ICs,X,V;

> SE:=-1/2*diff(psi(x),x$2)+V*psi(x)=lambda*psi(x);

> sol:=dsolve({SE,ICs},psi(x),numeric,output=listprocedure): unprotect(Psi):

> Psi:=subs(sol,psi(x)):

> Psi(X); end:

>

> fsolve('SEroot(E_t)',E_t,0.8..0.9);

Error, (in fsolve) invalid arguments

Unfortunately, it is not as easy as that...

We first write a basic bisection step: given a map f and a bracketing interval, we come up with a new bracket for the root which is half the size:

> Bisect:=proc(f,xL) local x1,x2,xh,f1,f2,fh; x1:=xL[1]: x2:=xL[2]: f1:=f(x1): f2:=f(x2): if f1*f2>0 then RETURN("Same sign on both ends of the bracket, cannot proceed") fi: xh:=0.5*(x1+x2): fh:=f(xh): if f1*fh<0 then [x1,xh] else [xh,x2] fi: end:

> Bisect(SEroot,[0.8,0.9]);

Now we need a driver, to drive the trial energy bracket to some desired tolerance.

In this procedure we limit ourselves to a maximum of 100 iterations. We define an internal variable rL to represent the bracket as a list. We cannot overwrite the list that was passed down originally, i.e.. xL (Maple syntax). When successful we return the midpoint of the bracket.

We also catch the incidence when Bisect returns with an error message and the result of rL is not of type list .

> SEsolve:=proc(xL,tol) local i,rL; rL:=xL: for i from 1 to 100 do: rL:=Bisect(SEroot,rL): if not type(rL,list) then RETURN(rL) fi: if abs(rL[2]-rL[1]) < evalf(tol) then RETURN(0.5*(rL[1]+rL[2])) else fi: od: RETURN("Did not converge in 100 iterations in SEsolve"); end:

> SEsolve([0.8,0.9],0.001);

We should keep in mind that the number is only accurate to the level determined by the variable tol :

> SEsolve([0.8,0.9],0.00001);

> SEsolve([0.8,0.9],0.0000001);

> E_t:=%;

> sol:=dsolve({SE,ICs},psi(x),numeric,output=listprocedure):

> Psi:=subs(sol,psi(x)):

> plot('Psi(x)',x=0..4,-0.5..1.5);

Note what can happen when the original bracket is too wide:

> SEsolve([0.1,1.0],0.001);

> SEsolve([0.1,2.0],0.001);

> SEsolve([0.1,3.0],0.001);

> SEsolve([0.1,4.0],0.001);

> SEsolve([0.1,5.0],0.001);

> SEsolve([4.,5.0],0.001);

> SEsolve([0.1,6.0],0.001);

The program still requires some attention, but the process of finding the eigenvalues is automated to some extent.

Note that in the present form the program calculates the symmetric eigenfunctions (and the corresponding eigenvalues) only.

How far does the particle in the ground state explore the classically forbidden region in this case?

> E0:=SEsolve([0.1,2.0],0.00001);

> E_t:=E0;

> sol:=dsolve({SE,ICs},psi(x),numeric,output=listprocedure):

> Psi:=subs(sol,psi(x)):

> P1:=plot(E0+'Psi(x)^2',x=0..3,0..2,color=green):

> P2:=plot([x,E0],x=0..3,color=[red,blue]):

> plots[display](P1,P2);

Note that the green curve does not show a properly normalized probability distribution. It is simply based on a solution that satisfies psi(0)=1.

Solved Problem:

Determine several low-lying eigenvalues and eigenfunctions for a 1d potential with a Rydberg series, namely -1/sqrt(1+x^2)

> V:=-1/sqrt(1+x^2);

> X:=10;

> E0:=SEsolve([-1.,-0.2],0.0001);

> E_t:=E0;

> SE:=-1/2*diff(psi(x),x$2)+V*psi(x)=E_t*psi(x);

> sol:=dsolve({SE,ICs},psi(x),numeric,output=listprocedure):

> Psi:=subs(sol,psi(x)):

> P1:=plot(E0+'Psi(x)^2',x=0..5,-1..1,color=green):

> P2:=plot([V,E0],x=0..5,color=[red,blue]):

> plots[display](P1,P2);

> evalf(Int('Psi(x)^2',x=0..5),5);

Error, (in evalf/int) cannot evaluate boolean: t-abs(-.156399739200000000+t) < 0

We can estimate the normalization integral:

> dx:=0.05: N:=5/dx: evalf(2*dx*add(Psi(i*dx)^2,i=0..N));

To obtain the actual probability to find a particle at a given location x one would need to divide the square of the wavefunction by that number.

Let us calculate the first excited state by defining a procedure with a choice of appropriate boundary condition:

> SEroot:=proc(lambda) local SE,IC,sol,Psi; global X,V,f_sym;

> SE:=-1/2*diff(psi(x),x$2)+V*psi(x)=lambda*psi(x);

> if f_sym=S then IC:=psi(0)=1,D(psi)(0)=0; elif f_sym=A then IC:=psi(0)=0,D(psi)(0)=1; fi;

> sol:=dsolve({SE,IC},psi(x),numeric,output=listprocedure): unprotect(Psi):

> Psi:=subs(sol,psi(x)):

> Psi(X); end:

> f_sym:=S;

> E0:=SEsolve([-1.,-0.2],0.0001);

> f_sym:=A;

> E1:=SEsolve([-1.,-0.2],0.0001);

> ICa:=psi(0)=0,D(psi)(0)=1;

> E_t:=E1;

> SE:=-1/2*diff(psi(x),x$2)+V*psi(x)=E_t*psi(x);

> sol:=dsolve({SE,ICa},psi(x),numeric,output=listprocedure):

> Psi:=subs(sol,psi(x)):

> P1:=plot(E1+'Psi(x)^2',x=0..8,-1..1,color=green):

> P2:=plot([V,E1],x=0..8,color=[red,blue]):

> plots[display](P1,P2);

Exercise 7:

Calculate the next two energye levels (one symmetric, and one antisymmetric eigenfunction). Allow for a sufficiently large value of X such that the solution is not constrained by forcing the zero value of the wavefunction at too small an x value. X has to be chosen large enough such that the final answer is independent of it.

Note that the above potential is one-dimensional, but shares some feature of the hydrogen atom problem:

there are bound and unbound solutions; the bound solutions are characterized by negative energies; the unbound E >0 solutions are not quantized.

there is an infinite series of bound states (Rydberg series)

>