PHYS 2063 Wave Physics:
Wave Mechanics of the One Dimensional Square Box with Finite Walls

Remove["Global`*"]

We are pretty much following Pain pages 434 and following. This notebook, though, actually carries through the numerical calculation of a few eigenvalues, and shows the shape of the resulting wave functions.

The Potential Energy Function

The potential energy function is zero for x between 0 and a, and V(>0) otherwise.
We won’t actually use this to solve the Schrodinger Equation (it’s just too easy to do piecewise),
but this will be useful for displaying the results.

pe=V(1-HeavisideTheta[x]+HeavisideTheta[x-a]);
avVals={V->2,a->3};
Plot[{pe/.avVals,0.8},{x,-3,6}]

The plot is drawn with a straight line at an arbitrary value to show a possible energy eigenvalue. Inside the well, the Schrodinger equation says the second derivative of the wave function is negative, so the wave function is sinusoidal.  Outside the well, the second derivative is positive, so the wave function is an exponential.

The Schrodinger Equation

The equation is

-(ℏ^2)/2m (d^2)ψ/(dx^2)+U(x)ψ=Eψ

where U(x)=0 for x between 0 and a, and U(x)=V otherwise. It is natural to scale x by a, that is x=ay, in which case

d^2 ψ/dy^2=-E/(ℏ^2/2(ma^2))ψ inside the well, or
d^2 ψ/dy^2=(V-E)/(ℏ^2/2(ma^2))ψ outside the well

So, if we just agree to measure x in units of a, and E and V in units of ℏ^2/2(ma^2), and then recognize that 0<E<V for bound states, then we have very simple solutions inside and outside the well. To do this use

E = (ℏ^2)(k^2)/2m and V-E = (ℏ^2)(κ^2)/2m where k and κ are both positive quantities.

Set up the solutions including the different regions

Region 1 is for 0<x<a: The wave function is sinusoidal

ψ1=A1 Cos[k x]+B1 Sin[k x];

Region 2 is for x>a: The wave function is a decaying exponential

ψ2=B2 Exp[-κ x];

Region 3 is for x<0: The wave function is an increasing exponential

ψ3=A3 Exp[κ x];

Match the boundary conditions to find equations relating the coefficients

Equate the functions and derivatives at x=0 and x=1

eq1=(ψ1/.x->0)==(ψ3/.x->0)

eq2=(D[ψ1,x]/.x->0)==(D[ψ3,x]/.x->0)

eq3=(ψ1/.x->a)==(ψ2/.x->a)

eq4=(D[ψ1,x]/.x->a)==(D[ψ2,x]/.x->a)

Can we solve these for the coefficients A1, A3, B1, and B2?

Solve[{eq1,eq2,eq3,eq4},{A1,A3,B1,B2}]

Of course. Four homogeneous equations for four unknowns gives all zeros. We must avoid the “boring” solution. Or, speaking more quantum mechanically, we have to find a normalizable solution, so all four coefficiencts cannot be zero.

Using linear algebra to find the relationship between k and κ

We’ve seen this before, in the ansatz we use to solve the multiple masses and spring problem. Think of the problem in terms of a matrix equation, where we’ll need to set the determinant equal to zero so that the four equations are no longer independent.

Mathematica makes it easy to form the matrix equation for the coefficients

CoefficientArrays[{eq1,eq2,eq3,eq4},{A1,A3,B1,B2}];
mat=Part[Normal[%],2]

Take a look at the determinant

Det[mat]

Setting the determinant equal to zero will give us some conditions on k and κ. We will scale the two wavenumbers by the width of the well, so they are replace by dimensionless quantitites b and β.

kβReps={k->b/a,κ->β/a};
eqx=Simplify[(Det[mat]/.kβReps)==0,{β>0,a>0}]

Examine the solution of the resulting equation

It will be easier to interpret the result if we write this equation in terms of r=β/b=κ/k

exCoef=Simplify[eqx/.β->r b,b>0]

There are two solutions for r in terms of b

Solve[exCoef,r]

Separte the two solutions for r, and call them r1 and r2

r1=r/.Solve[exCoef,r][[1]];
r2=r/.Solve[exCoef,r][[2]];

Now plot the solutions. Note that both b and β are positive (because k, κ, and a are positive), so r is also positive

Plot[{r1,r2},{b,0,5.1Pi},PlotRange->{0,10}]

(The vertical lines mark the points where tan and cot are not defined.) Note that r goes to ∞ when b=ka=π, 2π, 3π,...

Using scaling to interpret the solutions

In order for there to be normalizable solutions for A1, B1, A3, and B2, r=β/b=κ/k must be one of the following

r1
r2

Recall our definitions of the wavenumbers k and κ (with energies in units of ℏ^2/2(ma^2))

kκReps={k->Sqrt[energy]/a,κ->Sqrt[V-energy]/a};

Scale the energy eigenvalue using ϵ=energy/V and use this to rewrite r and b. Note that ϵ must be between zero and one.

