Quantum Wells -- An Eigenvalue Problem

Introduction

This notebook follows the lines of the previous notebook on the rectangular well. In units where ℏ = m = 1 (where m is the mass of the particle), Schrodinger's equation is

(((d^2)ψ)/(d (x^2))) + 2 (E - V(x)) ψ = 0

We will find the lowest energy level and corresponding wavefunction for each parity. We take the potential

V(x) = -3(sech^2)(x)

because, as shown in some of the more advanced texts, one can solve exactly for the bound states. 

The potential is an even function of x and so the eigenstates are either even or odd functions of x, ψ(x) = ±ψ(-x). The sign, +1 or -1, is called the parity of the state. It is shown in the textbooks on quantum mechanics  that (i)  the ground state is symmetric ψ(x) = ψ(-x) and has no zeroes (nodes), (ii) the first excited state is odd and has one node, and (iii) for each higher energy eigenvalue the the number of nodes increases by oneand the parity alternates. 

The method we will use requires that V(x) = 0 for |x| > L/2. Now our V(x) is never quite zero but it is very close to zero for |x| > 5. Hence we take L=10.  First we define and plot the potential:

Clear ["Global`*"]

L = 10;

v[x_] :=-3 Sech[x ]^2 /; Abs[x] <= L/2

v[x_] := 0 /; Abs[x] > L/2

Plot[{v[x], -3.0 + 100 * (x - 5), -3.0 + 100 * (x + 5)},{x,-0.75 L,0.75L},PlotStyle ->{ {Red, AbsoluteThickness[3]}, {Black, Dashing[{0.01,0.03}], AbsoluteThickness[3]}, {Black,Dashing[{0.01,0.03}], AbsoluteThickness[3]}},AxesLabel->{"x","V(x)"}, PlotRange->{-3.1,0.1},Epilog->{Text[Style["Region 1",FontSize->14],{-6.5,-1.2}],Text[Style["Region 3",FontSize->14],{6.5,-1.2}],Text[Style["Region 2",FontSize->14],{0,-1.2}]}]

Setting up the Problem

In region 1,  (Subscript[ψ, 1])= A (e^(κx))where κ = (Sqrt[2|E|]),and in region 3, (Subscript[ψ, 3])= B (e^(-κx)), so we set up these functions:

B = 1;

ψ1[x_] := A Exp[ Abs[2 en]^(1/2) x]

ψ3[x_] := B Exp[- Abs[2 en]^(1/2) x]

We also define the Schrodinger equation for region 2, using a delayed assignment, ":=", since we will only use it later:

eqn[en_] := ψ2''[x] + 2(en - v[x])ψ2[x]

We also set up the calculation of the wavefunction in region 2 matching the function and its derivative to (Subscript[ψ, 1])at x=-L/2. We will adjust the energy so that either the derivative of ψ vanishes at x=0 (for even eigenfunctions) or ψ vanishes (for odd eigenfunctions). This value for the energy will be an eigenvalue. Since the energy will be determined by a boundary condition at x = 0 we only need to integrate from x = -L/2 up to x = 0.

wavefunc2[energy_] := (en = energy; NDSolve[{ eqn[energy] == 0, ψ2[-L/2] == ψ1[-L/2], 
					ψ2'[-L/2] == ψ1'[-L/2]}, 
		ψ2, {x, -L/2, 0} ])

"wavefunc2[en_]"  is given as a replacement rule in the form "{{ψ2->InterpolatingFunction[{{-0.5,0.5}},"<>"]}}". In order to directly access the wavefunction in region 2 we define a function sol2[x, en], which applies the replacement rule, removing one of the sets of curly brackets, i.e.

sol2[x_?NumericQ, en_?NumericQ] := ψ2[x] /. wavefunc2[en][[1]]

We also do the same for the derivative of the wavefunction

sol2prime[x_?NumericQ,en_?NumericQ]:=ψ2'[x]/.wavefunc2[en][[1]]

Even Parity Solution

We are now in a position to start calculating. First of all we look for an even eigenfuction by making the amplitudes of the wavefunction in regions 1 and 3 equal:

A= B;

We require that the derivative of the wavefunction vanishes at x = 0. We give two initial starting guesses

evalue=energy /. FindRoot[sol2prime[0,energy],{energy,-5, -1}]

This result agrees with the exact value which is -2. This state is actually the only even parity bound state.

Now we want the eigenfunction coresponding to our eigenvalue. Since we now have the eigenvalue, we do not  want to keep recalculating  the wavefunction so we define a function "efunc2" with immediate assignment, where we input the eigenvalue for the energy

efunc2[x_]=ψ2[x]/.wavefunc2[evalue][[1]];

