Hi.I am trying to find the energy levels in a finite potential well through a numerical method. The only numerical method I have found is performing iterations on the odd and even parity equations (derived from the schrodinger equation) using the newton-raphson method. Is there a better method or is this a good method?

Thanks

 

Hi.I am trying to find the energy levels in a finite potential well through a numerical method. The only numerical method I have found is performing iterations on the odd and even parity equations (derived from the schrodinger equation) using the newton-raphson method. Is there a better method or is this a good method?

Thanks

Can you give details about the potential function, such as number of dimensions (1D, 2D or 3D) and functional shape.

 

Ok. It is one a one dimensional symmetric well spanning from -L to L with the well centred at 0. Potential in the well is zero and potential of barriers is V.

Broke the wavefunction into three pieces
ψ(x)= ψ1(x), x< -L
ψ2(x), -L<= x <= L
ψ3(x), x> L

By using schrodinger equation to solve for the potential inside and outside I end up with
even parity at energies that satisfy

√((V-E)/E) = tan(√((2mE)/h)L)

and odd parity at √((V-E)/E) = - cot(√((2mE)/h)L).

 

Ok. It is one a one dimensional symmetric well spanning from -L to L with the well centred at 0. Potential in the well is zero and potential of barriers is V.

Broke the wavefunction into three pieces
ψ(x)= ψ1(x), x< -L
ψ2(x), -L<= x <= L
ψ3(x), x> L

By using schrodinger equation to solve for the potential inside and outside I end up with
even parity at energies that satisfy

√((V-E)/E) = tan(√((2mE)/h)L)

and odd parity at √((V-E)/E) = - cot(√((2mE)/h)L).

OK, so it sounds like this is a case where you know the solution ahead of time, but you want to solve it using numerical techniques as a learning exercise. Am I correct about that?

You already know the energy levels and the wave function solutions are known as sinusoidal and exponential functions.

In this particular case you can use a simple technique of matching boundary conditions to find the energy levels. This is one I once used for waveguide solutions, and the two problems are essentially identical from a mathematical point of view.

Basically, the bound state solutions will be a sine or cosine function inside the well, and an exponential function outside the well. Bound states will have decaying exponentials. The two functions (inside vs. outside) must match in value and also in slope at the boundary. You can write equations that force the solutions to match in value, and then scan the range of energy levels to see where the derivative (slopes) also match. You can do a course scan of the range to identify where zero crossings occur. This will tell you the number of bound state modes. Then you can do a finer binary search near each solution.

I've solved many waveguide problems using many different numerical techniques, and I did this in the days when computers were pretty slow by today's standards. The above technique is the most efficient numerical technique for such problems. It can even be used when the potential function is not a simple step-like shape. Continuously changing shapes can be approximated by a staircase function with many narrow gaps, each with constant energy potential. Every step has sinusoidal or exponential solutions inside. The boundaries at each step can be matched. All boundaries are matched in value and slope, with the exception of one boundary which is matched only in value. Then the energy level range can be scanned to match the slope on that last remaining boundary.

 

Hi. Thanks for your help. I understand the mathematics behind the process and yes the solution you mentioned is the one I have come across but I need to model this and I guess it is not apparent to me as to how exactly i will obtain values for E.
I require the electronvolt values for further calculations.
The only numerical method I have come across says to perform iterations on those equations I stated earlier until I obtain valid values for the the different energies.
The link to what I have been looking at is below
http://en.wikibooks.org/wiki/Materials_in_Electronics/Confined_Particles/1D_Finite_Wells
if you scroll to the "numerical solution".

 

... it is not apparent to me as to how exactly i will obtain values for E ...

Doing this takes a little mathematical manipulation, but it is fairly simple in principle.

The exponential function (outside the well) and the sinusoidal function (inside the well) are functions of distance times a spatial frequency number. This spacial frequency is a function of the energy level value. Hence, the energy eigenvalue controls the shapes of those functions.

Graphically, the wavefunction must be smooth looking, which means that mathematically it must be continuous and differentiable across the boundary. This is only possible at discrete energy levels. So you can easily write equations that force the two pieces of the wavefunction to be equal at the boundary. Then, you search for values of energy that also make the slope of the curves match at the boundary.

The searching method is your choice. You and that article both mention the Newton-Raphson method. This is fine, but I always perferred computer search algorithms. You simple scan the range of energy levels using small steps, and you will be able to see when the difference in slopes (on each side of the boundary) changes sign. This tells you that you just went past one of the eigenvalues. After you do this course search, you then know how many bound state modes are allowed by the structure. Note that there will be a finite number, unless you have an infinite well.

The next step is to use a binary search to find the precise values of the energy levels. Binary searches are very fast and you get double precision accuracy in about 50 iterations. A binary search is possible because you can always tell whether your guess is too high or too low based on the sign of the slope difference.

 