z=Simplify[(κ/k)/.kκReps,m>0];
rScaled=Simplify[z/.energy->ϵ V,V>0]

bScaled=k a/.kκReps/.energy->ϵ V

Equating r and either tan(b/2) or -cot(b/2), for a given V (which is in units of ℏ^2/2(ma^2)) gives an equation for ϵ, which in turn gives the energy eigenvalue ϵV. Note that if ϵ is very small, that is the walls are very high, then r gets very large and b=ka=π, 2π, 3π,... as we saw on the last slide. For these values of k we indeed recover the eigenvalues for the infinite high walled box.

Find the energy eigenvalues

First examine the solution graphically. Use V=1 to begin.

Vval=1;
Show[
Plot[rScaled,{ϵ,0,1},PlotRange->{0,10},PlotStyle->Black],
Plot[r1/.b->bScaled/.V->Vval,{ϵ,0,1},PlotRange->{0,10},PlotStyle->Red],
Plot[r2/.b->bScaled/.V->Vval,{ϵ,0,1},PlotRange->{0,10},PlotStyle->Blue]]

Only r2 (the tangent solution) gives an answer. There is only one eigenvalue, and it is close to the top of the well.

FindRoot[rScaled-r2/.b->bScaled/.V->Vval,{ϵ,0.8}]

In fact, the eigenvalue only gets closer to the top of the well for smaller values of V

FindRoot[rScaled-r2/.b->bScaled/.V->0.2,{ϵ,0.8}]

For larger values of V there will be several eigenvalues

Vval=50;
Show[
Plot[rScaled,{ϵ,0,1},PlotRange->{0,10},PlotStyle->Black],
Plot[r1/.b->bScaled/.V->Vval,{ϵ,0,1},PlotRange->{0,10},PlotStyle->Red],
Plot[r2/.b->bScaled/.V->Vval,{ϵ,0,1},PlotRange->{0,10},PlotStyle->Blue]]

This time we find three eigenvalues

FindRoot[rScaled-r2/.b->bScaled/.V->Vval,{ϵ,0.1}]
FindRoot[rScaled-r1/.b->bScaled/.V->Vval,{ϵ,0.5}]
FindRoot[rScaled-r2/.b->bScaled/.V->Vval,{ϵ,0.9}]

Plot a wave function

First go back to V=1

Vval=1;
ϵV1=FindRoot[rScaled-r2/.b->bScaled/.V->Vval,{ϵ,0.8}];

Get the values of k and κ

kVal=(Sqrt[ϵ V]/a)/.ϵV1/.V->Vval;
κVal=(Sqrt[(1-ϵ) V]/a)/.ϵV1/.V->Vval;

Solve our boundary matching equations for A3, B1, and B2 in terms of A1=A

solA3=Solve[eq1/.A1->A,A3][[1]];
solB1=Solve[eq2/.solA3,B1][[1]];
solB2=Solve[eq3/.solB1/.A1->A,B2][[1]];

Now build the wave function pieces

ψ1V1=ψ1/.A1->A/.solA3/.solB1/.solB2/.{k->kVal,κ->κVal};
ψ2V1=ψ2/.A1->A/.solA3/.solB1/.solB2/.{k->kVal,κ->κVal};
ψ3V1=ψ3/.A1->A/.solA3/.solB1/.solB2/.{k->kVal,κ->κVal};

Finally, plot the wave function

aVal=3;
aMin=-3aVal;
aMax=4aVal;
A=1
Show[
Plot[ψ1V1/.a->aVal,{x,0,aVal},PlotRange->{{aMin,aMax},{0,2}}],
Plot[ψ2V1/.a->aVal,{x,aVal,aMax}],
Plot[ψ3V1/.a->aVal,{x,aMin,0}]]

The wave function has very long tails, because the energy eigenvalue is close to the top of the well

Plot[{pe/.{a->aVal,V->Vval},ϵ Vval/.ϵV1},{x,aMin,aMax},PlotRange->{0,Vval}]

Plot another wave function

Pick a somewhat higher well

Vval=5;
ϵV1=FindRoot[rScaled-r2/.b->bScaled/.V->Vval,{ϵ,0.8}];

kVal=(Sqrt[ϵ V]/a)/.ϵV1/.V->Vval;
κVal=(Sqrt[(1-ϵ) V]/a)/.ϵV1/.V->Vval;

ψ1V1=ψ1/.A1->A/.solA3/.solB1/.solB2/.{k->kVal,κ->κVal};
ψ2V1=ψ2/.A1->A/.solA3/.solB1/.solB2/.{k->kVal,κ->κVal};
ψ3V1=ψ3/.A1->A/.solA3/.solB1/.solB2/.{k->kVal,κ->κVal};

aVal=3;
aMin=-3aVal;
aMax=4aVal;
A=1
Show[
Plot[ψ1V1/.a->aVal,{x,0,aVal},PlotRange->{{aMin,aMax},{0,2}}],
Plot[ψ2V1/.a->aVal,{x,aVal,aMax}],
Plot[ψ3V1/.a->aVal,{x,aMin,0}]]

