Quantum Pendulum

The problem of quantizing the nonlinear mathematical pendulum was taken up by E U Condon (Phys Rev 31, 891). It seems of academic interest at first, and is not pursued in quantum texts. It has had some attention in American Journal of Physics articles: a visualization of the solutions in terms of Mathieu functions was given in AJP 71, 233; applications in molecular oscillations and hindered rotations due to electrostatic restoring torques were discussed in AJP 70, 525; semiclassical quantization was considered in AJP 48, 660; further pedagogical AJP references can be found in these articles.

The pendulum problem is interesting as it connects two distinct limits: for small-amplitude oscillations the well-known harmonic oscillator is obtained; for large eigenenergies the states should approach those of a free rotator.

These solutions are also of interest in the problem of periodic potential wells with wrap-around boundary conditions, such as discussed, e.g., in the QM text of Griffiths. Bloch's theorem can be seen at work in this context when replacing the angular variable of the quantum pendulum by a one-dimensional Cartesian coordinate.

The quantum pendulum Schroedinger equation for an object of rotational inertia m*L^2 is cast into dimensionless form. and is then connected with the Mathieu differential equation. Our aim is to find the allowed eigenenergies from a matrix representation in a Fourier basis. The eigenfunctions are then obtained from the Mathieu functions.

>

>	restart; Digits:=15: with(LinearAlgebra):

>

a:=Pi;

We selected a standard interval for the dynamical variable, i.e., [-Pi,Pi). The re-scaled potential energy is chosen such that in the oscillator limit the eigenenergies are given as the odd integers.

>

V0:=2;

>	V:=x->V0*(1-cos(Pi*x/a));

>	Vmax:=V(a);

We pick a single well (our program allows to create multiple regions adjacent to each other):

>

Nw:=1;

>	Vmw:=piecewise(seq(op([x>a*(2*j-1) and x<a*(2*j+1),V(x-2*a*j)]),j=-(Nw-1)/2..(Nw-1)/2));

>	PTS:=[seq((2*j-1)*a,j=-(Nw-1)/2..(Nw+1)/2)];

We want to have periodicity with 2a

>	Pmw:=plot(Vmw,x=PTS[1]..PTS[nops(PTS)],numpoints=500):

A periodic basis would allow for left- and right-travelling states with periodicity conditions at a, and -a. We would solve within a single well, but the states would be of travelling plane-wave type.

>	psi:=n->eval(exp(I*Pi*n*x/(Nw*a)));

>	int(conjugate(psi(0))*psi(2),x=PTS[1]..PTS[nops(PTS)]);

>	L:=int(conjugate(psi(1))*psi(1),x=PTS[1]..PTS[nops(PTS)]);

>	int(-diff(psi(3),x$2)*conjugate(psi(3)),x=PTS[1]..PTS[nops(PTS)])/L;

>	4*3^2*Pi^2/(L)^2;

>

We need positive and negative n, including 0.

>

>	n_max:=20; N:=2*n_max+1;

>	VME:=unapply(int(V(x)*exp(Pi*I*dn*x/(Nw*a)),x=x1..x2),x1,x2,dn);

>	#VME(PTS[1],PTS[2],21);

>	#limit(VME(PTS[1],PTS[2],s),s=0);

>	kh:=1; # for the cosine potential we put the strength into V(x).

>	H:=Matrix(N,N,shape=hermitian,datatype=complex[8]):

We should have a 1/sqrt(2*a) = 1/sqrt(L) normalization factor! We can fix this by dividing the matrix by L=n*2a, where n is the number of wells.

>	for i1 from 1 to N do: n1:=-n_max+i1-1: for i2 from 1 to i1 do: n2:=-n_max+i2-1: if i2=i1 or n1-n2=Nw then

>	H[i1,i2]:=evalf(1/L*kh*add(limit(VME(PTS[i0],PTS[i0+1],s),s=n1-n2),i0=1..Nw)): else H[i1,i2]:=evalf(1/L*kh*add(VME(PTS[i0],PTS[i0+1],n1-n2),i0=1..Nw)): fi: od: H[i1,i1]:=H[i1,i1]+evalf(4*n1^2*Pi^2/(L)^2): od:

>	Eigenvalues(H);

>	EV:=convert(%,list);

We could construct the eigenfunctions from the Eigenvectors command. The calculations with these states are tedious (slow). We prefer to connect with the Mathieu differential equation. This requires us to rescale theta=2*x.

>

q:=2*V0;

>	MDE:=a->diff(u(x),x$2)+(a+2*q*cos(2*x))*u(x)=0;

Here a is the eigenvalue scaled according to a=4*(EV-V0), and the potential strength parameter q=2*V0.

>	aM:=e->4*(e-V0);

>	sol:=dsolve(MDE(aM(EV[1])));