Hi. Thanks for your help.I'm trying to get to grips with this and I think it is starting to make sense. So if lets say I plug in the respective equations for the sinusoidal wavefunction in the well and the exponential wavefunction outside the well I would then have to obtain values from both equations through inputting values of L (distance), and when continually subtracting the value of the exponential from the sinusoid a change in sign will indicate that there is an eigen value there. Is that correct?

From there I know the points at which aenergy levels exists and can carry on to use a binary search to find the exact values?

Are these steps correct?

 

Are these steps correct?

It's hard to be sure if this is correct based on your descirption. I would recommend that you start writing out the equations and posting them. It will then be easier for me and others to make comment and guide you. It sounds like you are on the right track.

 

Cool. Ok so for the bound state where E<V the particle is bound in the well the wavefunctions for the three regions ψ1(x), ψ2(x) and ψ3(x) are

ψ1(x)= Ge^(αx)
ψ2(x)= Asin(kx) + Bcos(kx)
ψ3(x)= He^(-αx)

where a=√(2m(V-E)/h) and k=√((2mE)/h) h- reduced plancks constant

For even wavefunctions there must be no sine components

ψ1(x)= Ge^(αx)
ψ2(x)= Bcos(kx)
ψ3(x)= He^(-αx)

At the boundary the condition for coninuity is

Ge^(-αL)=Bcos(-kL) for the left side of the well


and for continuity of derivatives gives

αGe^(-αL) = -kBsin(-kL) (1)


and so if I evaluate the two sides of equation (1) for various values of L and continually subtract one side from the other then when a sign change occurs it indicates that an eigen value has just been passed. Is this correct up till here?

 

and so if I evaluate the two sides of equation (1) for various values of L and continually subtract one side from the other then when a sign change occurs it indicates that an eigen value has just been passed. Is this correct up till here?

Everything was correct, with the exception for this part. The value of L is fixed by the problem. Instead you scan the possible values of E. Perhaps this was just a typo, but I want to make sure you understand that.

 

Ok. So the way I picture it is that by doing this we are pivoted at distance L and we are moving upwards along the boundary. The sinusoidal waveforms only exist at specific points along the boundary whereas the exponential components exist all the way along the boundary.
So I perform the scanning using a range of E values such that I will get something like the following assuming these are values obtained as we increase E (moving up along the boundary) from E1 to E6. With the e values being the value from the exponential and s from the sinusoidal wavefunction.


(Energy) (exponential result) (sinusoidal result)

E1 e1 0
E2 e2 0
E3 e3 0
E4 e4 s4
E5 e5 0
E6 e6 0

And then I can do the e values minus the s values. That set of data will then let me know where the energy levels exist. Yes?

Sorry for dragging this, your help has been immense.

 

Ok. So the way I picture it is that by doing this we are pivoted at distance L and we are moving upwards along the boundary. The sinusoidal waveforms only exist at specific points along the boundary whereas the exponential components exist all the way along the boundary.
So I perform the scanning using a range of E values such that I will get something like the following assuming these are values obtained as we increase E (moving up along the boundary) from E1 to E6. With the e values being the value from the exponential and s from the sinusoidal wavefunction.


(Energy) (exponential result) (sinusoidal result)

E1 e1 0
E2 e2 0
E3 e3 0
E4 e4 s4
E5 e5 0
E6 e6 0

And then I can do the e values minus the s values. That set of data will then let me know where the energy levels exist. Yes?

Sorry for dragging this, your help has been immense.

I'm a little confused your explanation, but it sounds like you are not quite visualizing it correctly.

You can match the value of the exponential and the value of the sinusoidal function for any value of energy. So if you try to graph the function, you will see a continuous function, but the slopes will look disjointed (or kinked) at the boundary most of the time. This lack of smoothness is not allowed. So, you need to find the discrete energy levels that allow both the sinusoidal function and the exponential function to have matching slopes.

The clue to knowing where the solutions of E are, is to look at the kink at the boundary. You will see a kink that faces upward, suddenly become a kink that faces downward. This occurs when you skip over the correct E that would have removed the kink entirely. The testing of slopes is done numerically by subtracting one slope from the other slope. Then your solutions will occur when the difference in slope is equal to zero. Your solutions then become more obvious because you will see zero crossings and sign changes as you scan the range of E.

 

Ok. So if I take a set of E values within a realistic range and with fairly good precision and plug the values in both sides of the equation and do the subtraction then I will get results with the signs as follows

- - - - - + + + + + - - - - - - + + + + + + - - - - - - -

and so I need to investigate the location of the changes in sign. Correct?

Then I need to perform more in depth scans centred around the rough E values located in order to get good precision?

 

Ok. So if I take a set of E values within a realistic range and with fairly good precision and plug the values in both sides of the equation and do the subtraction then I will get results with the signs as follows

