that both one or both Ψ (+∞) →0, Ψ (-∞) → 0.Or the function is finite at one boundary and zero at the other. The particle in the box, though simple, exhibits the features of all eigenvalue problems.

Using Schrodinger's equation solutions ( plural) are sought that satisfy the boundary conditions . In this case they become zero at the boundaries. From this condition a set of discrete energy eigenvalues E₁ , E ₂... E_n results accompanied by a set of corresponding eigenfunctions Ψ₁ , Ψ₂ ....Ψ_n .

One starts with the general statement H Ψ = E Ψ and ends with a set { E_i } , { Ψ_n } such that ,

H Ψ_n = E_n Ψ_n .

The infinite one dimensional well

Fig 1. Square well with infinte potential at the walls.

Let V= ∞ when x < 0 and x > L . V(x)=0 , 0<x<L and the boundary condition be

Ψ (0) = Ψ (L) =0. The walls are impenetrable because they have an infinite height.

Schrodinger equation is

d² Ψ /dx² = -(2mE/h'² ) Ψ = -k²Ψ . (1)

The two independent solutions are

Ψ (x) = A sin(kx) + B cos(kx) . (2)

Applying the boundary condition at the origin ,gives , Ψ (0) =0 = B . The cosine part drops out.

Applying the boundary condition at x= L , gives Ψ (L) ==0 = A sin(kL) which implies kL= nπ , which is the quantization condition.

The constant A will be determined from the normalization condition. It is not necessary to find the energy since

E_n = (h'² /(2m) ) (π/L)² n² , n=1,2,3..... (3)

The normalization condition requires that ∫ Ψ² dx = 1 = A² ∫ sin (kx)² dx , 0 < x < L

= A² ∫ ( 1-cos (2kx) )/2 dx = A² L/2 .

Thus the normalization constant is simply A = (2/L)^1/2 .The normalized eigenfunctions are

Ψ _n (x) = (2/L)^1/2 sin ( n π x/L) . (4)

The set Ψ _n (x) is orthonormalized ∫ Ψ _n (x) Ψ _m (x) dx = δ _{n,m .} ( δ _n,m=1 ,if n=m ;δ _n,m=0 for n≠m ).

The delta function potential

Some properties of the Delta function in one dimension are

∫ f(x) δ(x- x₀ ) dx = f (x ₀ ) , -∞ < x < +∞ (5)

∫ δ(x- x₀ ) dx =1 , (6)

and ∫ δ( a x ) dx =1/a . (7)

The dimensions of δ( x ) are 1/length .

A dirac sequence is established when certain functions approach the delta function with the help of an adjustable parameter. Two example are given. First consider the gaussian function

P(x) = (1/σ) ( 1/(2π)^1/2 ) exp[ -(x-μ)² /(2σ² )] , (8)

μ is the mean and σ² the variance. In other words

∫ x P(x) dx = μ , -∞ < x < +∞ (9)

and ∫ x² P(x) dx = σ² . . (10)

Fig 2. The gaussian function P(x) will resemble Dirac delta function as σ→ 0 .

Figure 2 shows that P(x) becomes sharply peaked as the variance ( σ ) goes to zero.

A second sequence that converges to the Dirac delta function is

(11)

It resembles the dirac delta function as n becomes large .

Its integral ∫ δ_n (x) dx = 1 , - ∞ < x < + ∞ .

Let the potential be V(x) =- η δ(x) where η is the strength of the delta function. This peculiar potential produces a single energy eigenvalue given by E = - {h'² /(2m)} η² . To show this start with

d² Ψ /dx² - (2m/h'² ) V(x) Ψ = - (2m/h'² ) E Ψ . (12)

Integrate the whole expression around the origin from -ε to +ε . We get

∫ (d² Ψ /dx² ) dx = ( d Ψ( +ε ) /dx - d Ψ( -ε ) /dx ) = 2 ( d Ψ( +ε ) /dx ) , (13)

see next figure 3.

Fig 3. Shape of the eigenfunction of the delta function.

The potential term integration gives

- (2m/h'² )∫ V(x) Ψ dx = - (2m/h'² )∫ -η δ(x) Ψ dx = (2m/h'² ) η Ψ(0) . (14)

The last term is zero , i.e.

- (2m/h'² ) E ∫ Ψ dx = - (2m/h'² ) E Ψ( 0) (2ε ) → 0 . (15)

From (13) and (140 it follows that

d Ψ( +ε ) /dx ) = (-1/2)(2m/h'² ) η Ψ(0) = -(m/h'² ) η Ψ(0) . (16)

The function near the origin ( x > 0) is thus of the form

Ψ ~ exp{-(m η /h'² ) x } . (17)

To get the energy go back to (12) at x near the origin (x>0) and use (17)

E Ψ = - {h'²/(2m) } d² Ψ /dx² = - (1/2) (m/h'² ) η² Ψ , (18)

E = -m η² /(2 h'² ) ,

or E = - η² /2 in units where m=1,h'=1.

This is the only allowed energy eigenvalue by an attractive delta function of strength η.

Figure 4 shows a delta like potential. Its eigenvalue is E = - η² /2 = -50.

Sage code for the plot

var('x')
mu=0 ; sigma=1.e-3 ;eta=10;
gauss= (1/sigma)*(1/(2*pi)^(1/2))*exp(-(x-mu)^2/(2*sigma^2));
v=-eta*gauss;
#v =x;
y=plot(v,x,-3*sigma,3*sigma)
show(y)

eta,-eta**2/2.= 0.1000E+02 -0.5000E+02

Fig 4. The potential V(x)= - η (1/σ) ( 1/(2π)^1/2 ) exp[ -(x-μ)² /(2σ² )] , η=10, μ=0 ,σ = 1.E-3.

We integrated Schrodinger equation numerically ,starting at x=0, with initial conditions Ψ(0)=1 , dΨ(0) /dx =-1. The integration is carried out up to x_final= 250 σ for different values of the energy. The eigenvalue becomes apparent when Ψ ( x_final) changes sign.

e,psif= -0.9667E+02 0.4139E+01
e,psif= -0.9333E+02 0.3697E+01
e,psif= -0.9000E+02 0.3279E+01
e,psif= -0.8667E+02 0.2885E+01
e,psif= -0.8333E+02 0.2513E+01
e,psif= -0.8000E+02 0.2162E+01
e,psif= -0.7667E+02 0.1832E+01
e,psif= -0.7333E+02 0.1522E+01
e,psif= -0.7000E+02 0.1230E+01
e,psif= -0.6667E+02 0.9578E+00
e,psif= -0.6333E+02 0.7027E+00
e,psif= -0.6000E+02 0.4646E+00
e,psif= -0.5667E+02 0.2428E+00
e,psif= -0.5333E+02 0.3659E-01
e,psif= -0.5000E+02 -0.1546E+00
e,psif= -0.4667E+02 -0.3315E+00
e,psif= -0.4333E+02 -0.4946E+00
e,psif= -0.4000E+02 -0.6446E+00
e,psif= -0.3667E+02 -0.7820E+00
e,psif= -0.3333E+02 -0.9073E+00
e,psif= -0.3000E+02 -0.1021E+01
e,psif= -0.2667E+02 -0.1124E+01
e,psif= -0.2333E+02 -0.1216E+01
e,psif= -0.2000E+02 -0.1299E+01
e,psif= -0.1667E+02 -0.1371E+01
e,psif= -0.1333E+02 -0.1435E+01
e,psif= -0.1000E+02 -0.1490E+01
e,psif= -0.6667E+01 -0.1537E+01
e,psif= -0.3333E+01 -0.1576E+01
e,psif= 0.0000E+00 -0.1608E+01

c FORTRAN code for the numericalintegration of delta function potential
implicit real*8(a-h,o-z)
real*8 k
dimension psif(200),energy(200),steps(200)
data ne,epsi,nx/30,1.d-2,2000/
gauss(x) =(1.d0/(sigma*dsqrt(2.d0*pi)))*
$ dexp(-(x-amu)**2/(2.d0*sigma**2))
V(x)=-eta*gauss(x)
c v(u)=-eta*delta(u,epsi)
pi=2.d0*asin(1.d0)
amu=0.d0
sigma=1.d-3
eta=10.d0
c ei= -eta/(2.d0*epsi)
ei=-eta**2
ef=0.d0
de=(ef-ei)/dfloat(ne)
do 20 ie=1,ne
e=ei + de*dfloat(ie)
k=dsqrt(-2.d0*e)
xi=0.d0
c xf= -(1.d0/k)*log(1.d-2)
xf=250.d0*sigma
c dx=epsi/10.d0
dx=(xf-xi)/dfloat(nx)
c epsi=dx
kp=int(float(nx)/40)
kount=kp
psi0=1.d0
deriv=-1.d0
psi1=psi0+dx*deriv
do 10 i=2,nx
x=xi+dx*dfloat(i)
psi2=2.d0*psi1-psi0 -dx**2*(2.d0*(E-V(x-dx))*psi1)
c psi2=2.d0*psi1-psi0 -dx**2*(2.d0*(E-V(x-dx,eta,epsi))*psi1)
c if(i.eq.kount)then
c print 130,x,psi2
c kount=kount+kp
c endif
psi0=psi1
psi1=psi2
10 continue
energy(ie)=e
psif(ie)=psi2
20 continue
print 120 ,eta,-eta**2/2.d0
print*,' '
do 30 i=1,ne
print 100,energy(i),psif(i)
30 continue
100 format('e,psif=',3(3x,d11.4))
120 format('eta,-eta**2/2.=',2(3x,d11.4))
130 format('x,ps=',2(3x,d10.3))
stop
end

Fig 5. Eigenfunction (not normalized) for the delta function potential.

The Linear Oscillator

The problem of the harmonic oscillator is one of paramount importance in quantum mechanics.

For the one dimensional oscillator the hamiltonian is

H = - {h'²/(2m)} d²/dx² + (1/2) k x² , (19)

accordingly

- {h'²/(2m)} d² Ψ/dx² + (1/2) k x² Ψ = E Ψ . (20)

scale units h'²/(mL²) ~ k L² , L_scale = ( h'² /(mk) )^1/4 (21)

E_scale ~ k (L_scale )² = h' (k/m)^1/2 = h' ω , (22)

where ω is the classical frequency of the oscillator.

It was shown in chapter two that the energy dependence on the quantum number n has the form E_n ~ n^2α/(2+α) if V(x) ~ x^α .

For the harmonic oscillator α =2 and E_n ~ n , or more explicitly E_n ~ n h^' ω + constant term.

We cannot physically have simply E_n ~ n h^' ω , because n=0 would be satisfied only by Ψ(x)=0. Thus a constant term must be present in the formula for E_n .

( See Ref http://www1.uprh.edu/rbaretti/LQMch2.htm for a discussion of the dependency of E upon the quantum number n.)

Replacing h'=1,m=1,k=1 in (20) transforms the equation in

- (1/2) d² Ψ/dx² + (1/2) x² Ψ = E Ψ . (23)

The set of eigenfunctions Ψ_nis an orthogonal set. Let H Ψ_n = E_n Ψ_n and H Ψ_m = E_mΨ_m . Multiply the first by

Ψ_m and the second by Ψ_{n .}

Substract ∫ ( Ψ_m H Ψ_n -Ψ_n H Ψ_m ) dx = 0 = ( E_n - E_m)∫ Ψ_m Ψ_n dx . Therefore

∫ Ψ_m Ψ_n dx =0 .

We will seek the first eigenvalues and eigenfunctions using what is called the creation operator.

When (1/2)x² >>E, the asymptotic form of (23) is - (1/2) d² Ψ/dx² + (1/2) x² Ψ ≈ 0 and the approximate solution for x large (x>>1) is

Ψ_asymptotic = exp(-x²/ 2) . (24)

- (1/2) d² Ψ_asymptotic /dx² + (1/2) x² Ψ_asymptotic = (1/2) exp(-x² /2)→ 0 as x → +∞ or -∞ .

Consider now the substitution of the normalized function Ψ = (1/π)^1/4 exp(-x² /2) ≡ N exp(-x² /2) in eq (23).

The left side is

- (1/2) d² Ψ/dx² + (1/2) x² Ψ = (1/2) (1/π)^1/4 exp(-x² /2) =(1/2) Ψ . (25)

Equated to E Ψ gives E₀ =1/2 the energy of the ground state. Hence the complete expression for the eigenvalues is

E _n = (n +1/2) h'ω , or

E_n = (n+1/2) (26)

Figure 6 shows the normalized ground state wave function Ψ₀ = (1/π)^1/4 exp(-x² /2).

Fig 6 . Ground state wave function of the harmonic oscillator.

Part a: Creation and annihilation operators

The first excited state has n=1 , E= 3/2 and can be obtained with the rising or creation operator defined by

a+ = (1/2^1/2 ) ( -i p +x) = (1/2^1/2 ) ( - d/dx + x ) . (27)

Operating with with the rising operator a⁺ upon the ground state ,

a+ Ψ₀ = (1/2^1/2 ) ( - d/dx + x ) (1/π)^1/4 exp(-x² /2) = ( 2^1/2 / π^1/4 ) x exp(-x²/2) = Ψ_plus , (28)

the subscript label (plus) attached to psi is provisional.

SAGE CODE

psi0= (1/pi^(1/4))*exp(-x^2/2);
psiplus=(1/2^(1/2) )*(-diff(psi0,x)+x*psi0 );

psiplus

sqrt(2)*x*e^(-1/2*x^2)/pi^(1/4)

We integrate next psi squared square, (Ψ_plus )² , finding ∫ (Ψ_plus )² dx = 1. , -∞ < x < + ∞ .

Sage code

psiplus =sqrt(2)*x*e^(-1/2*x^2)/pi^(1/4)

integral(psiplus^2,x,-oo,+oo)

Evaluate now the average energy E = ∫ Ψ_plus H Ψ _plusdx / ∫ (Ψ_plus )² dx

= ∫ Ψ_plus{- (1/2) d² /dx² + (1/2) x² } Ψ _plusdx

= 3/2 . (30)

When (30) is equated to E_n we have (n+1/2) = 3/2 ; n=1.

Function Ψ_plus is renamed Ψ₁ and corresponds to the first excited state of the quantum oscillator,

Ψ₁ = ( 2^1/2 / π^1/4 ) x exp(-x²/2) , E₁ =3/2 . (31)

In conclusion, the effect of the rising operator upon Ψ₀ creates a function proportional ,or equal , to Ψ₁ .

Fig 7. Normalized Ψ₁ .

Evidently ∫ Ψ₁Ψ₀ dx= 0 since the integrand is an odd function in the interval -∞ < x < + ∞ .

Applying a⁺ on Ψ₁ , produces a⁺ Ψ₁ = Ψ_plus = (1/π^1/4 )*(2x² -1)*exp(-x²/2) . To normalize Ψ_plus , we find

∫ (Ψ_plus )² dx = ∫ {(1/π^1/4 )*(2x² -1)*exp(-x²/2) } ² dx= 2 .

It turns out that the function obtained from Ψ_plus / 2^1/2 = (1/2^1/2) (1/π^1/4 )*(2x² -1)*exp(-x²/2) , is Ψ₂ as we verify next. Its energy is

E = ∫ Ψ₂ {- (1/2) d² /dx² + (1/2) x² } Ψ₂ dx

= 5/2 = ( 2 + 1/2) , (32)

i.e. n=2 .

So Ψ₂ = (1/2^1/2) (1/π^1/4 )*(2x² -1)*exp(-x²/2) , E₂ = (2 +1/2) . (33)

Sage code

f=sqrt(2)*x*e^(-1/2*x^2)/pi^(1/4);
(1/2^(1/2))*(-diff(f,x) +x*f)

psiplus= (1/pi^(1/4))*(2*x^2 -1)*exp(-x^2/2);
integral(psiplus^2,x,-oo,+oo)

Sage code energy eigenvalue (n=2)

psi2=(1/2^(1/2))*(1/pi^(1/4))*(2*x^2 -1)*exp(-x^2/2);
integral(psi2*( (-1/2)*diff(psi2,x,2) +(1/2)*x^2*psi2),x,-oo,+oo)

psi2=(1/2^(1/2))*(1/pi^(1/4))*(2*x^2 -1)*exp(-x^2/2);
psi3=(1/3^(1/2))*(1/2^(1/2))*(-diff(psi2,x) +x*psi2);
#integral(psi3^2,x,-oo,+oo)
# 1
integral(psi3*( (-1/2)*diff(psi3,x,2) +(1/2)*x^2*psi3),x,-oo,+oo)
# 7/2

Sage code

f=sqrt(2)*x*e^(-1/2*x^2)/pi^(1/4);
(1/2^(1/2))*(-diff(f,x) +x*f)

psiplus= (1/pi^(1/4))*(2*x^2 -1)*exp(-x^2/2);
integral(psiplus^2,x,-oo,+oo)

Sage Code

psiplus=sqrt(2)*x*e^(-1/2*x^2)/pi^(1/4);
hpsiplus= (-1/2)*diff(psiplus,x,2) +(1/2)*x^2*psiplus;
integral(psiplus*hpsiplus,x,-oo,+oo)

3/2

Do it one more time. Apply a⁺ upon Ψ₂ . Guided by the previous example, we suppose that

a⁺ Ψ₂ = (3)^1/2 Ψ₃ = Ψ_plus . (34)

Hence Ψ₃ = (1/3^1/2 )(1/π^1/4 ) *exp( -x²/2)*( 2 x³ - 3x ) . (35)

We verify that ∫ Ψ² dx =1 ; it is normalized.

psi3= (1/3^(1/2) )*(1/pi^(1/4))*exp( -x^2/2)*( 2*x^3 - 3*x );
integral(psi3^2,x,-oo,oo)

The energy eigenvalue is

E₃ = ∫ Ψ₃ {- (1/2) d² /dx² + (1/2) x² } Ψ₃ dx

= 7/2 = ( 3+1/2) . (36)

d2=-(1/2)*diff(psi3,x,2)
integral(psi3*(d2+(1/2)*x^2*psi3),x,-oo,oo)

Thus applying the creation operator on a normalized eigenfunction gives,

a⁺ Ψ_n = (n)^1/2 Ψ_{n+1
,}(37)

or in a more abstract form

a⁺ │ n> = (n+1)^1/2 │ n+1 > , (38)

The vectors │ n> are an abstraction whose projection in x space is written as

Ψ_n (x) = < x│n> . Alternatively the function in p -space (momentum ) is

φ_n (p) = < p│n> .

From the Fourier transform defintion φ_n (p) = ∫ dx exp(-px) Ψ_n (x) , this can be rewritten as

φ_n (p) = ∫ dx <p│ x> <x │n > .

The term <p│ x> =exp(-ipx) and <x│ p> = exp( +ipx) ( remember we use h'=1) .

The annihilation operator , a , is the conjugate of a⁺ ,

a ≡ (1/2^1/2 ) ( i p +x) = (1/2^1/2 ) ( d/dx + x ) . (39)

One readily verifies that

a Ψ₀ =(1/2^1/2 ) ( d/dx + x )(1/π)^1/4 exp(-x² /2)=0 , (40)

psi0=(1/pi^(1/4))*exp(-x^2/2) ;
(1/2^(1/2))*(diff(psi0,x)+x*psi0)

and a Ψ₁ = (1/2^1/2 ) (d/dx + x ) ( 2^1/2 / π^1/4 ) x exp(-x²/2)

= (1/π)^1/4 exp(-x² /2) = Ψ₀ . (41)

psi1=(2^(1/2)/pi^(1/4))*x*exp(-x^2/2);
(1/2^(1/2))*(diff(psi1,x)+x*psi1)

e^(-1/2*x^2)/pi^(1/4)

Also a Ψ₂ = (2 /π^1/4 ) x exp(-x²/2) = 2^1/2 Ψ₁ (42)

and in general a Ψ_n = n^1/2 Ψ_n-1 , (43)

or a│ n> = (n)^1/2 │ n-1 > . (44)

2*x*e^(-1/2*x^2)/pi^(1/4)

Consider the average < n │a+ a│n> = < n │a+ n^1/2 │n-1 > = n^1/2 < n │a+ │n-1 >

= n^1/2 < n │ (n-1+1)^1/2 │n-1+1 >

= n < n │n > = n (45)

Thus

E_n = n + 1/2 = <n │ a+ a +1/2│ n> , (46)

= <n │ H│ n> . (47)

= H_{n,n ,}all non diagonal values are zero.

The physical dimensions of these particular operators , a⁺ and a, are energy^1/2 ~ (h' ω)^1/2 = h^'1/2(k/m)^1/4 .

In many fields of physics the use of creation and annihilation operators is preferred over the use of wave functions. The technique is also called second quantization.

Example: Find the matrix elements of a , a + , x , p , x² and p².

a_n,m= m^1/2 δ _n,m-1 . (48)

The only non zero elements are those when the column index is one more than the row index.

Thus

m = 0 1 2 3 4

a_{n,m =} 0 1 0 0 0 .......

0 0 2^1/2 0 0

0 0 0 3^/2 0

and <n │a⁺│ m > = <n│(m+1)^1/2 │ m+1 >

a⁺_n,m = (m+1)^1/2 δ _n,m+1 . (49)

m = 0 1 2 3 4

a⁺_{n,m =} 0 0 0 0 0 .......

1 0 0 0 0

0 2^1/2 0 0 0

0 0 3^1/2 0 0

Evidently <n │a+ a │ m> = 0 0 0 0 0 ......

0 1 0 0 0

0 0 2 0 0

0 0 0 3 0

The matrix elements of x are

<n │x │ m> = (1/2)^1/2 ( a+ + a ) _n,m

0 1 0 0 0 .......

= (1/2)^1/2 1 0 2^1/2 0 0

0 2^1/2 0 3^1/2 0

0 0 3^1/2 0 4^1/2

The matrix elements of p are

<n │p │ m> = (i/2)^1/2 ( a+ - a ) _n,m

0 -1 0 0 0

= (i/2)^1/2 1 0 -2^1/2 0 0

0 2^1/2 0 - 3^1/2 0

0 0 3^1/2 0 - 4^1/2

<n │x² │ m> = <n │p² │ m> = (1/2) 1 0 0 0.......

0 3 0 0....

0 0 5 0......

0 0 0 7....

For more details on creation and annihilation operators see ,

http://en.wikipedia.org/wiki/Creation_and_annihilation_operators

Part b: Hermite polynomials

We saw above that exp(-x²/2) is both the asymptotic solution and the ground state wave function ( except for a

normalization constant). Now write

Ψ(x) = f(x) exp(-x²/2) , (50)

and seek the differential equation for the function f(x) . The solution of f (x) is sought by series method with the proviso that it is a finite polynomial. Susbstitution of (50) in (23) gives ,

d² f/dx² -2x df/dx -( 1-2E)f = 0 . (51)

Suppose first f (x) = x^s ∑_n=0 a_n xⁿ . Near the origin , x^s is the leading power. Let f ~ x^s in (51) ,

s(s-1) x^s-2 - 2 sx^s -( 1-2E)x^s ~ 0 .

The leading power near the origin is x^s-2 thus either s=0 or s= +1 . As will be shown consequence this gives rise to even or odd power series.Let s=0 , then f (x) = ∑_n=0 a_n xⁿ and

d² f/dx² = ∑_n a_n n(n-1) x^n-2 (52-a)

-2x df/dx =∑_n -2 n a_n xⁿ (52-b)

-( 1-2E)f= -(1-2E)∑a_n xⁿ . (52-c)

Alternatively setting s=1 and f(x) = ∑_n=0 a_n xⁿ⁺¹ leads essentially to the same algebraic relations as above.

The factors of the xⁿ power must ad up to zero , i.e. letting n→ n+2 in term (52-a) we get

a_n+2(n+2)(n+1) + ( - 2n -(1-2E) ) a_n = 0 , producing the recursion formula

a_n+2 = {2n +(1-2E)} a_n /((n+2)(n+1)) . (53)

If a₀ ≠ 0 the recursion relation (53 ) , yields the even power series a₂ , a₄ , a₆ ..... .

If a₁≠ 0 the recursion relation , yields an odd power series a₃ , a₅ , a_{7 .....}

To terminate the series at the nth power one must set 2n +(1-2E) = 0 , producing the quantization condition,

E_n = n +1/2.

Example : Let n=0 , a₀=1 , E =1/2 ; then a₂ = 0 . The function f(x) = a₀ and the un normalized Ψ is

Ψ ₀ = exp( -x²/2) . (54)

Let n=1 , a₁ =1 , E=3/2 ;then a₃ = 0 , f(x) = a₁x = x . The un normalized Ψ is

Ψ ₁= x exp( -x²/2) . (55)

Let n=2 , a₀ =1 , E=5/2 ;then a₂ = -2 , f(x) = a₀+ a₂ x² = 1-2x² . The un normalized Ψ is

Ψ ₂= ( 1-2x² ) exp( -x²/2) . (56)

The first four Hermite polynomials (see e.g Ref 2) , in one of two possible conventions, are written as

H₀(x) =1 ; H₁(x) = 2x ; H₂(x) = 4 x² -2 , H₃ (x)= 8 x³-12 x (57)

which can be obtained from the formula

H_n(x) = (-1)ⁿ exp( x²) { dⁿ exp(-x²) /dxⁿ } . (58)

Our function f(x) of these examples is proportional to a Hermite polynomial. The wavefunctions of the harmonic oscillator are of the form

Ψ_n(x) = N_n H_n(x) exp(-x²/2) , (59)

where one can find by induction that the normalization constant is N_n= (1/π^1/4) {1/(2ⁿ n! )}^1/2 .

Exercise : Justify the expression for N_n by normalizing successively Ψ ₀ ,Ψ ₁ ,Ψ ₂ , Ψ ₃ etc.

First calculate the integrals ∫ (Ψ _n )² dx to find out what pattern is established.

The normalization constant of Ψ₀ , is obtained from,

∫ (Ψ ₀ )² dx = ∫ (H₀ exp(-x^2/2) )² dx = π^1/2 = 1/ (N₀ )² , thus N₀ = 1/π^1/4

SAGE code

h0=1; psi0=h0*exp(-x^2/2);
s0=integral(psi0^2,x,-oo,oo);
n0=1/s0^(1/2);
s0 , n0

The normalization constant of Ψ₁ , is obtained from ,

∫ (Ψ ₁ )² dx = ∫ (H₁ exp(-x^2/2) )² dx π^1/2 = S₁ = 2 π^1/2 = 1/(N₁)² ; thus N₁ =(1/ π^1/4 ) 1/(2)^1/2 .

The normalization constant of Ψ₂ , is obtained from ,

∫ (Ψ ₂ )² dx = ∫ (H₂ exp(-x^2/2) )² dx π^1/2 = S₂ = 8 π^1/2 = 1/(N₂)² ; thus N₂=(1/ π^1/4 ) 1/( 2² *2)^1/2 .

One last normalization constant .The normalization constant of Ψ₃ , is obtained from

∫ (Ψ ₃ )² dx = ∫ (H₃ exp(-x^2/2) )² dx π^1/2 = S₃ = 48 π^1/2 = 1/(N₃)² ; thus N₃ =(1/ π^1/4 ) 1/( 2³ *3!)^1/2 .

Only the last two constants finally give away the formula

N_n = (1/ π^1/4 ) 1/( 2ⁿ * n!)^1/2 . (58)

The normalized eigenfunctions in adimensional units are

Ψ ₂ (x) = {(1/ π^1/4 ) 1/( 2ⁿ * n!)^1/2 } H_n (x) exp(-x²/2) . (59)

The following SAGE code can be used to plot some solutions (see http://www.sagenb.org/home/pub/590/ ).

x = var('x')
def u(n):
fac = ((2*factorial(n))^(-1/2))*(2^(1/4))
return fac * hermite(n, sqrt(2*pi)*x) * exp(-pi*x^2)

u_0(x) = u(0)
u_2(x) = u(2)
u_4(x) = u(4)

graph1 = plot(u_0(x), -3, 3, rgbcolor=(1, 0, 0))
graph2 = plot(u_2(x), -3, 3, rgbcolor=(0, 1, 0))
graph3 = plot(u_4(x), -3, 3, rgbcolor=(0, 0, 1))

show(graph1 + graph2 + graph3)

Part c: Numerical solution ( last but not least)

Equation (23) becomes a finite difference equation when the second derivative is approximated by

d² Ψ /dx² ≈ (Ψ_n -2Ψ_n-1+ Ψ_n-2 )/ (∆x)² . (60)

The unit of length is L_scale = ( h'² /(mk) )^1/4 ≡ 1 . In the integration one must choose ∆x << L_scale =1 .

The solution to

- (1/2) d² Ψ/dx² + (1/2) x² Ψ = E Ψ is,

Ψ_n = 2Ψ_n-1- Ψ_n-2 - 2 (∆x)² { E- (1/2) (x_n-1)² }Ψ_{n-1
,}(61)

where n is not an eigenvalue number but a supscript, Ψ_n = Ψ(x_n) , Ψ_n-1 = Ψ(x_n-1) = Ψ(x_n - ∆x) ,

Ψ_n-2 = Ψ(x_n-2) = Ψ(x_n - 2∆x) .

The initial values are arbitrary ,

Ψ₀ = Ψ(0) = Ψ₁ =Ψ( ∆x) =1 , ie. (dΨ/dx)_x=0 =0 for even functions like the ground state , the second excited state ,the fourth excited state, etc.

For odd functions Ψ₀ = Ψ(0) = 0 , (dΨ/dx)_x=0 = 1 , Ψ₁ = Ψ₀ + ∆x .

A FORTRAN code is provided so that Eq (61) is integrated for a variable E until x_asymtotic >> x_{classical
turning point} . The classical turning point is defined from E = V(x) = (1/2) x² ; x_{classical turning point} = (2 E) ^1/2 . The "asymptotic " Ψ (x_asymtotic ) flips from positive to negative or vice versa when it passes through an energy eigenvalue.

Fig 8 . Even states energy eigenvalues are found at E= 0.5 , 2.5, 4.5, 6.5 .

Fig 8 . Odd states energy eigenvalues are found at E= 1.5 , 3.5, 5.5 .

FORTRAN code

c Quantum harmonic oscillator
dimension psinf(100),energy(100)
data niter,nstep/80,2000/
v(x)=.5*x**2
e=.05
efinal=7.
de=(efinal-e)/float(niter)
do 110 iter=1,niter
xi=0.
xlim=2.5*sqrt(2.*e)
dx=xlim/float(nstep)
kp=int(float(nstep)/60.)
kount=kp
c IC for even functions
psi0=1.
psi1=psi0
c IC for odd functions
c psi0=0.
c psi1=psi0+dx
c print 200 ,0.,psi0,psi0/gauss(0.)
do 100 i=2,nstep
x=xi+dx*float(i)
psi2=2.*psi1-psi0 -dx**2*2.*(E-v(x-dx))*psi1
if(i.eq.kount)then
c print 200 ,x, psi2,psi2/gauss(x)
c print 200 ,-x,-psi2,-psi2/gauss(x)
kount=kount+kp
endif
psi0=psi1
psi1=psi2
100 continue
energy(iter)=e
psinf(iter)=psi2
20 continue
e=e + de
110 continue
do 30 i=1,niter
print 120 ,energy(i),psinf(i)/abs(psinf(i))
30 continue
120 format(3x,'energy,psinf/abs(psinf)=',2(3x,e10.3))
200 format(3x,'x,psi,h=',3(3x,e10.3))
stop
end

Finite Square Well Potential

A finite square well (next figure) can be defined by

V(x) =0 , -L < x < L

V(x) = V₀ , x>L , x < -L.

The solutions can be classified according to their parity.

Starting with the ground state a set of solutions are symmetric Ψ(x) = Ψ(-x). Another set starting with the first excited state are anti symmetric , Ψ(x) = - Ψ(-x) .

The symmetric states ( even functions) will "resemble" cos (k x) where k =(2E)^1/2 and k L ≈ (2n + 1)π/2 . It would be an equality if we had an infinite well. The anti symmetric states ( odd functions) will resemble sin(kx) and kL ≈ nπ .In the regions /x/ >L the functions are decaying exponentials. Depending on the height V₀ only a finite number of energy eigenvalues is permitted.

Fig 9 .Sample of class notes about the finite square potential.

For even states

Ψ ₁ (x) = A exp (k₁x) , x< -L , k₁ = ( 2 (V₀ -E) )^1/2 , (62)

Ψ ₂ (x) = B cos (k₂x) , -L < x < L , k₂ = ( 2 E )^1/2 , (63)

Ψ ₂ (x) = A exp (-k₁x) , x > L .

From the continuity conditions on Ψ and dΨ /dx , at x= -L, we get ,

A exp(-k₁L) = B cos (k₂L) (64)

k₁ A exp(-k₁L) = k₂ B sin (k₂L) (65)

Dividing eq (65) by eq (64) gives the transcendental equation ,

k₁ = k₂ tan ( k₂L) . (66)

The solutions can be obtained by graphical or numerical means .

For odd states the continuity equations ,(at x = -L) , lead to

A exp(-k₁L) = - B sin (k₂L)

k₁ A exp(-k₁L) = k₂ B cos (k₂L)

and the transcendental equation is now

k₁ = - k₂ cot (k₂L) . (67)

Figures (10) and (11) show the allowed energies (three) for a finite square well with V₀ = 8.

Fig 10 . Even state energies are found near E = 0.80 , 6.4 .

A refinement of the energies is given by Newton's formula for the roots. Define

f (E) = k₁(E) - k₂ (E) tan ( k₂(E)*L) .

The roots are found from the iteration formula

E_n = E_n-1 - f (E_n) / ( df/dE)_n . (68)

Ten iterations for the ground state E =0.784194

i, e= 1 0.5
i, e= 2 0.939939141
i, e= 3 0.838132858
i, e= 4 0.790814519
i, e= 5 0.784336388
i, e= 6 0.784195423
i, e= 7 0.784194291
i, e= 8 0.78419435
i, e= 9 0.784194231
i, e= 10 0.784194291

Fortran code for Newton's roots formula

c roots of k1 = k2 tan ( k2*L)
real L , lscale ,k1,k2
data L ,V0 /1.,8. /
k1(e) = (2.* (V0 -E) )**(1./2.)
k2(e)=(2.*e)**(1./2.)
c f(e) for even states
f(e)=k1(e)-k2(e)*tan(k2(e)*L)
c f(e) for odd states
c f(e)= k1(e)+ k2(e)/tan(k2(e)*L)
df(e)= (f(e+epsi)-f(e))/epsi
Lscale=1/v0**.5
epsi=lscale/100.
c initial E must be selected near an eigenvalue
e=0.5
do 10 i=1,10
e1=e - f(e)/df(e)
print*,'i, e= ', i, e
e=e1
10 continue
stop
end

Fig 11. Odd states , only one energy is allowed ; E ≈ 3.0 .

The first three energy eigenvalues for the infinite square well of size 2L (letting L=1) are ,

E_n = { π² /(8 L² ) } n² = 1.23 n² = 1.23 , 4.92 , 11.07 . Since V₀ is finite , the allowed penetration of Ψ into the classically forbidden region represents a larger lambda and consequently a lower energy.

c energies for the finite square well
real k1,k2, L
data L,ne /1.,80/
k1(e) = (2.*(V0 -E))**(1./2.)
k2(e)=(2.*E)**(1./2.)
v0=8./L**2
print*,'Eisberg,.098*v0,.8*v0=',.098*v0,.8*v0
print*,' '
ei=0.01
ef=v0
de=(ef-ei)/float(ne)
e=ei
do 10 i=0,ne
even=k2(e)*tan(k2(e)*L)
odd= k2(e)/tan(k2(e)*L)
if(abs(even).lt.8.) print 100,e, k1(e) ,even
e=e+de
10 continue
100 format('e,k1,even=',3(3x,e10.3))
stop
end

A refinement of the energies

Finite square well by numerical integration

Fig 12. Approximate energies of the first two and only even states are E ≈ 0.8 , 6.76 . Compare with the estimate in figure 10.

Fig 13 . Energy of the odd state ≈ 3.075 . The four digits are read from the data results.

Applying Newton's method ,eq. (68) we get E= 3.061.

i, e= 1 2.5
i, e= 2 3.18810344
i, e= 3 3.06919932
i, e= 4 3.06180215
i, e= 5 3.06176519
i, e= 6 3.06176543
i, e= 7 3.06176496
i, e= 8 3.06176519
i, e= 9 3.06176543
i, e= 10 3.06176496

c numerical integration of finite square well
dimension psinf(100),energy(100)
real L
data L,v0 , niter,nstep/1.,8.,80,2000/
e=.05
efinal=v0
de=(efinal-e)/float(niter)
do 110 iter=1,niter
xi=0.
xlim=1.6*L
dx=xlim/float(nstep)
kp=int(float(nstep)/60.)
kount=kp
c IC for even functions
psi0=1.
psi1=psi0
c IC for odd functions
c psi0=0.
c psi1=psi0+dx
c print 200 ,0.,psi0,psi0/gauss(0.)
do 100 i=2,nstep
x=xi+dx*float(i)
psi2=2.*psi1-psi0 -dx**2*2.*(E-v(x-dx,L,v0))*psi1
c if(i.eq.kount)then
c kount=kount+kp
c endif
psi0=psi1
psi1=psi2
100 continue
energy(iter)=e
psinf(iter)=psi2
20 continue
e=e + de
110 continue
do 30 i=1,niter
print 120 ,energy(i),psinf(i)/abs(psinf(i))
30 continue
120 format(3x,'energy,psinf/abs(psinf)=',2(3x,e10.3))
200 format(3x,'x,psi,h=',3(3x,e10.3))
stop
end

function v(x,L,v0)
real L
if(x.lt.L)v=0.
if(x.ge.L)v=v0
return
end

END OF CHAPTER