]> harmonic_oscil_1d_vs_3d.html harmonic_oscil_1d_vs_3d.mws

The harmonic oscillator in one and in three dimensions

Marc R. Roussel

Department of Chemistry and Biochemistry

University of Lethbridge

One dimension

I'm going to do a few calculations with the ground-state harmonic oscillator wavefunction. This wavefunction is

>    psi0_x := x -> A1*exp(-alpha*x^2/2);

psi0_x := proc (x) options operator, arrow; A1*exp(-1/2*alpha*x^2) end proc

>    assume(alpha>0);

We need to normalize the wavefunction.

>    solve(int(psi0_x(x)^2,x=-infinity..infinity)=1,A1);

1/Pi^(1/2)*(Pi^(1/2)*alpha^(1/2))^(1/2), -1/Pi^(1/2)*(Pi^(1/2)*alpha^(1/2))^(1/2)

>    A1 := %[1];

A1 := 1/Pi^(1/2)*(Pi^(1/2)*alpha^(1/2))^(1/2)

Suppose that we want to calculate the average kinetic energy. This would be

>    -hbar^2/(2*mu)*int(psi0_x(x)*diff(psi0_x(x),x$2),x=-infinity..infinity);


Similarly, the average of x ^2 is

>    avg_x2_1d := int(psi0_x(x)^2*x^2,x=-infinity..infinity);

avg_x2_1d := 1/(2*alpha)

It's all pretty simple in one dimension.

Three dimensions


>    psi0_r := r -> A3*exp(-alpha*(r-req)^2/2);

psi0_r := proc (r) options operator, arrow; A3*exp(-1/2*alpha*(r-req)^2) end proc

>    assume(req>0);

>    solve(int(psi0_r(r)^2,r=0..infinity)=1,A3);

1/(Pi^(1/2)+erf(alpha^(1/2)*req)*Pi^(1/2))*2^(1/2)*(Pi^(1/2)*(1+erf(alpha^(1/2)*req))*alpha^(1/2))^(1/2), -1/(Pi^(1/2)+erf(alpha^(1/2)*req)*Pi^(1/2))*2^(1/2)*(Pi^(1/2)*(1+erf(alpha^(1/2)*req))*alpha^(1...

>    A3:=%[1];

A3 := 1/(Pi^(1/2)+erf(alpha^(1/2)*req)*Pi^(1/2))*2^(1/2)*(Pi^(1/2)*(1+erf(alpha^(1/2)*req))*alpha^(1/2))^(1/2)

Note that the normalization factor is different  from the one-dimensional normalization factor. This isn't too surprising because in the simple one-dimensional case, we arbitrarily extend the range of integration over the entire real line. The real question is whether it makes a significant difference to our estimation of (e.g.) expectation values. Let's start with <x^2>:

>    avg_x2_3d := int(psi0_r(r)^2*(r-req)^2,r=0..infinity);

avg_x2_3d := 1/2*(alpha^(3/2)*exp(alpha*req^2)*Pi+alpha^(3/2)*exp(alpha*req^2)*Pi*erf(alpha^(1/2)*req)-2*Pi^(1/2)*req*alpha^2)/(1+erf(alpha^(1/2)*req))/alpha^(5/2)/exp(alpha*req^2)/Pi
avg_x2_3d := 1/2*(alpha^(3/2)*exp(alpha*req^2)*Pi+alpha^(3/2)*exp(alpha*req^2)*Pi*erf(alpha^(1/2)*req)-2*Pi^(1/2)*req*alpha^2)/(1+erf(alpha^(1/2)*req))/alpha^(5/2)/exp(alpha*req^2)/Pi

That's pretty complicated. The best way to decide whether or not this number is significantly different from the simple 1d result is by substituting in some values.

First, let's estimate alpha. The reduced mass for a diatomic molecule is typically about 10^(-26) kg (less for HX molecules). Force constants ( k ) are usually a few hundred N/m. For the sake of argument, take 200 N/m. Typical values for alpha are therefore of order

>    alpha_typical := sqrt(1e-26*200)/1e-34;

alpha_typical := .1414213562e23

Let's estimate < x ^2> for the simple 1-d case:

>    subs(alpha=alpha_typical,avg_x2_1d);


Typical bond lengths are of order 10^(-10) m. In 3-d we therefore have

>    subs(alpha=alpha_typical,req=1e-10,avg_x2_3d);


>    evalf(%);


Hardly different at all...

Let's try it again with a high- n  wavefunction, just to make sure that there's nothing special about the low-energy states. The general formula for the 1-d harmonic oscillator wavefunctions is

>    psi_1d := (n,x) -> (alpha/Pi)^(1/4)/sqrt(2^n*n!)*exp(-alpha*x^2/2)*HermiteH(n,sqrt(alpha)*x);

psi_1d := proc (n, x) options operator, arrow; (alpha/Pi)^(1/4)/sqrt(2^n*n!)*exp(-1/2*alpha*x^2)*HermiteH(n,sqrt(alpha)*x) end proc

Let's check that these are normalized properly:

>    simplify(int(psi_1d(0,x)^2,x=-infinity..infinity));


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



Now let's calculate < x ^2> for (say) n =5 and n =6. (We could in principle work with higher values of n , but the integrals get hard, even for Maple, in the 3-d case.)

>    avg_x2_1d_5 := simplify(int(psi_1d(5,x)^2*x^2,x=-infinity..infinity));

avg_x2_1d_5 := 11/2/alpha

>    avg_x2_1d_6 := simplify(int(psi_1d(6,x)^2*x^2,x=-infinity..infinity));

avg_x2_1d_6 := 13/2/alpha

The 3-d results are going to be a little work to get because we'll have to work out the normalization constant on the fly:

>    zeta_unnorm := (n,r) -> exp(-alpha*(r-req)^2/2)*HermiteH(n,sqrt(alpha)*(r-req))/r;

zeta_unnorm := proc (n, r) options operator, arrow; exp(-1/2*alpha*(r-req)^2)*HermiteH(n,sqrt(alpha)*(r-req))/r end proc

>    NFzeta := n -> 1/sqrt(int(zeta_unnorm(n,r)^2*r^2,r=0..infinity));

NFzeta := proc (n) options operator, arrow; 1/sqrt(int(zeta_unnorm(n,r)^2*r^2,r = 0 .. infinity)) end proc

>    zeta := (n,r) -> NFzeta(n)*zeta_unnorm(n,r);

zeta := proc (n, r) options operator, arrow; NFzeta(n)*zeta_unnorm(n,r) end proc


>    int(zeta(0,r)^2*r^2,r=0..infinity);


>    int(zeta(5,r)^2*r^2,r=0..infinity);


We can now calculate the expectation values:

>    avg_x2_3d_5 := int(zeta(5,r)^2*(r-req)^2*r^2,r=0..infinity):

>    avg_x2_3d_6 := int(zeta(6,r)^2*(r-req)^2*r^2,r=0..infinity):

We can now do the comparison using our typical values of alpha and req:

>    subs(alpha=alpha_typical,avg_x2_1d_5);


>    evalf(subs(alpha=alpha_typical,req=1e-10,avg_x2_3d_5));


>    subs(alpha=alpha_typical,avg_x2_1d_6);


>    evalf(subs(alpha=alpha_typical,req=1e-10,avg_x2_3d_6));


The agreement is even better than for n =0!