- - - - - + + + + + - - - - - - + + + + + + - - - - - - -

and so I need to investigate the location of the changes in sign. Correct?

Then I need to perform more in depth scans centred around the rough E values located in order to get good precision?

Yes, just make sure that your initial scan is fine enough. Once you have bounds with a single solution in between, then use a binary search to converge to the final values. It's a very primitive technique.

 

Hi. I guess understanding the procedure doesn't mean I know how to put it into code.
What I was thinking was to do
1) Eq1 - Eq2 as we discussed earlier for various values of V and then create an array holding the results.
2)Scan the results array to find where there are sign changes and backtrack the index to determine at what energy that occured
3) Do a binary search to refine energy value based on the value obtained in step 2

If all this is correct can you explain to me the binary search?

 

Hi. I apologise for having generated several threads on this I was unable to locate the previous but have found it now.
So I have written a matlab code which takes the difference between the derivatives of the wavefunction inside the well and outside the well. I then check for sign changes and save the the two energy values which bracket the root in order to use them later in a finer search to get a precise value for each energy level.
If this procedure is right then it means the function for derivative of the wavefunction I am using for even parity may be incorrect.

I have the following wavefunction

(outside well) a/k = tan(kL) (inside well) (1)

a=√(2m(V-E)/hbar) and k=√((2mE)/hbar) and energy levels exist when the E values satisfy the equation (1)

the derivatives of eq (1) I am using are shown below

((-V)/((2*E*E*sqrt((V-E)/E)))) (Outside the well) and then

(inside well)
(((L*me)/(hbar*sqrt(2*me*E)))*(sec2((L*(sqrt(2*me*E)))/hbar))

I evaluate the difference of the derivatives for a range of E but all result in the same sign so I don't know where I am going wrong.

Matlab code below

%This program calculates the energy levels in a fintite potential
%quantum well.

hbar=6.582e-016; %h/2pi [eV.s]
me=9.11e-31; %electron effective mass [Kg]
V=75; %Barrier potential in eV
L=0.35*10^-9; %Well width
i=1; %Initialising i
j=2; %Initialising j
E=0.001; %initialising E
dE=0.1; %Increment E by this value


%k=((sqrt(2*me*E))/(hbar));
%a=((sqrt((2*me*(V-E))))/(hbar));
%Evenfun=(a/k)-tan(k*L) %The even parity function


while ((V)>=(E)) %'While' loop to continue until the new energy level reaches the
%the barrier energy level



Ef0=((-V)/((2*E*E*sqrt((V-E)/E))))-(((L*me)/(hbar*sqrt(2*me*E)))*(sec((L*(sqrt(2*me*E)))/hbar))*(sec((L*(sqrt(2*me*E)))/hbar))) %Evaluate ederivative for E
E=E+dE; % Increment E by dE
Ef1=((-V)/((2*E*E*sqrt((V-E)/E))))-(((L*me)/(hbar*sqrt(2*me*E)))*(sec((L*(sqrt(2*me*E)))/hbar))*(sec((L*(sqrt(2*me*E)))/hbar))) % Evaluate even derivative for new value of E

if ((Ef0*Ef1)<0) %If a root is found (sign change between the 2 values)
Bracketroot(i)=E-dE %Place value of V in array
Bracketroot(j)=E %Place incremented V value in next position in array
i=i+1; %increment array position
j=j+1; %increment array position
end
end

 

If this procedure is right then it means the function for derivative of the wavefunction I am using for even parity may be incorrect.

For the even parity case, there is always at least one solution. So, you need to find out why you don't see at least one sign change. One thing to check is whether your starting value for the energy is not already past the first and only solution. Try different initial values of energy.

Another thing to check is the units. I notice you are using eV instead of J, this is OK as long as you are careful, but double check to make sure you don't have a unit mismatch somewhere.

Also, have you chosen the well energy and thickness to be physically reasonable, or are they just random choices? You want to start with usual values that won't create numerical issues.

Here is a link to an article on the problem, and this shows a different approach to solve using graphical techniques. You can create similar graphs, and you will be able to see what the issue is.

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

 

Following on from the previous threads do you know how I can go about integrating properties such as the barrier and well doping levels in order to account for how they influence energy levels?

 

The PDEs used in Finance and Physics have many similarities. To solve for option prices with complicated, but smooth, payoff functions we used the Crank-Nicolson method.

http://en.wikipedia.org/wiki/Crank%E2%80%93Nicolson_method

It may or may not be applicable in your case. I know the the Shroedinger Equation is not the same as the Heat Equation, but it might be worh a look. Let me know what you find out.

 

Well not exactly sure if its applicable. I am trying to calculate the energy levels but now take into account the impurity concentration in the barrier. Does anyone know anything about this?