Double infinite square well, energies

[Clever-Name]

Homework Statement

We're modelling an ammonia maser with a double infinite square well defined by:
<br /> V(x) = <br /> \begin{cases}<br /> V_{0} &amp; |x| &lt; b - \frac{a}{2}\\<br /> 0 &amp; b-\frac{a}{2} &lt; |x| &lt; b+\frac{a}{2}\\<br /> \infty &amp; |x| \geq b + \frac{a}{2}<br /> \end{cases}<br />

I have had no trouble with the assignment up until I'm asked to solve for the transcendental equations for energy and solve for the lowest odd and even state using MATLAB

Homework Equations

We are told to use the following form for the wavefunction as defined below:

Even states:
<br /> \Psi(x)_{even} =<br /> \begin{cases}<br /> Acosh(\kappa x) &amp; 0 &lt; x &lt; b-\frac{a}{2}\\<br /> Bsin(k(b + \frac{a}{2} - x)) &amp; b-\frac{a}{2} &lt; x &lt; b+\frac{a}{2}<br /> \end{cases}<br />
and for the Odd states:
<br /> \Psi(x)_{odd} =<br /> \begin{cases}<br /> Asinh(\kappa x) &amp; 0 &lt; x &lt; b-\frac{a}{2}\\<br /> Bsin(k(b + \frac{a}{2} - x)) &amp; b-\frac{a}{2} &lt; x &lt; b+\frac{a}{2}<br /> \end{cases}<br />

The Attempt at a Solution

Since the well is symmetric we can ignore the left-side classically allowed region, the function will just be defined as plus or minus Bsin(k(b + \frac{a}{2} - x)) depending on whether we're looking at an even or odd state.