The tails are not so long, because the wave function is not so close to the top of the well

Plot[{pe/.{a->aVal,V->Vval},ϵ Vval/.ϵV1},{x,aMin,aMax},PlotRange->{0,Vval}]

Plot a set of wave functions

Return to the example of V=50. We found three eigenvalues

Vval=50;
ϵV1=FindRoot[rScaled-r2/.b->bScaled/.V->Vval,{ϵ,0.1}];
ϵV2=FindRoot[rScaled-r1/.b->bScaled/.V->Vval,{ϵ,0.5}];
ϵV3=FindRoot[rScaled-r2/.b->bScaled/.V->Vval,{ϵ,0.9}];

kVal=(Sqrt[ϵ V]/a)/.ϵV1/.V->Vval;
κVal=(Sqrt[(1-ϵ) V]/a)/.ϵV1/.V->Vval;
ψ1V1=ψ1/.A1->A/.solA3/.solB1/.solB2/.{k->kVal,κ->κVal};
ψ2V1=ψ2/.A1->A/.solA3/.solB1/.solB2/.{k->kVal,κ->κVal};
ψ3V1=ψ3/.A1->A/.solA3/.solB1/.solB2/.{k->kVal,κ->κVal};

kVal=(Sqrt[ϵ V]/a)/.ϵV2/.V->Vval;
κVal=(Sqrt[(1-ϵ) V]/a)/.ϵV2/.V->Vval;
ψ1V2=ψ1/.A1->A/.solA3/.solB1/.solB2/.{k->kVal,κ->κVal};
ψ2V2=ψ2/.A1->A/.solA3/.solB1/.solB2/.{k->kVal,κ->κVal};
ψ3V2=ψ3/.A1->A/.solA3/.solB1/.solB2/.{k->kVal,κ->κVal};

kVal=(Sqrt[ϵ V]/a)/.ϵV3/.V->Vval;
κVal=(Sqrt[(1-ϵ) V]/a)/.ϵV3/.V->Vval;
ψ1V3=ψ1/.A1->A/.solA3/.solB1/.solB2/.{k->kVal,κ->κVal};
ψ2V3=ψ2/.A1->A/.solA3/.solB1/.solB2/.{k->kVal,κ->κVal};
ψ3V3=ψ3/.A1->A/.solA3/.solB1/.solB2/.{k->kVal,κ->κVal};

aVal=3;
aMin=-aVal;
aMax=2aVal;
A=1;
Show[
Plot[ψ1V1/.a->aVal,{x,0,aVal},PlotRange->{{aMin,aMax},{-3,3}},PlotStyle->Black],
Plot[ψ2V1/.a->aVal,{x,aVal,aMax},PlotStyle->Black],
Plot[ψ3V1/.a->aVal,{x,aMin,0},PlotStyle->Black],
Plot[ψ1V2/.a->aVal,{x,0,aVal},PlotRange->{{aMin,aMax},{-3,3}},PlotStyle->Red],
Plot[ψ2V2/.a->aVal,{x,aVal,aMax},PlotStyle->Red],
Plot[ψ3V2/.a->aVal,{x,aMin,0},PlotStyle->Red],
Plot[ψ1V3/.a->aVal,{x,0,aVal},PlotRange->{{aMin,aMax},{-3,3}},PlotStyle->Blue],
Plot[ψ2V3/.a->aVal,{x,aVal,aMax},PlotStyle->Blue],
Plot[ψ3V3/.a->aVal,{x,aMin,0},PlotStyle->Blue]]

The wave functions have alternating “parity”, and the tails get longer as the energy level gets closer to the top of the well

Plot[{pe/.{a->aVal,V->Vval},ϵ Vval/.ϵV1,ϵ Vval/.ϵV2,ϵ Vval/.ϵV3},
{x,aMin,aMax},PlotRange->{0,Vval}]

Putting in some numbers

Let’s look at an electron in a box the size of an atom. We’ll use the well height and the three eigenvalues we found on the last slide to calculate the “binding energy” of the electrons for each of those three states. This is the amount of energy needed to lift the electron out of the well.

be1=Vval(1-ϵ/.ϵV1)ℏ^2/(2m a^2);
be2=Vval(1-ϵ/.ϵV2)ℏ^2/(2m a^2);
be3=Vval(1-ϵ/.ϵV3)ℏ^2/(2m a^2);

Get the values for the various constants

ℏ=QuantityMagnitude[UnitConvert[Quantity["ReducedPlanckConstant"]]];
m=QuantityMagnitude[UnitConvert[Quantity["ElectronMass"]]];
a=10^(-10);
e=QuantityMagnitude[UnitConvert[Quantity["ElementaryCharge"]]];

Print out the binding energies in electron volts

be1/e
be2/e
be3/e

These are not atypical for an atom, with inner electrons bound by hundreds of eV, but with an ionization potential in the neighborhood of 10eV. Of course, it really just comes down to the value of (ℏ^2)/2(ma^2) which you can get from simple scaling relations, but at least it is good to see the numbers work out to something reasonable.