We have now obtained the wavefunction in all three regions, so let's collect these into a single function ψ [x_] , which can then easily be plotted. Remember that we only computed efunc2 [x] for -L/2 ≤ x < 0, and so we also have to specify it for 0 ≤  x ≤ L/2, noting that the eigenfunction is even):

ψ [x_] := efunc2 [x] /;-L/2 ≤ x ≤ 0

ψ [x_] := efunc2 [-x] /;0 ≤ x ≤ L/2

ψ [x_] := ψ1[x] /; x < -L/2

ψ [x_] := ψ3[x]/; x > L/2

We first normalize the wavefunction:

normconst = Sqrt[ NIntegrate [ ψ[x]^2, {x, -Infinity, Infinity} ] ];

ψnorm[x_] := ψ[x] / normconst

and then plot it. We first define the plot for ψnorm (but do not yet display it because of thesemicolon after the comand). Then we add the labels with the "Show[Graphics[... ]]" command. This looks a bit complicated but has the advantage that the current values of the well depth and width, and the eigenvalue are automatically printed without needing to modify this command each time.

fig = Plot[ψnorm[x], {x, -0.55L, 0.55L}, AxesLabel -> {"x", "ψ"}, PlotStyle -> {Red, AbsoluteThickness[2]} ];

figsech = Plot[Sqrt[3/4] Sech[x]^2, {x, -0.55L, 0.55L}, PlotRange -> {0, 1},PlotStyle -> {Blue, Dashing[{0.01, 0.03}], AbsoluteThickness[6]} ];

Show[fig,figsech, Graphics [ { 
			Text[ evalue , {3, 0.7}, {-1, 0} ] ,
Text[ "E = ", {3, 0.7}, {1, 0}  ],
			Text[ "L = 10", {3.5, 0.6}, {1, 0}] } ]]

The curve tracks precisely the exact solution,  ψ(x) = (Sqrt[3/4])(sech^2)(x), shown by the thick dashed line. The absence of nodes (zeroes) in the wavefunction confirms, since we are in one dimension, that it is the ground state.  You can easily verify by hand that ψ(x) is a normalized eigenstate of the Hamiltonian with eigenvalue -2.  Mathematica can also do this:

en = -2; ψexact[x_] = Sqrt[3/4] Sech[x]^2; pot[x_] = -3 Sech[x]^2;

ψexact''[x] + 2(en - pot[x])ψexact[x] == 0 //FullSimplify

Integrate[ψexact[x]^2, {x, -Infinity, Infinity}]

Odd Parity Solution

We now look for an odd parity solution.

A= -B;

We use the "FindRoot" command to locate the eigenvalue and give it two starting values, which we take to be -0.2 and -2. The boundary condition is that the wavefunction vanishes at the origin

evalue= en/.FindRoot[sol2[0,en],{en,-0.2, -2}]

This agrees with the exact value of -1/2. This is actually the only odd parity bound state.

Now we want the eigenfunction coresponding to our eigenvalue.We define a function "efunc2" with immediate assignment, where we input the eigenvalue for the energy:

efunc2[x_]=ψ2[x]/.wavefunc2[evalue][[1]];

We have now obtained the wavefunction in all three regions, so let's collect these into a single function ψ [x_] , which can then easily be plotted

ψ [x_] := -efunc2 [-x] /;0 ≤ x ≤ L/2

We first normalize the wavefunction:

normconst = Sqrt[ NIntegrate [ ψ[x]^2, {x, -Infinity, Infinity} ] ];

and then plot it.

fig = Plot[Evaluate[ψnorm[x]], {x, -0.7L, 0.7L},PlotStyle -> {Red, AbsoluteThickness[2]}, AxesLabel -> {"x", "ψ"}];

figsech = Plot[Sqrt[3/2] Sinh[x] Sech[x]^2, {x, -0.7L, 0.7L}, PlotStyle -> {Blue, Dashing[{0.01, 0.03}], AbsoluteThickness[4]} ];

Show[fig,figsech,  Graphics [ { 
			Text[evalue, {3.5, 0.5}, {-1, 0} ] ,
Text[ "E = ", {3.5, 0.5}, {1, 0}  ],
			Text[ "L = 10", {4, 0.6}, {1, 0}] } ]]

The curve tracks precisely the exact solution, ψ(x) = (Sqrt[3/2])sinh(x) (sech^2)(x), show by the thick dashed line. There is one node which, since we are in one dimension, shows that this is the 1st excited state.

Again Mathematica confirms that this is the exact, normalized, solution:

en = -1/2; ψexact[x_] = Sqrt[3/2] Sinh[x] Sech[x]^2;

ψexact''[x] + 2(en - pot[x])ψexact[x] == 0 //FullSimplify

Integrate[ψexact[x]^2, {x, -Infinity, Infinity}]