Now, applying the boundary conditions at x = b-a/2 gives me the following equations:
edit - Found one mistake, negative signs, but that still doesn't fix my 2nd problem.
Even:
<br /> \sqrt{\frac{V_{0} - E}{E}}tan(ka) = -coth(\kappa(b-\frac{a}{2})<br />

Odd:

<br /> \sqrt{\frac{V_{0} - E}{E}}tan(ka) = -tanh(\kappa(b-\frac{a}{2}))<br />

with

<br /> k = \sqrt{\frac{2mE}{\hbar^2}}<br />
and
<br /> \kappa = \sqrt{\frac{2m(V_{0}-E)}{\hbar^2}}<br />

and I should note that V_{0} &lt; E

Typing out all of my work would have been a huge pain, if these don't look right let me know and I can type out what I did.

Now, my first question is: Do these equations for the energy look right?

If they are right my second question will be how the heck do I solve for the energies?! I'm using the fzero function, with a given initial guess of ~0.00661eV, in MATLAB and am getting ~0.041eV for the even state but 3.7e-16 for the odd state. They should be almost identical!

If anyone wants to tinker around with it and try and test my answers we are given:
V_{0} = 0.130 eV
b = 0.070 nm
a = 0.112 nm

and to get the energy guess we are told
k \approx \frac{\pi}{a}

 

[vela]

Your equations look fine, but you want V₀ > E so that κ is real. In terms of k and κ, your equations are
$$k \cot ka = \begin{cases}
-\kappa\tanh \kappa(b-a/2) & \text{even} \\
-\kappa\coth \kappa(b-a/2) & \text{odd}
\end{cases}$$and $\kappa^2 = k_v^2 - k^2$ where $k_v^2 = 2mV_0/\hbar^2$. Try plotting both sides of the equation as a function of k and see where the graphs intersect. That'll give you a better first guess for k. When I used $k=\pi/a$ for the initial guess, Mathematica had trouble converging on the solution.

I used $m = 1.304 \times 10^{10}~\text{eV/c}^2$ (the mass of a nitrogen atom), and with your numbers, I found the energy of the states to be both around $1.103 \times 10^{-3}\text{ eV}$, which is much smaller than V₀. The difference in the energies was $6.79 \times 10^{-8}\text{ eV}$.

 

[Clever-Name]

Ah yes, I meant to say V > E, typed it in wrong in my initial post and now I can't edit it for some reason.

We are asked to use a reduced mass of
<br /> m = \frac{(3m_{H})*m_{N}}{3m_{H} + m_{N}}<br />

and so that yields m = 2.31634 \times 10^{9} eV/c^{2}

Using this with the approximate value for k though and plugging into:

<br /> <br /> E = \frac{k^{2}\hbar^{2}}{2m}<br /> <br />

I'm getting a value for the energy on the order of 10^{-20} eV

That can't be right!

Even with plotting and finding k values then converting to energy I'm still in that range. What am I doing wrong?

Edit -

Here's the code I'm using to solve for the zeros in MATLAB:

function energies = even(E)
m = 2.316346*10^9;
hbar = 6.58211928*10^-16; 
a = 0.122*10^-9;
b = 0.070*10^-9;
V = 0.130;
alpha = sqrt((2*m*V)/(hbar^2));
k = sqrt((2*m*E)/(hbar^2));
K = sqrt(alpha^2 - k^2);

energies = k*cot(k*a) + K*tanh(K*(b - (a/2)));

fzero(@(E) even(E), E_guess)

 

[vela]

Tip: A useful combination to know is $\hbar c = 197~\text{eV nm}$.

If $k = \pi/a \approx 28\text{ nm}^{-1}$, you get
$$ E = \frac{(\hbar c)^2 k^2}{2 mc^2} = \frac{(197\text{ eV nm})^2 (28\text{ nm}^{-1})^2}{2(2.31634 \times 10^{9}\text{ eV})} = 0.00659\text{ eV}$$

EDIT: Ah, I just saw your edit. You're just having a problem with units, I think — the speed of light in particular. Using that value for the mass implicitly assumes c=1. You can probably fix this by changing your value for $\hbar$ to 197.

 

[Clever-Name]

Ah, of course. Finally! It was the darn factors of c^2 that was screwing me up.

Now I finally get a solution of

<br /> E_{even} = 0.005573391601792 eV<br />

<br /> E_{odd} = 0.005573391601791 eV<br />

My difference is in the 10^-{15} eV range, not quite the same result as you were getting, not sure why either.

 

[vela]

I recalculated the energies with the reduced mass and found
\begin{align*}
E_\text{even} &= 5.638\times 10^{-3}\text{ eV} \\
E_\text{odd} &= 5.782\times 10^{-3}\text{ eV} \\
\Delta E &= 0.099\times 10^{-3}\text{ eV}
\end{align*}

 

[Clever-Name]

What are you using as your initial energy guess? Because after changing my code to:

m = 2.316346*10^9; %eV/c^2
h = 197.326972; %eV nm    hbar*c
a = 0.122; %nm
b = 0.070; %nm
V = 0.130; %eV
c = 299792458; %m/s

alpha = sqrt((2*m*V)/((h)^2));
k = sqrt((2*m*E)/((h)^2));
K = sqrt(alpha^2 - k^2);

energies = k*cot(k*a) + K*tanh(K*(b - (a/2)));

and odd is the same except the tanh is a coth

I am still getting the values in my last post

 

[vela]

I used k=26, which corresponds to E=5.663x10^-3 eV.

I attached a plot of k cot ka and -κ tanh κ(b-a/2). The graphs intersect at k=26.04 nm^-1.

Oh, I see the difference. In your OP, you said a=0.112 nm. In your code, you have a=0.122 nm. Which one is correct? (I recalculated with a=0.122 nm, and my results still don't match yours.)

 

[Clever-Name]

Ah, yes 0.112 is right, not 0.122, just a typo on my end.

Since the prof explicitly tells use to use k = pi/a as our guess when evaluating this I am still using that in my initial guess (E \approx 0.006613 eV. With the change in a in my code I get, GAH, the exact same value

E_{even} = E_{odd} = 0.006613070853476 \text{ eV}

With your initial guess I get:

<br /> E_{even} = 0.005699839962893 \text{ eV}<br />
<br /> E_{odd} = 0.005799847037600 \text{ eV}<br />
<br /> {\Delta}E \approx 1 \times 10^{-4} \text{ eV}<br />

Which is much closer to what you ended up getting initially and is a little more dramatic.

I'm not sure what to do now :S

 

[vela]

I added a few more digits onto my value for $\hbar c$, and now my results match yours.

 

[Clever-Name]

Ok that's good enough for me then. I think the issue with taking k ~ pi/a is that it occurs right around the asymptote just to the right of the plot you posted. I used the better k value and explained why in my assignment. Man this question was a headache. Thanks a lot for your help!