The MathieuC and MathieuS functions are symmetric and antisymmetric respectively, i.e., behave like sine and cosine. We invoke a physics argument which states that the eigenstates of the problem are either symmetric or antisymmetric, as the Hamiltonian is symmetric. The ground state is symmetric.

>	plot([eval(rhs(sol),{_C1=1,_C2=0}),eval(rhs(sol),{_C2=1,_C1=0})],x=-Pi/2..Pi/2,color=[red,blue]);

The antisymmetric solution is not acceptable, it is not periodic. The ground state is unique, it has a symmetric eigenfunction.

Let us plot a probability using the eigenenergy as a baseline, and superimpose it on the potential as a function of theta.

>	P0:=plot([EV[1],EV[1]+(eval(rhs(sol),{_C1=1,_C2=0,x=s/2}))^2],s=-Pi..Pi,color=[black,green]):

>	plots[display](Pmw,P0);

>	sol:=dsolve(MDE(aM(EV[2])));

We expect the first excited state is anti-symmetric.

>	plot([eval(rhs(sol),{_C1=1,_C2=0}),eval(rhs(sol),{_C2=1,_C1=0})],x=-Pi/2..Pi/2,color=[red,blue]);

Now the symmetric solution is not acceptable, it is not periodic. The first excited state is unique, it has an anti-symmetric eigenfunction.

>	P1:=plot([EV[2],EV[2]+(eval(rhs(sol),{_C1=0,_C2=3,x=s/2}))^2],s=-Pi..Pi,color=[black,magenta]):

>	plots[display](Pmw,P0,P1);

This looks quite analogous to the harmonic oscillator result except that we can see the considerable lowering of the first excited state (cf E_HO=3) as the potential turns around.

Based upon the harmonic oscillator we might expect the next level to be a rotator state as E_HO=5 is above the potential barrier. Our numerically obtained eigenvalue is below the barrier. Let us see what happens.

>	sol:=dsolve(MDE(aM(EV[3])));

We find a bound state which is symmetric. The antisymmetric function is not periodic.

>	plot([eval(rhs(sol),{_C1=1,_C2=0}),eval(rhs(sol),{_C2=1,_C1=0})],x=-Pi/2..Pi/2,color=[red,blue]);

Exercise 1 (trivial):

Complete the graph of the probability density for all bound states superimposed on the potential energy curve.

Now we get into the rotator regime (E>V0):

>	sol:=dsolve(MDE(aM(EV[4])));

It looks like both functions are periodic.

>	plot([eval(rhs(sol),{_C1=1,_C2=0}),eval(rhs(sol),{_C2=1,_C1=0})],x=-Pi/2..Pi/2,color=[red,blue]);

>

>	fS:=eval(rhs(sol),{_C1=1,_C2=0});

>	fA:=eval(rhs(sol),{_C2=1,_C1=0});

>	evalf(limit(fS,x=-Pi/2));

>	evalf(limit(fS,x=Pi/2));

>	evalf(limit(diff(fS,x),x=-Pi/2));

>	evalf(limit(diff(fS,x),x=Pi/2));

The symmetric function is not really periodic!

>	evalf(limit(fA,x=-Pi/2));

>	evalf(limit(fA,x=Pi/2));

>	evalf(limit(diff(fA,x),x=-Pi/2));

>	evalf(limit(diff(fA,x),x=Pi/2));

The anti-symmetric function is perfectly acceptable! Now look at the next eigenvalue:

>	sol:=dsolve(MDE(aM(EV[5])));

Again, it looks like both functions are periodic. However, this is probably not true. The antisymmetric function is not matching at the boundaries, only almost so.

>	plot([eval(rhs(sol),{_C1=1,_C2=0}),eval(rhs(sol),{_C2=1,_C1=0})],x=-Pi/2..Pi/2,color=[red,blue]);

>	fS:=eval(rhs(sol),{_C1=1,_C2=0});

>	fA:=eval(rhs(sol),{_C2=1,_C1=0});

>	evalf(limit(fS,x=-Pi/2));

>	evalf(limit(fS,x=Pi/2));

>	evalf(limit(fA,x=-Pi/2));

>	evalf(limit(fA,x=Pi/2));

What is our finding?

In the regime where we expect rotator solutions we are getting closely spaced anti-symmetric and symmetric solutions. These look like bound states. We cannot form linear combinations that would correspond to travelling waves as the energies are not degenerate. Rotation is inhibited, because in quantum scattering both reflection and transmission are probable for energies not too high above the potential barrier.

Exercise 2:

Explore what happens for the higher energies. Do true rotator states appear? If the symmetric and anti-symmetric solutions appeat at identical energies, then one can form linear combinations along the lines of exp(I*s) = cos(s) + I*sin(s). Both signs are possible resulting in counterclockwise and clockwise rotator states at the same energy.

Exercise 3:

What is the energy spectrum of the pure rotator. How close are the higher states to the expected m-dependence, where m counts the states.

>

>

>

>