# Fermions and bosons in the infinite square well

Shankar, R. (1994), Principles of Quantum Mechanics, Plenum Press. Chapter 10, Exercise 10.3.4.

Suppose we have two identical particles in an infinite square well. The energy levels in a well of width ${L}$ are

$\displaystyle E=\frac{\left(\pi n\hbar\right)^{2}}{2mL^{2}} \ \ \ \ \ (1)$

where ${n=1,2,3,\ldots}$ The corresponding wave functions are given by

$\displaystyle \psi_{n}\left(x\right)=\sqrt{\frac{2}{L}}\sin\frac{n\pi x}{L} \ \ \ \ \ (2)$

If the total energy of the two particles is ${\pi^{2}\hbar^{2}/mL^{2}}$, the only possible configuration is for both particles to be in the ground state ${n=1}$. This means the particles must be bosons, so the state vector is

$\displaystyle \left|x_{1},x_{2}\right\rangle =\frac{2}{L}\sin\frac{\pi x_{1}}{L}\sin\frac{\pi x_{2}}{L} \ \ \ \ \ (3)$

If the total energy is ${5\pi^{2}\hbar^{2}/2mL^{2}}$, then one particle is in the state ${n=1}$ and the other is in ${n=2}$. Since the states are different, the particles can be either bosons or fermions. For bosons, the state vector is

 $\displaystyle \left|x_{1},x_{2}\right\rangle$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{2}}\left[\frac{2}{L}\sin\frac{\pi x_{1}}{L}\sin\frac{2\pi x_{2}}{L}+\frac{2}{L}\sin\frac{2\pi x_{1}}{L}\sin\frac{\pi x_{2}}{L}\right]\ \ \ \ \ (4)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{\sqrt{2}}{L}\left[\sin\frac{\pi x_{1}}{L}\sin\frac{2\pi x_{2}}{L}+\sin\frac{2\pi x_{1}}{L}\sin\frac{\pi x_{2}}{L}\right] \ \ \ \ \ (5)$

For fermions, the state must be antisymmetric, so we have

$\displaystyle \left|x_{1},x_{2}\right\rangle =\frac{\sqrt{2}}{L}\left[\sin\frac{\pi x_{1}}{L}\sin\frac{2\pi x_{2}}{L}-\sin\frac{2\pi x_{1}}{L}\sin\frac{\pi x_{2}}{L}\right] \ \ \ \ \ (6)$

# Infinite square well – expanding well

References: Shankar, R. (1994), Principles of Quantum Mechanics, Plenum Press. Section 5.2, Exercise 5.2.1.

Shankar’s treatment of the infinite square well is similar to that of Griffiths, which we’ve already covered, so we won’t go through the details again. The main difference is that Shankar places the potential walls at ${x=\pm\frac{L}{2}}$ while Griffiths places them at ${x=0}$ and ${x=a}$. As a result, the stationary states found by Shankar are shifted to the left, with the result

$\displaystyle \psi_{n}\left(x\right)=\begin{cases} \sqrt{\frac{2}{L}}\cos\frac{n\pi x}{L} & n=1,3,5,7,\ldots\\ \sqrt{\frac{2}{L}}\sin\frac{n\pi x}{L} & n=2,4,6,\ldots \end{cases} \ \ \ \ \ (1)$

These results can be obtained from the form given by Griffiths (where we take the width of the well to be ${L}$ rather than ${a}$):

 $\displaystyle \psi_{n}\left(x\right)$ $\displaystyle =$ $\displaystyle \sqrt{\frac{2}{L}}\sin\frac{n\pi\left(x+\frac{L}{2}\right)}{L}\ \ \ \ \ (2)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \sqrt{\frac{2}{L}}\left[\sin\frac{n\pi x}{L}\cos\frac{n\pi}{2}+\cos\frac{n\pi x}{L}\sin\frac{n\pi}{2}\right] \ \ \ \ \ (3)$

Choosing ${n}$ to be even or odd gives the results in 1.

The specific problem we’re solving here involves a particle that starts off in the ground state (${n=1}$) of a square well of width ${L}$. The well then suddenly expands to a width of ${2L}$ symmetrically, that is, it now extends from ${x=-L}$ to ${x=+L}$. We are to find the probability that the particle will be found in the ground state of the new well.

We solved a similar problem before, but in that case the well expanded by moving its right-hand wall to the right while keeping the left-hand wall fixed, so that the particle found itself in the left half of the new, expanded well. In the present problem, the particle finds itself centred in the new expanded well. You might think that this shouldn’t matter, but it turns out to make quite a difference. To calculate this probability, we need to express the original wave function in terms of the stationary states of the expanded well, which we’ll refer to as ${\phi_{n}\left(x\right)}$. That is

$\displaystyle \psi_{1}\left(x\right)=\sum_{n=1}^{\infty}c_{n}\phi_{n}\left(x\right) \ \ \ \ \ (4)$

Working with Shankar’s functions 1 we find ${\phi_{n}}$ by replacing ${L}$ by ${2L}$:

$\displaystyle \phi_{n}\left(x\right)=\begin{cases} \frac{1}{\sqrt{L}}\cos\frac{n\pi x}{2L} & n=1,3,5,7,\ldots\\ \frac{1}{\sqrt{L}}\sin\frac{n\pi x}{2L} & n=2,4,6,\ldots \end{cases} \ \ \ \ \ (5)$

Using the orthonormality of the wave functions, we have

 $\displaystyle c_{1}$ $\displaystyle =$ $\displaystyle \int_{-L}^{L}\psi_{1}\left(x\right)\phi_{1}\left(x\right)dx\ \ \ \ \ (6)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \int_{-L/2}^{L/2}\sqrt{\frac{2}{L}}\cos\frac{\pi x}{L}\frac{1}{\sqrt{L}}\cos\frac{\pi x}{2L}dx\ \ \ \ \ (7)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{\sqrt{2}}{L}\int_{-L/2}^{L/2}\cos\frac{\pi x}{L}\cos\frac{\pi x}{2L}dx\ \ \ \ \ (8)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{\sqrt{2}}{L}\int_{-L/2}^{L/2}\left(1-2\sin^{2}\frac{\pi x}{2L}\right)\cos\frac{\pi x}{2L}dx\ \ \ \ \ (9)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{8}{3\pi} \ \ \ \ \ (10)$

The limits of integration are reduced in the second line since ${\psi_{1}\left(x\right)=0}$ if ${x>\left|\frac{L}{2}\right|}$.

Thus the probability of finding the particle in the new ground state is

$\displaystyle \left|c_{1}\right|^{2}=\frac{64}{9\pi^{2}} \ \ \ \ \ (11)$

Note that in the earlier problem where the well expanded to the right, the probability was ${\frac{32}{9\pi^{2}}}$, so the new probability is twice as much when the wave function remains centred in the new well.

We could have also done the calculation using Griffiths’s well which extended from ${x=0}$ to ${x=L}$. If this well expands symmetrically, it now runs from ${x=-\frac{L}{2}}$ to ${x=\frac{3L}{2}}$, and the stationary states of this new well are obtained by replacing ${L\rightarrow2L}$ and ${x\rightarrow x+\frac{L}{2}}$, so we have

$\displaystyle \phi_{n}\left(x\right)=\frac{1}{\sqrt{L}}\sin\frac{n\pi\left(x+\frac{L}{2}\right)}{2L} \ \ \ \ \ (12)$

We then get

 $\displaystyle c_{1}$ $\displaystyle =$ $\displaystyle \int_{-L/2}^{3L/2}\psi_{1}\left(x\right)\phi_{1}\left(x\right)dx\ \ \ \ \ (13)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{\sqrt{2}}{L}\int_{0}^{L}\sin\frac{\pi x}{L}\sin\frac{\pi\left(x+\frac{L}{2}\right)}{2L}dx\ \ \ \ \ (14)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{8}{3\pi} \ \ \ \ \ (15)$

The integral can be done by expanding the second sine using the sine addition formula. (I just used Maple.)

# Infinite square well with variable delta function barrier: location of the particle

Reference: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 10.8.

Continuing our example of the infinite square well with the growing delta function barrier, we can now find the probability of the particle being found in the right section of the well. The wave function is

$\displaystyle \psi\left(x\right)=\begin{cases} Ae^{ikx}+Be^{-ikx} & 0\le x<\frac{a}{2}+\epsilon\\ Ce^{ikx}+De^{-ikx} & \frac{a}{2}+\epsilon

With the boundary conditions we can solve for 3 of the coefficients in terms of the fourth and we got

 $\displaystyle B$ $\displaystyle =$ $\displaystyle -A\ \ \ \ \ (2)$ $\displaystyle A$ $\displaystyle =$ $\displaystyle e^{-iz}D\frac{\sin\left[\frac{1}{2}z\left(1-\delta\right)\right]}{\sin\left[\frac{1}{2}z\left(1+\delta\right)\right]}\ \ \ \ \ (3)$ $\displaystyle C$ $\displaystyle =$ $\displaystyle -De^{-2iz} \ \ \ \ \ (4)$

where ${z\equiv ka}$.

We can call the left and right wave functions ${\psi_{l}}$ and ${\psi_{r}}$ and we get

 $\displaystyle \psi_{l}$ $\displaystyle =$ $\displaystyle 2ie^{-iz}D\frac{\sin\left[\frac{1}{2}z\left(1-\delta\right)\right]}{\sin\left[\frac{1}{2}z\left(1+\delta\right)\right]}\sin kx\ \ \ \ \ (5)$ $\displaystyle \psi_{r}$ $\displaystyle =$ $\displaystyle 2ie^{-iz}D\sin\left(z-kx\right) \ \ \ \ \ (6)$

The probability of being in the right well is

$\displaystyle P_{r}=\frac{\int_{\frac{a}{2}+\epsilon}^{a}\left|\psi_{r}\right|^{2}dx}{\int_{0}^{\frac{a}{2}+\epsilon}\left|\psi_{l}\right|^{2}dx+\int_{\frac{a}{2}+\epsilon}^{a}\left|\psi_{r}\right|^{2}dx} \ \ \ \ \ (7)$

With ${\delta\equiv2\epsilon/a}$ this is

$\displaystyle P_{r}=\frac{\int_{\frac{a}{2}\left(1+\delta\right)}^{a}\left|\psi_{r}\right|^{2}dx}{\int_{0}^{\frac{a}{2}\left(1+\delta\right)}\left|\psi_{l}\right|^{2}dx+\int_{\frac{a}{2}\left(1+\delta\right)}^{a}\left|\psi_{r}\right|^{2}dx} \ \ \ \ \ (8)$

We get

 $\displaystyle X_{-}\equiv\int_{\frac{a}{2}\left(1+\delta\right)}^{a}\left|\psi_{r}\right|^{2}dx$ $\displaystyle =$ $\displaystyle \frac{D^{2}}{4k}\left[z\left(1-\delta\right)-\sin\left(z\left(1-\delta\right)\right)\right]\ \ \ \ \ (9)$ $\displaystyle X_{+}\equiv\int_{0}^{\frac{a}{2}\left(1+\delta\right)}\left|\psi_{l}\right|^{2}dx$ $\displaystyle =$ $\displaystyle \frac{D^{2}}{4k}\frac{\sin^{2}\left[\frac{1}{2}z\left(1-\delta\right)\right]}{\sin^{2}\left[\frac{1}{2}z\left(1+\delta\right)\right]}\left[\left[z\left(1+\delta\right)-2\sin\left[\frac{z}{2}\left(1+\delta\right)\right]\cos\left[\frac{z}{2}\left(1+\delta\right)\right]\right]\right]\ \ \ \ \ (10)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{D^{2}}{4k}\frac{\sin^{2}\left[\frac{1}{2}z\left(1-\delta\right)\right]}{\sin^{2}\left[\frac{1}{2}z\left(1+\delta\right)\right]}\left(z\left(1+\delta\right)-\sin\left(z\left(1+\delta\right)\right)\right) \ \ \ \ \ (11)$

So

 $\displaystyle P_{r}$ $\displaystyle =$ $\displaystyle \frac{X_{-}}{X_{+}+X_{-}}\ \ \ \ \ (12)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{1}{1+X_{+}/X_{-}} \ \ \ \ \ (13)$

with

 $\displaystyle \frac{X_{+}}{X_{-}}$ $\displaystyle =$ $\displaystyle \frac{z\left(1+\delta\right)-\sin\left(z\left(1+\delta\right)\right)}{z\left(1-\delta\right)-\sin\left(z\left(1-\delta\right)\right)}\times\frac{\sin^{2}\left[\frac{1}{2}z\left(1-\delta\right)\right]}{\sin^{2}\left[\frac{1}{2}z\left(1+\delta\right)\right]}\equiv\frac{I_{+}}{I_{-}}\ \ \ \ \ (14)$ $\displaystyle I_{+}$ $\displaystyle \equiv$ $\displaystyle \left[z\left(1+\delta\right)-\sin\left(z\left(1+\delta\right)\right)\right]\sin^{2}\left[\frac{1}{2}z\left(1-\delta\right)\right]\ \ \ \ \ (15)$ $\displaystyle I_{-}$ $\displaystyle \equiv$ $\displaystyle \left[z\left(1-\delta\right)-\sin\left(z\left(1-\delta\right)\right)\right]\sin^{2}\left[\frac{1}{2}z\left(1+\delta\right)\right] \ \ \ \ \ (16)$

We can work out this probability for several values of ${T}$, with ${\delta=0.01}$:

 ${T}$ ${z}$ ${P_{r}}$ 1 3.673 0.487 5 4.760 0.471 20 5.720 0.401 100 6.135 0.147 1000 6.215 0.00248

The probability of the particle being to the right of the barrier drops to zero as the barrier strength becomes infinite, which is just what we saw above. Because the left well has a lower ground state energy, the particle favours that region.

We can see this by plotting the wave functions 5 and 6 for these same values of ${T}$. We’ve used ${\delta=0.01}$, ${a=1}$ and omitted the common factor of ${2ie^{-iz}D}$ in the plots, so the wave functions aren’t normalized, but all we’re interested in is the relative wave functions in the two regions of the well.

We can see that the wave function drops off as ${T}$ increases.

# Infinite square well with variable delta function barrier: ground state energy

Reference: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 10.8.

Here’s another example of the adiabatic theorem. This time, we have an infinite square well in which a delta function barrier is inserted slowly at a position that is slightly off centre, so that for ${0 we have the potential

$\displaystyle V\left(x\right)=f\left(t\right)\delta\left(x-\frac{a}{2}-\epsilon\right) \ \ \ \ \ (1)$

where ${f\left(t\right)}$ is a function that rises slowly from 0 to ${\infty}$. The adiabatic theorem says that the system will remain in the ground state of the time-varying hamiltonian.

First, we’ll look at what the state is when the barrier has attained infinite strength, so that ${t\rightarrow\infty}$. [OK, the delta function itself is always infinite at a single point, but it can have a constant ‘strength’ factor multiplying it. We’ve looked at the case of the infinite square well with a constant delta function barrier and we’ve seen that increasing the strength factor to ${\infty}$ effectively divides the well into two wells that are isolated from each other, while a finite strength barrier does allow the wave function to communicate across the barrier.]

For an infinitely strong delta function barrier then, we have one well of width ${\frac{a}{2}+\epsilon}$ and one well of width ${\frac{a}{2}-\epsilon}$. The wave functions in both wells must be zero at their boundaries, so we get for the ground state (${n=1}$):

$\displaystyle \psi\left(x\right)=\begin{cases} A\sin\frac{\pi}{\frac{a}{2}+\epsilon}x & 0\le x<\frac{a}{2}+\epsilon\\ A\sin\left[\frac{\pi}{\frac{a}{2}-\epsilon}\left(x-\frac{a}{2}-\epsilon\right)\right] & \frac{a}{2}+\epsilon

 $\displaystyle E_{l}$ $\displaystyle =$ $\displaystyle \frac{\pi^{2}\hbar^{2}}{2m\left(\frac{a}{2}+\epsilon\right)^{2}}\ \ \ \ \ (3)$ $\displaystyle E_{r}$ $\displaystyle =$ $\displaystyle \frac{\pi^{2}\hbar^{2}}{2m\left(\frac{a}{2}-\epsilon\right)^{2}} \ \ \ \ \ (4)$

Thus ${E_{l} so the ground state confines the particle to the left well. The wave function for the ground state is

$\displaystyle \psi\left(x\right)=\begin{cases} \sqrt{\frac{2}{\frac{a}{2}+\epsilon}}\sin\frac{\pi}{\frac{a}{2}+\epsilon}x & 0\le x<\frac{a}{2}+\epsilon\\ 0 & \frac{a}{2}+\epsilon

The plot looks like this:

Now for the general case where ${f\left(t\right)}$ is finite. In this case we can write the wave functions as

$\displaystyle \psi\left(x\right)=\begin{cases} Ae^{ikx}+Be^{-ikx} & 0\le x<\frac{a}{2}+\epsilon\\ Ce^{ikx}+De^{-ikx} & \frac{a}{2}+\epsilon

where

$\displaystyle k\equiv\frac{\sqrt{2mE}}{\hbar} \ \ \ \ \ (7)$

The barriers at ${x=0}$ and ${x=a}$ are still infinite so the wave function must be zero there, giving

 $\displaystyle A$ $\displaystyle =$ $\displaystyle -B\ \ \ \ \ (8)$ $\displaystyle C$ $\displaystyle =$ $\displaystyle -De^{-2ika} \ \ \ \ \ (9)$

The wave function must be continuous at the barrier, so we get

$\displaystyle A\left(e^{ik\left(\frac{a}{2}+\epsilon\right)}-e^{-ik\left(\frac{a}{2}+\epsilon\right)}\right)=D\left(-e^{-2ika}e^{ik\left(\frac{a}{2}+\epsilon\right)}+e^{-ik\left(\frac{a}{2}+\epsilon\right)}\right) \ \ \ \ \ (10)$

Finally, we can analyze the derivative at the barrier in the same way we did for the delta function well and we get

 $\displaystyle -\frac{\hbar^{2}}{2m}\int_{\frac{a}{2}+\epsilon-\beta}^{\frac{a}{2}+\epsilon+\beta}\frac{d^{2}\psi}{dx^{2}}dx+f\left(t\right)\int_{\frac{a}{2}+\epsilon-\beta}^{\frac{a}{2}+\epsilon+\beta}\delta(x)\psi dx$ $\displaystyle =$ $\displaystyle E\int_{\frac{a}{2}+\epsilon-\beta}^{\frac{a}{2}+\epsilon+\beta}\psi dx\ \ \ \ \ (11)$ $\displaystyle -\frac{\hbar^{2}}{2m}\left.\frac{d\psi}{dx}\right|_{\frac{a}{2}+\epsilon-\beta}^{\frac{a}{2}+\epsilon+\beta}+f\left(t\right)\psi\left(\frac{a}{2}+\epsilon\right)$ $\displaystyle =$ $\displaystyle E\int_{\frac{a}{2}+\epsilon-\beta}^{\frac{a}{2}+\epsilon+\beta}\psi dx \ \ \ \ \ (12)$

The integral on the RHS goes to zero as ${\beta\rightarrow0}$ since ${\psi}$ is finite, so

 $\displaystyle A\frac{2mf\left(t\right)}{\hbar^{2}}\left(e^{ik\left(\frac{a}{2}+\epsilon\right)}-e^{-ik\left(\frac{a}{2}+\epsilon\right)}\right)$ $\displaystyle =$ $\displaystyle -ikD\left(e^{-2ika}e^{ik\left(\frac{a}{2}+\epsilon\right)}+e^{-ik\left(\frac{a}{2}+\epsilon\right)}\right)-\nonumber$ $\displaystyle$ $\displaystyle$ $\displaystyle ikA\left(e^{ik\left(\frac{a}{2}+\epsilon\right)}+e^{-ik\left(\frac{a}{2}+\epsilon\right)}\right) \ \ \ \ \ (13)$

If we now define

 $\displaystyle z$ $\displaystyle \equiv$ $\displaystyle ka\ \ \ \ \ (14)$ $\displaystyle \delta$ $\displaystyle \equiv$ $\displaystyle \frac{2\epsilon}{a}\ \ \ \ \ (15)$ $\displaystyle k\left(\frac{a}{2}+\epsilon\right)$ $\displaystyle =$ $\displaystyle \frac{1}{2}z\left(1+\delta\right) \ \ \ \ \ (16)$

we get, transforming the complex exponentials to trig functions

$\displaystyle A\frac{4imf\left(t\right)}{\hbar^{2}}\sin\left[\frac{1}{2}z\left(1+\delta\right)\right]=-2ikDe^{-iz}\cos\left[\frac{1}{2}z\left(1-\delta\right)\right]-2ikA\cos\left[\frac{1}{2}z\left(1+\delta\right)\right] \ \ \ \ \ (17)$

Multiplying through by ${a}$ and defining

$\displaystyle T\equiv\frac{maf\left(t\right)}{\hbar^{2}} \ \ \ \ \ (18)$

we get

$\displaystyle 2AT\sin\left[\frac{1}{2}z\left(1+\delta\right)\right]=-zDe^{-iz}\cos\left[\frac{1}{2}z\left(1-\delta\right)\right]-zA\cos\left[\frac{1}{2}z\left(1+\delta\right)\right] \ \ \ \ \ (19)$

We can write 10 as

 $\displaystyle 2iA\sin\left[\frac{1}{2}z\left(1+\delta\right)\right]$ $\displaystyle =$ $\displaystyle 2ie^{-iz}D\sin\left[\frac{1}{2}z\left(1-\delta\right)\right]\ \ \ \ \ (20)$ $\displaystyle A$ $\displaystyle =$ $\displaystyle e^{-iz}D\frac{\sin\left[\frac{1}{2}z\left(1-\delta\right)\right]}{\sin\left[\frac{1}{2}z\left(1+\delta\right)\right]} \ \ \ \ \ (21)$

Substituting this into 19, multiplying through by ${\sin\left[\frac{1}{2}z\left(1+\delta\right)\right]}$ and cancelling terms we get

 $\displaystyle 2T\sin\left[\frac{1}{2}z\left(1-\delta\right)\right]\sin\left[\frac{1}{2}z\left(1+\delta\right)\right]$ $\displaystyle =$ $\displaystyle -z\left[\cos\left[\frac{1}{2}z\left(1-\delta\right)\right]\sin\left[\frac{1}{2}z\left(1+\delta\right)\right]+\sin\left[\frac{1}{2}z\left(1-\delta\right)\right]\cos\left[\frac{1}{2}z\left(1+\delta\right)\right]\right]\ \ \ \ \ (22)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -z\sin\left[\frac{1}{2}z\left(1-\delta\right)+\frac{1}{2}z\left(1+\delta\right)\right]\ \ \ \ \ (23)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -z\sin z \ \ \ \ \ (24)$

The LHS can be transformed using

 $\displaystyle 2\sin\left[\frac{1}{2}z\left(1-\delta\right)\right]\sin\left[\frac{1}{2}z\left(1+\delta\right)\right]$ $\displaystyle =$ $\displaystyle \cos\left[\frac{1}{2}z\left(1-\delta\right)-\frac{1}{2}z\left(1+\delta\right)\right]-\cos\left[\frac{1}{2}z\left(1-\delta\right)+\frac{1}{2}z\left(1+\delta\right)\right]\ \ \ \ \ (25)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \cos z\delta-\cos z \ \ \ \ \ (26)$

Putting it together we get the transcendental equation for the ground state energy

$\displaystyle z\sin z=T\left(\cos z-\cos z\delta\right) \ \ \ \ \ (27)$

The energy is

$\displaystyle E=\frac{\hbar^{2}k^{2}}{2m}=\frac{\hbar^{2}z^{2}}{2ma^{2}} \ \ \ \ \ (28)$

We can solve 27 graphically or numerically. For ${\delta=0}$, here are plots of ${z\sin z}$ (red) and ${T\left(\cos z-1\right)}$ (green) for 3 values of ${T}$:

Discounting the trivial solutions at ${z=0}$ and ${z=2\pi}$ which are valid for all ${T}$, we see that the intersection point moves from ${z=\pi}$ out to ${z=2\pi}$ as ${T}$ increases. This is equivalent to the energy moving from ${\frac{\hbar^{2}\pi^{2}}{2ma^{2}}}$ up to ${\frac{4\hbar^{2}\pi^{2}}{2ma^{2}}}$. The first energy is that of a well of width ${a}$ while the second energy is that of a well of width ${\frac{a}{2}}$, which is what we’d expect. As the barrier gets stronger, the ground state energy approaches that of a well of half the width of the original.

Using Maple’s fsolve command we can work out some other values of ${z}$ for ${\delta=0.01}$:

 ${T}$ ${z}$ 1 3.673 5 4.760 20 5.720 100 6.135 1000 6.215

# Field operators for the infinite square well

Reference: Tom Lancaster and Stephen J. Blundell, Quantum Field Theory for the Gifted Amateur, (Oxford University Press, 2014), Section 4.1.

The creation and annihilation operators ${\hat{a}_{\mathbf{p}}^{\dagger}}$ and ${\hat{a}_{\mathbf{p}}}$ create and annihilate a particle in a specific momentum state so, because of the uncertainty principle, such states are completely unlocalized in position. We can construct analogous operators that create and annihilate a particle at a specific position. Such operators are called field operators. As you would expect, a field operator, operating at a precise position, must be completely unlocalized in momentum.

The creation and annihilation field operators for a particle in a 3-d infinite square well are defined as

 $\displaystyle \hat{\psi}^{\dagger}\left(\mathbf{x}\right)$ $\displaystyle \equiv$ $\displaystyle \frac{1}{\sqrt{\mathcal{V}}}\sum_{\mathbf{p}}\hat{a}_{\mathbf{p}}^{\dagger}e^{-i\mathbf{p}\cdot\mathbf{x}}\ \ \ \ \ (1)$ $\displaystyle \hat{\psi}\left(\mathbf{x}\right)$ $\displaystyle \equiv$ $\displaystyle \frac{1}{\sqrt{\mathcal{V}}}\sum_{\mathbf{p}}\hat{a}_{\mathbf{p}}e^{i\mathbf{p}\cdot\mathbf{x}} \ \ \ \ \ (2)$

where ${\mathcal{V}}$ is the volume of the square well.

To see that they actually do create and annihilate a particle at position ${\mathbf{x}}$ we’ll have a look at their effect when operating on particular states.

First, we’ll apply ${\hat{\psi}^{\dagger}\left(x\right)}$ to the vacuum state.

 $\displaystyle \hat{\psi}^{\dagger}\left(\mathbf{x}\right)\left|0\right\rangle$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{\mathcal{V}}}\sum_{\mathbf{p}}e^{-i\mathbf{p}\cdot\mathbf{x}}\hat{a}_{\mathbf{p}}^{\dagger}\left|0\right\rangle \ \ \ \ \ (3)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{\mathcal{V}}}\sum_{\mathbf{p}}e^{-i\mathbf{p}\cdot\mathbf{x}}\left|\mathbf{p}\right\rangle \ \ \ \ \ (4)$

We can now insert the unit operator in the form

$\displaystyle 1=\int d^{3}\mathbf{y}\left|\mathbf{y}\right\rangle \left\langle \mathbf{y}\right| \ \ \ \ \ (5)$

and we get

style=”text-align:center”> $\displaystyle =$
 $\displaystyle \hat{\psi}^{\dagger}\left(\mathbf{x}\right)\left|0\right\rangle$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{\mathcal{V}}}\int d^{3}\mathbf{y}\sum_{\mathbf{p}}e^{-i\mathbf{p}\cdot\mathbf{x}}\left|\mathbf{y}\right\rangle \left\langle \mathbf{y}\left|\mathbf{p}\right.\right\rangle \ \ \ \ \ (6)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{\mathcal{V}}}\int d^{3}\mathbf{y}\sum_{\mathbf{p}}e^{-i\mathbf{p}\cdot\mathbf{x}}\left|\mathbf{y}\right\rangle \left[\frac{1}{\sqrt{\mathcal{V}}}e^{i\mathbf{p}\cdot\mathbf{y}}\right]\ \ \ \ \ (7)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \int d^{3}\mathbf{y}\left[\frac{1}{\mathcal{V}}\sum_{\mathbf{p}}e^{i\mathbf{p}\cdot\left(\mathbf{y}-\mathbf{x}\right)}\right]\left|\mathbf{y}\right\rangle \ \ \ \ \ (8)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \int d^{3}\mathbf{y}\delta^{\left(3\right)}\left(\mathbf{y}-\mathbf{x}\right)\left|\mathbf{y}\right\rangle \ \ \ \ \ (9)$ $\displaystyle$ $\displaystyle \left|\mathbf{x}\right\rangle \ \ \ \ \ (10)$

So the creation operator ${\hat{\psi}^{\dagger}\left(\mathbf{x}\right)}$ operating on the vacuum state creates a single particle at position ${\mathbf{x}}$.

We can now try the annihilation operator ${\hat{\psi}\left(\mathbf{y}\right)}$ operating on a one-particle state ${\left|\mathbf{x}\right\rangle }$. We get

 $\displaystyle \hat{\psi}\left(\mathbf{y}\right)\left|\mathbf{x}\right\rangle$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{\mathcal{V}}}\sum_{\mathbf{p}}\hat{a}_{\mathbf{p}}e^{i\mathbf{p}\cdot\mathbf{y}}\left|\mathbf{x}\right\rangle \ \ \ \ \ (11)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{\mathcal{V}}}\sum_{\mathbf{p}.\mathbf{q}}\hat{a}_{\mathbf{p}}e^{i\mathbf{p}\cdot\mathbf{y}}\left|\mathbf{q}\right\rangle \left\langle \mathbf{q}\left|\mathbf{x}\right.\right\rangle \ \ \ \ \ (12)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{\mathcal{V}}}\sum_{\mathbf{p}.\mathbf{q}}\hat{a}_{\mathbf{p}}e^{i\mathbf{p}\cdot\mathbf{y}}\left|\mathbf{q}\right\rangle \left[\frac{1}{\sqrt{\mathcal{V}}}e^{-i\mathbf{q}\cdot\mathbf{x}}\right]\ \ \ \ \ (13)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{1}{\mathcal{V}}\sum_{\mathbf{p},\mathbf{q}}e^{i\left(\mathbf{p}\cdot\mathbf{y}-\mathbf{q}\cdot\mathbf{x}\right)}\hat{a}_{\mathbf{p}}\left|\mathbf{q}\right\rangle \ \ \ \ \ (14)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{1}{\mathcal{V}}\sum_{\mathbf{p},\mathbf{q}}e^{i\left(\mathbf{p}\cdot\mathbf{y}-\mathbf{q}\cdot\mathbf{x}\right)}\delta_{\mathbf{pq}}\left|0\right\rangle \ \ \ \ \ (15)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{1}{\mathcal{V}}\sum_{\mathbf{p}}e^{i\mathbf{p}\cdot\left(\mathbf{y}-\mathbf{x}\right)}\left|0\right\rangle \ \ \ \ \ (16)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \delta^{\left(3\right)}\left(\mathbf{y}-\mathbf{x}\right)\left|0\right\rangle \ \ \ \ \ (17)$

where in the fourth line the momentum annihilation operator ${\hat{a}_{\mathbf{p}}}$ operating on ${\left|\mathbf{q}\right\rangle }$ produces the vacuum state only if ${\mathbf{p}=\mathbf{q}}$; otherwise it produces zero.

Thus ${\hat{\psi}\left(\mathbf{y}\right)\left|\mathbf{x}\right\rangle }$ produces the vacuum state if ${\mathbf{y}=\mathbf{x}}$ and zero otherwise.

# The adiabatic approximation in quantum mechanics

Reference: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 10.1.

The adiabatic approximation in quantum mechanics is a method by which approximate solutions to the time dependent Schrödinger equation can be found. The method works in cases where the hamiltonian changes slowly by comparison with the natural, internal frequency of the wave function. For example, the general solution to the time-independent Schrödinger equation can be written as

$\displaystyle \Psi\left(x,t\right)=\sum_{n}c_{n}\psi_{n}\left(x\right)e^{-iE_{n}t/\hbar} \ \ \ \ \ (1)$

where the ${\psi_{n}\left(x\right)}$ are the eigenfunctions of the hamiltonian with eigenvalues (energies) ${E_{n}}$. The internal frequency of eigenfunction ${n}$ is ${E_{n}/\hbar}$, so if the variations in the hamiltonian have frequency components much lower than this, we can use the adiabatic approximation.

The adiabatic theorem states that in a system with non-degenerate energy levels, if we start with the system in level ${n}$ of the original hamiltonian and then undergo an adiabatic process that takes us to some final hamiltonian, then the system will be in level ${n}$ of the final hamiltonian, although the wave function may pick up a phase factor along the way.

A common example of an adiabatic process is an infinite square well that starts with width ${a}$ that increases slowly over time with constant speed ${v}$ so that the width is given by

$\displaystyle w\left(t\right)=a+vt \ \ \ \ \ (2)$

for ${t\ge0}$. If we let the wall expand to twice its original width, then we can take as the ‘external’ time ${T_{e}}$ the time it takes to complete its expansion:

$\displaystyle T_{e}=\frac{a}{v} \ \ \ \ \ (3)$

The ‘internal’ time ${T_{i}}$ can be the period of the phase factor ${e^{-iE_{n}t/\hbar}}$ in the starting state, so we would have

 $\displaystyle \frac{E_{n}T_{i}}{\hbar}$ $\displaystyle =$ $\displaystyle 2\pi\ \ \ \ \ (4)$ $\displaystyle T_{i}$ $\displaystyle =$ $\displaystyle \frac{2\pi\hbar}{E_{n}} \ \ \ \ \ (5)$

However, it turns out that the time dependent Schrödinger equation can be solved exactly in this case, with the ${n}$th eigenfunction given by

$\displaystyle \Phi_{n}\left(x,t\right)=\sqrt{\frac{2}{w}}\sin\left(\frac{n\pi}{w}x\right)e^{i\left(mvx^{2}-2E_{n}^{i}at\right)/2\hbar w} \ \ \ \ \ (6)$

where ${E_{n}^{i}}$ is the energy of level ${n}$ in the starting well, with width ${a}$:

$\displaystyle E_{n}^{i}=\frac{n^{2}\pi^{2}\hbar^{2}}{2ma^{2}} \ \ \ \ \ (7)$

We can verify by direct differentiation that 6 satisfies the time dependent Schrödinger equation

$\displaystyle -\frac{\hbar^{2}}{2m}\frac{\partial^{2}\Phi_{n}}{\partial x^{2}}=i\hbar\frac{\partial\Phi_{n}}{\partial t} \ \ \ \ \ (8)$

The calculation gets very messy (remember ${w}$ depends on ${t}$) so it’s best to use Maple to do it, and we get

$\displaystyle i\hbar\frac{\partial\Phi_{n}}{\partial t}=e^{i\left(mvx^{2}-2E_{n}^{i}at\right)/2\hbar w}\frac{\sqrt{2}}{w^{5/2}}\left[\frac{1}{2}\sin\left(\frac{n\pi}{w}x\right)\left(\frac{\left(n\pi\hbar\right)^{2}}{m}+m\left(vx\right)^{2}-i\hbar vw\right)-i\cos\left(\frac{n\pi}{w}x\right)\pi n\hbar vx\right] \ \ \ \ \ (9)$

Fortunately, we get the same expression for ${-\frac{\hbar^{2}}{2m}\frac{\partial^{2}\Phi_{n}}{\partial x^{2}}}$ so the Schrödinger equation is satisfied.

6 is the wave function for a single energy level, so we can get a general solution that is a superposition of energy levels in the usual way:

$\displaystyle \Psi\left(x,t\right)=\sum_{n=1}^{\infty}c_{n}\Phi_{n}\left(x,t\right) \ \ \ \ \ (10)$

In this case, all the time dependence is included in the ${\Phi_{n}}$, so the coefficients ${c_{n}}$ are true constants, independent of both space and time. The ${\Phi_{n}}$ are orthonormal at each instant of time, since

 $\displaystyle \int_{0}^{w}\Phi_{j}^*\Phi_{n}dx$ $\displaystyle =$ $\displaystyle \frac{2}{w}e^{2iat\left(E_{j}^{i}-E_{n}^{i}\right)/2\hbar w}\int_{0}^{w}\sin\left(\frac{j\pi}{w}x\right)\sin\left(\frac{n\pi}{w}x\right)dx\ \ \ \ \ (11)$ $\displaystyle$ $\displaystyle =$ $\displaystyle e^{2iat\left(E_{j}^{i}-E_{n}^{i}\right)/2\hbar w}\delta_{jn}\ \ \ \ \ (12)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \delta_{jn} \ \ \ \ \ (13)$

We can therefore use orthonormality to get an expression for ${c_{n}}$ by multiplying both sides by ${\Phi_{j}^*}$ and integrating. Since the ${c_{n}}$ are independent of time, we can do this at ${t=0}$ when ${w=a}$.

$\displaystyle c_{j}=\int_{0}^{a}\Phi_{j}^*\left(x,0\right)\Psi\left(x,0\right)dx \ \ \ \ \ (14)$

If the particle starts out in the ground state, then

 $\displaystyle \Psi\left(x,0\right)$ $\displaystyle =$ $\displaystyle \sqrt{\frac{2}{a}}\sin\left(\frac{\pi}{a}x\right)\ \ \ \ \ (15)$ $\displaystyle c_{n}$ $\displaystyle =$ $\displaystyle \frac{2}{a}\int_{0}^{a}e^{-imvx^{2}/2\hbar a}\sin\left(\frac{n\pi}{a}x\right)\sin\left(\frac{\pi}{a}x\right)dx \ \ \ \ \ (16)$

Substituting ${z=\pi x/a}$ we get

 $\displaystyle c_{n}$ $\displaystyle =$ $\displaystyle \frac{2}{\pi}\int_{0}^{\pi}e^{-i\alpha z^{2}}\sin\left(nz\right)\sin\left(z\right)dz\ \ \ \ \ (17)$ $\displaystyle \alpha$ $\displaystyle \equiv$ $\displaystyle \frac{mva}{2\pi^{2}\hbar} \ \ \ \ \ (18)$

So far, all this is exact. To use the adiabatic approximation, we need estimates of ${T_{e}}$ and ${T_{i}}$. We can get ${T_{e}}$ from 3. For ${T_{i}}$ we can look at 6 at ${x=0}$ and find the value of ${t}$ that makes the argument of the exponential advance by ${2\pi}$. Actually, because of the signs, the phase actually goes backwards as ${t}$ gets larger, but the principle is the same; we just need to find the value of ${t}$ at which the argument changes by ${2\pi}$.

 $\displaystyle \frac{2E_{1}^{i}aT_{i}}{2\hbar\left(a+vT_{i}\right)}$ $\displaystyle =$ $\displaystyle 2\pi\ \ \ \ \ (19)$ $\displaystyle T_{i}$ $\displaystyle =$ $\displaystyle \frac{2\pi\hbar a}{E_{1}^{i}a-2\pi\hbar v} \ \ \ \ \ (20)$

Dividing top and bottom by ${v}$ and using 3 we get

$\displaystyle T_{i}=\frac{2\pi\hbar T_{e}}{E_{1}^{i}T_{e}-2\pi\hbar} \ \ \ \ \ (21)$

To satisfy the adiabatic condition, we need ${T_{e}\gg T_{i}}$ so, using 7

 $\displaystyle \frac{2\pi\hbar}{E_{1}^{i}T_{e}-2\pi\hbar}$ $\displaystyle \ll$ $\displaystyle 1\ \ \ \ \ (22)$ $\displaystyle E_{1}^{i}T_{e}$ $\displaystyle \gg$ $\displaystyle 4\pi\hbar\ \ \ \ \ (23)$ $\displaystyle \frac{8mav}{\pi\hbar}$ $\displaystyle \ll$ $\displaystyle 1 \ \ \ \ \ (24)$

Apart from the numerical factors which don’t differ too drastically, we see from 18 that this effectively requires ${\alpha\ll1}$. Using this approximation, we can evaluate the integral in 17 to get

 $\displaystyle c_{n}$ $\displaystyle \approx$ $\displaystyle \frac{2}{\pi}\int_{0}^{\pi}\sin\left(nz\right)\sin\left(z\right)dz\ \ \ \ \ (25)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \delta_{1n} \ \ \ \ \ (26)$

Thus the system remains in the ground state ${n=1}$ as the wall expands, which is what the adiabatic theorem predicts. The wave function is therefore

$\displaystyle \Phi_{1}\left(x,t\right)\approx\sqrt{\frac{2}{w}}\sin\left(\frac{\pi}{w}x\right)e^{i\left(mvx^{2}-2E_{1}^{i}at\right)/2\hbar w} \ \ \ \ \ (27)$

The phase factor in this wave function is the exponential factor that doesn’t depend on ${x}$, that is

$\displaystyle \theta\left(t\right)=-\frac{2E_{1}^{i}at}{2\hbar w}=-\frac{\pi^{2}\hbar t}{2ma\left(a+vt\right)} \ \ \ \ \ (28)$

The instantaneous eigenvalue in the ground state is the original eigenvalue with the well width ${a}$ replaced by the dynamic width ${w}$:

$\displaystyle E_{1}\left(t\right)=\frac{\pi^{2}\hbar^{2}}{2mw^{2}} \ \ \ \ \ (29)$

If we integrate this over the time the wall has moved so far we get

 $\displaystyle \int_{0}^{t}E_{1}\left(t'\right)dt'$ $\displaystyle =$ $\displaystyle \frac{\pi^{2}\hbar^{2}}{2m}\int_{0}^{t}\frac{dt'}{\left(a+vt'\right)^{2}}\ \ \ \ \ (30)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{\pi^{2}\hbar^{2}t}{2ma\left(a+vt\right)}\ \ \ \ \ (31)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -\hbar\theta\left(t\right) \ \ \ \ \ (32)$

# Time-dependent perturbation of the infinite square well

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Chapter 9, Post 18.

Here’s another example of a perturbation in a multi-state time dependent system. We start with a particle of mass ${m}$ in the ground state ${\psi_{1}}$ of the one dimensional infinite square well. At time ${t=0}$ the potential is changed to

$\displaystyle V\left(x\right)=\begin{cases} V_{0} & 0\le x\le\frac{a}{2}\\ 0 & \frac{a}{2}

and at time ${t=T}$ the potential reverts to the unperturbed square well. We’d like to find the probability that the particle makes a transition to state ${\psi_{2}}$.

To do this, we can use first-order perturbation theory to find the coefficient ${c_{2}\left(T\right)}$ in the wave function

$\displaystyle \Psi\left(t\right)=\sum_{n}c_{n}\left(t\right)\psi_{n}e^{-iE_{n}t/\hbar} \ \ \ \ \ (2)$

The first order coefficient is given by

$\displaystyle c_{2}\left(t\right)=-\frac{i}{\hbar}\int_{0}^{T}H'_{12}e^{i\left(E_{2}-E_{1}\right)t/\hbar}dt \ \ \ \ \ (3)$

To get the matrix element ${H'_{12}}$ we need the unperturbed wave functions:

$\displaystyle \psi_{n}\left(x\right)=\sqrt{\frac{2}{a}}\sin\frac{n\pi x}{a} \ \ \ \ \ (4)$

We get

 $\displaystyle H'_{12}$ $\displaystyle =$ $\displaystyle \frac{2V_{0}}{a}\int_{0}^{a/2}\sin\frac{\pi x}{a}\sin\frac{2\pi x}{a}dx\ \ \ \ \ (5)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{4V_{0}}{3\pi}\ \ \ \ \ (6)$ $\displaystyle c_{2}$ $\displaystyle =$ $\displaystyle \frac{4V_{0}}{3\pi}\frac{e^{i\left(E_{2}-E_{1}\right)T/\hbar}-1}{E_{2}-E_{1}}\ \ \ \ \ (7)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{8iV_{0}}{3\pi\left(E_{2}-E_{1}\right)}e^{i\left(E_{2}-E_{1}\right)T/2\hbar}\sin\frac{\left(E_{2}-E_{2}\right)T}{2\hbar} \ \ \ \ \ (8)$

Using the energy levels of the infinite square well

$\displaystyle E_{n}=\frac{\left(n\pi\hbar\right)^{2}}{2ma^{2}} \ \ \ \ \ (9)$

we get

 $\displaystyle E_{2}-E_{1}$ $\displaystyle =$ $\displaystyle \frac{3\pi^{2}\hbar^{2}}{2ma^{2}}\ \ \ \ \ (10)$ $\displaystyle \left|c_{2}\left(T\right)\right|^{2}$ $\displaystyle =$ $\displaystyle \left(\frac{16ma^{2}V_{0}}{9\pi^{3}\hbar^{2}}\right)^{2}\sin^{2}\left(\frac{3\pi^{2}\hbar T}{4ma^{2}}\right) \ \ \ \ \ (11)$

# WKB approximation

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 8.1.

We now look at another approximation technique used in quantum mechanics. This is the WKB (named for the German physicist Gregor Wentzel (1898 – 1978), the Dutch physicist Hendrik Kramers (1894 – 1952) and the French physicist Léon Brillouin (1889 – 1969)) approximation, which is a mathematical technique applicable to one-dimensional differential equations.

$\displaystyle -\frac{\hbar^{2}}{2m}\frac{d^{2}\psi}{dx^{2}}+V\left(x\right)\psi=E\psi \ \ \ \ \ (1)$

and rewrite it as

 $\displaystyle \frac{d^{2}\psi}{dx^{2}}$ $\displaystyle =$ $\displaystyle -\frac{2m\left[E-V\left(x\right)\right]}{\hbar^{2}}\psi\ \ \ \ \ (2)$ $\displaystyle$ $\displaystyle \equiv$ $\displaystyle -\frac{p^{2}}{\hbar^{2}}\psi \ \ \ \ \ (3)$

where ${p=\sqrt{2m\left[E-V\left(x\right)\right]}}$ is the classical formula for the momentum of a particle with total energy ${E}$ moving in a one-dimensional potential ${V\left(x\right)}$, provided we assume ${E\ge V\left(x\right)}$ for all ${x}$.

In general, the wave function ${\psi\left(x\right)}$ is a complex function so we can write it in complex exponential form as

$\displaystyle \psi\left(x\right)=A\left(x\right)e^{i\phi\left(x\right)} \ \ \ \ \ (4)$

where ${A\left(x\right)}$ is the amplitude and ${\phi\left(x\right)}$ is the phase, both of which are real functions. In this form, we have (using a prime to denote a derivative with respect to ${x}$):

 $\displaystyle \psi'$ $\displaystyle =$ $\displaystyle \left(A'+iA\phi'\right)e^{i\phi}\ \ \ \ \ (5)$ $\displaystyle \psi^{\prime\prime}$ $\displaystyle =$ $\displaystyle \left(A^{\prime\prime}+iA'\phi'+iA\phi^{\prime\prime}+iA'\phi'-A\left(\phi'\right)^{2}\right)e^{i\phi}\ \ \ \ \ (6)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left(A^{\prime\prime}+2iA'\phi'+iA\phi^{\prime\prime}-A\left(\phi'\right)^{2}\right)e^{i\phi} \ \ \ \ \ (7)$

Inserting this into 3 and cancelling off ${e^{i\phi}}$ we get

$\displaystyle A^{\prime\prime}+2iA'\phi'+iA\phi^{\prime\prime}-A\left(\phi'\right)^{2}=-\frac{p^{2}}{\hbar^{2}}A \ \ \ \ \ (8)$

Everything in this equation is real apart from ${i}$ itself, so we can separate this equation into its real and imaginary parts to get two differential equations:

 $\displaystyle A^{\prime\prime}-A\left(\phi'\right)^{2}$ $\displaystyle =$ $\displaystyle -\frac{p^{2}}{\hbar^{2}}A\ \ \ \ \ (9)$ $\displaystyle 2A'\phi'+A\phi^{\prime\prime}$ $\displaystyle =$ $\displaystyle 0 \ \ \ \ \ (10)$

The second of these equations can always be solved as

 $\displaystyle 2A'\phi'+A\phi^{\prime\prime}$ $\displaystyle =$ $\displaystyle \frac{1}{A}\left[A^{2}\phi'\right]'=0\ \ \ \ \ (11)$ $\displaystyle \left[A^{2}\phi'\right]'$ $\displaystyle =$ $\displaystyle 0\ \ \ \ \ (12)$ $\displaystyle A^{2}\phi'$ $\displaystyle =$ $\displaystyle C^{2} \ \ \ \ \ (13)$

where ${C}$ is a constant, which may be complex (since ${\phi'}$ could be negative). We can therefore write the amplitude as

$\displaystyle A\left(x\right)=\frac{C}{\sqrt{\phi'}} \ \ \ \ \ (14)$

We’ve taken the positive square root, since any difference in sign can be accounted for by changing the integration constant ${C}$.

Equation 9 can’t be solved in general since it depends on ${V\left(x\right)}$ which could in principle be anything. However, we can rewrite it as

$\displaystyle \frac{A^{\prime\prime}}{A}=\left(\phi'\right)^{2}-\frac{p^{2}}{\hbar^{2}} \ \ \ \ \ (15)$

This is where the approximation comes in. If ${V\left(x\right)}$ were constant then we can solve the Schrödinger equation exactly (as we did with the finite square well and finite square barrier). If ${E>V}$ we get a travelling wave with a constant amplitude ${A}$ and constant wavelength ${\lambda=2\pi\hbar/p}$, while if ${E we get an exponential decay with a characteristic decay length of ${\ell=\hbar/\sqrt{2m\left[V\left(x\right)-E\right]}}$. Now suppose that ${V\left(x\right)}$ is not constant, but that it varies slowly compared to either ${\lambda}$ or ${\ell}$. In this case, we’d expect that the wave function is close to ${\psi}$ for a constant potential, except that its amplitude and phase will vary slightly. The approximation comes by assuming that the variation in amplitude is small enough that the derivatives of ${A\left(x\right)}$ are negligible compared to ${A\left(x\right)}$ itself. That is, we can set the LHS of 15 to zero, which allows us to write down a solution of the RHS:

 $\displaystyle \phi'$ $\displaystyle =$ $\displaystyle \pm\frac{p}{\hbar}\ \ \ \ \ (16)$ $\displaystyle \phi\left(x\right)$ $\displaystyle =$ $\displaystyle \pm\frac{1}{\hbar}\int p\; dx \ \ \ \ \ (17)$

This is written as an indefinite integral, since the constant of integration (${K}$ say) can be absorbed into ${C}$:

 $\displaystyle \phi\left(x\right)$ $\displaystyle =$ $\displaystyle \pm\frac{1}{\hbar}\int p\; dx+K\ \ \ \ \ (18)$ $\displaystyle \psi\left(x\right)$ $\displaystyle =$ $\displaystyle A\left(x\right)e^{i\phi\left(x\right)}\ \ \ \ \ (19)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{C}{\sqrt{\phi'}}e^{iK}e^{\pm i\int p\; dx/\hbar}\ \ \ \ \ (20)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{C\sqrt{\hbar}e^{iK}}{\sqrt{p\left(x\right)}}e^{\pm i\int p\; dx/\hbar}\ \ \ \ \ (21)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{C_{1}}{\sqrt{p\left(x\right)}}e^{\pm i\int p\; dx/\hbar} \ \ \ \ \ (22)$

where

$\displaystyle C_{1}\equiv C\sqrt{\hbar}e^{iK} \ \ \ \ \ (23)$

The WKB approximation for the wave function is thus

$\displaystyle \boxed{\psi\left(x\right)\approx\frac{C_{1}}{\sqrt{p\left(x\right)}}e^{\pm i\int p\; dx/\hbar}} \ \ \ \ \ (24)$

Provided we can do the integral in the exponent, we can find the wave function and, by imposing boundary conditions, we can often find the allowed energies as well.

Example We have an infinite square well with a shelf in the well, with the potential given by

$\displaystyle V\left(x\right)=\begin{cases} V_{0} & 0a \end{cases} \ \ \ \ \ (25)$

The classical momentum is

$\displaystyle p\left(x\right)=\begin{cases} \sqrt{2m\left(E-V_{0}\right)} & 0

so we have

$\displaystyle \int_{0}^{x}p\; dx'=\begin{cases} x\sqrt{2m\left(E-V_{0}\right)} & 0

We’ve put limits on the integral with the understanding that any constant of integration is absorbed into the overall constant that multiplies the amplitude. This constant will be determined by normalization as usual.

We thus have

$\displaystyle \phi\left(x\right)=\begin{cases} \pm\frac{x}{\hbar}\sqrt{2m\left(E-V_{0}\right)} & 0

We can now write the wave function as

$\displaystyle \psi\left(x\right)=\begin{cases} \frac{1}{\left[2m\left(E-V_{0}\right)\right]^{1/4}}\left(C_{+}e^{i\phi\left(x\right)}+C_{-}e^{-i\phi\left(x\right)}\right) & 0

where the constants ${C_{+}}$ and ${C_{-}}$ can be determined from boundary conditions and normalization as usual. (Actually, it’s worth pointing out here that although ${\phi\left(x\right)}$ is continuous at ${x=\frac{a}{2}}$, the WKB wave function ${\psi\left(x\right)}$ is not continuous at this point. This violates one of the central conditions imposed on wave functions. The discontinuity arises from the discontinuity in the potential at this point, which violates one of the assumptions of WKB theory: that the potential varies slowly. Thus at the boundaries, the WKB assumptions don’t hold, so it’s a bit surprising that WKB gives a reasonable result for this problem since there are infinite discontinuities at ${x=0}$ and ${x=a}$.)

To get the allowed energies, we can impose continuity at the end points of the well. That is

$\displaystyle \psi\left(0\right)=\psi\left(a\right)=0 \ \ \ \ \ (30)$

To use this condition, it’s easier to write the wave function using trigonometric functions:

$\displaystyle \psi\left(x\right)=\begin{cases} \frac{1}{\left[2m\left(E-V_{0}\right)\right]^{1/4}}\left(C_{1}\sin\phi\left(x\right)+C_{2}\cos\phi\left(x\right)\right) & 0

The condition ${\psi\left(0\right)=0}$ gives us ${C_{2}=0}$ while the other condition ${\psi\left(a\right)=0}$ tells us

 $\displaystyle \sin\left(\phi\left(a\right)\right)$ $\displaystyle =$ $\displaystyle 0\ \ \ \ \ (32)$ $\displaystyle \phi\left(a\right)$ $\displaystyle =$ $\displaystyle n\pi\ \ \ \ \ (33)$ $\displaystyle \left[\sqrt{2m\left(E_{n}-V_{0}\right)}+\sqrt{2mE_{n}}\right]\frac{a}{2\hbar}$ $\displaystyle =$ $\displaystyle n\pi \ \ \ \ \ (34)$

Doing a bit of algebra we get

$\displaystyle E-\frac{V_{0}}{2}+\sqrt{E\left(E-V_{0}\right)}=\frac{n^{2}\pi^{2}\hbar^{2}}{ma^{2}} \ \ \ \ \ (35)$

The quantity on the RHS is twice the allowed energies ${E_{n}^{0}}$ of the ordinary infinite square well. Solving this for ${E}$ we get

$\displaystyle E=E_{n}^{0}+\frac{V_{0}}{2}+\frac{V_{0}^{2}}{16E_{n}^{0}} \ \ \ \ \ (36)$

Griffiths shows in his Example 6.1 that first-order perturbation theory gives a result of

$\displaystyle E=E_{n}^{0}+\frac{V_{0}}{2} \ \ \ \ \ (37)$

This agrees with the WKB result when ${V_{0}\ll E_{n}^{0}}$, which occurs either when the height of the shelf is very small, or when ${E_{n}^{0}}$ is very large; the latter case occurs when ${n}$ is large and we’re approaching the classical zone.

# Electron gas: a crude model of a solid

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 5.15.

Solids are complex structures, and a full quantum theory of solids is a vast subject. In this post, we’ll have a look at the simplest model of solids, known as the electron gas.

In a solid, typically the outer electrons of an atom become detached from their parent atom and wander more or less freely thoughout the solid. Obviously, the potential felt by such a wandering electron is complex, as it depends on the positions of the various nuclei and on the positions of all the other electrons. A very crude model can be constructed, in which all of these influences on the free electrons are ignored, and the free electrons move in what is effectively a 3-d infinite square well, where the potential is ${V=0}$ everywhere inside the solid and ${V=\infty}$ at the boundaries.

If we consider the 3-d box to be rectangular with side lengths ${l_{x}}$, ${l_{y}}$ and ${l_{z}}$, then the solution is a straightforward generalization of the 1-d case (since the Schrödinger equation separates into 3 equations, one for each dimension), giving a wave function of

$\displaystyle \psi_{n_{x}n_{y}n_{z}}=\sqrt{\frac{8}{l_{x}l_{y}l_{z}}}\sin\left(\frac{n_{x}\pi x}{l_{x}}\right)\sin\left(\frac{n_{y}\pi y}{l_{y}}\right)\sin\left(\frac{n_{z}\pi z}{l_{z}}\right) \ \ \ \ \ (1)$

where the ${n}$s are positive integers. The energies are also a generalization of the 1-d case:

$\displaystyle E_{n_{x}n_{y}n_{z}}=\frac{\pi^{2}\hbar^{2}}{2m}\left(\frac{n_{x}^{2}}{l_{x}^{2}}+\frac{n_{y}^{2}}{l_{y}^{2}}+\frac{n_{z}^{2}}{l_{z}^{2}}\right) \ \ \ \ \ (2)$

For any macroscopic quantity of solid, the number of free electrons is enormous (remember that one mole of a substance contains Avogadro’s number of ${6.02\times10^{23}}$ atoms, and each atom can contribute more than one electron in many cases). Since electrons are fermions, no more than 2 electrons can occupy any single spatial state (one with spin up and the other spin down). Thus the number of quantum numbers (the number of values of the three ${n}$s) required to accommodate all the electrons is going to be huge.

$\displaystyle k^{2}\equiv\frac{\pi^{2}n_{x}^{2}}{l_{x}^{2}}+\frac{\pi^{2}n_{y}^{2}}{l_{y}^{2}}+\frac{\pi^{2}n_{z}^{2}}{l_{z}^{2}}\equiv k_{x}^{2}+k_{y}^{2}+k_{z}^{2} \ \ \ \ \ (3)$

In the ground state, the electrons will fill the individual states starting from ${n_{x}=n_{y}=n_{z}=1}$ and working upwards. We can think of these states as occupying the lattice points in a 3-d rectangular coordinate system with the ${k_{i}}$s as the axes. Each particular combination of the ${k_{i}}$s can be assigned the rectangular cell adjacent to it. Since the distance between lattice points is ${\pi/l_{i}}$ the volume of each cell is ${\pi^{3}/l_{x}l_{y}l_{z}=\pi^{3}/V}$, where ${V}$ is the volume of the solid.

The states will be filled out to a maximum value of ${k_{F}}$ and since ${k_{F}^{2}=k_{x}^{2}+k_{y}^{2}+k_{z}^{2}}$ is the equation of a sphere of radius ${k_{F}}$, and we’re taking all the ${k_{i}}$s to be positive, the filled states in the ground state of the solid occupy the first octant of this sphere. The number of lattice points within this octant is the volume of the octant divided by the volume of each cell, so

$\displaystyle N_{points}=\frac{1}{8}\frac{4\pi}{3}k_{F}^{3}\frac{V}{\pi^{3}}=\frac{Vk_{F}^{3}}{6\pi^{2}} \ \ \ \ \ (4)$

If ${N}$ is the total number of atoms in the solid and ${q}$ is the number of free electrons per atom, the number of spatial states required for these electrons is ${Nq/2}$, so

 $\displaystyle \frac{Nq}{2}$ $\displaystyle =$ $\displaystyle N_{points}\ \ \ \ \ (5)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{Vk_{F}^{3}}{6\pi^{2}}\ \ \ \ \ (6)$ $\displaystyle k_{F}$ $\displaystyle =$ $\displaystyle \left(\frac{3\pi^{2}Nq}{V}\right)^{1/3} \ \ \ \ \ (7)$

If we define ${\rho\equiv Nq/V}$ as the density of free electrons, then we get

$\displaystyle k_{F}=\left(3\pi^{2}\rho\right)^{1/3} \ \ \ \ \ (8)$

with corresponding energy, known as the Fermi energy:

$\displaystyle E_{F}=\frac{\hbar^{2}k_{F}^{2}}{2m}=\frac{\hbar^{2}}{2m}\left(3\pi^{2}\rho\right)^{2/3}=\frac{\hbar^{2}}{2m}\left(3\pi^{2}Nq\right)^{2/3}V^{-2/3} \ \ \ \ \ (9)$

This is the maximum energy that a free electron in the ground state will have.

The total energy can be found from an integral. All states with a given value of ${k}$ have the same energy ${E=\hbar^{2}k^{2}/2m}$, and these states lie in a spherical shell in the first octant. The number of states with this energy in a volume element is the volume element ${k^{2}\sin\theta dkd\theta d\phi}$ divided by the volume of a cell ${\pi^{3}/V}$ times 2 spin states per spatial state, so

 $\displaystyle E_{tot}$ $\displaystyle =$ $\displaystyle \frac{2\hbar^{2}V}{2m\pi^{3}}\int_{0}^{k_{F}}\int_{0}^{\pi/2}\int_{0}^{\pi/2}k^{4}\sin\theta d\phi d\theta dk\ \ \ \ \ (10)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{\hbar^{2}V}{10\pi^{2}m}k_{F}^{5}=\frac{\hbar^{2}\left(3\pi^{2}Nq\right)^{5/3}}{10\pi^{2}m}V^{-2/3} \ \ \ \ \ (11)$

The average energy per electron is then

 $\displaystyle E_{av}$ $\displaystyle =$ $\displaystyle \frac{E_{tot}}{Nq}\ \ \ \ \ (12)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{\hbar^{2}\left(3\pi^{2}\right)^{5/3}\left(Nq\right)^{2/3}}{10\pi^{2}m}V^{-2/3}\ \ \ \ \ (13)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{3}{5}\frac{\hbar^{2}}{2m}\left(3\pi^{2}Nq\right)^{2/3}V^{-2/3}\ \ \ \ \ (14)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{3}{5}E_{F} \ \ \ \ \ (15)$

# Exchange force: infinite square well

Required math: calculus

Required physics: Schrödinger equation

Reference: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Section 5.1.2 & Problem 5.6.

We’ve seen that distinguishable particles and identical particles must be treated differently in quantum mechanics, resulting in different combinations of the single-particle wave functions in 2-particle systems. It’s useful to work out what this means for some of the observables in a 2-particle system.

We can begin by looking at possibly the simplest case: the average positions of the two particles. If the particles are distinguishable, then the wave function is ${\psi(x_{a},x_{b})=\psi_{1}\left(x_{a}\right)\psi_{2}\left(x_{b}\right)}$ and

 $\displaystyle \left\langle x_{a}\right\rangle$ $\displaystyle =$ $\displaystyle \left\langle \psi\left|x_{a}\right|\psi\right\rangle \ \ \ \ \ (1)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left\langle \psi_{1a}\left|x_{a}\right|\psi_{1a}\right\rangle \left\langle \left.\psi_{2b}\right|\psi_{2b}\right\rangle \ \ \ \ \ (2)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left\langle \psi_{1a}\left|x_{a}\right|\psi_{1a}\right\rangle \ \ \ \ \ (3)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left\langle x\right\rangle _{1} \ \ \ \ \ (4)$

where the notation ${\left|\psi_{1a}\right\rangle \equiv\psi_{1}\left(x_{a}\right)}$ and so on.

That is, ${\left\langle x\right\rangle }$ is the mean value of ${x}$ in state ${\psi_{1}}$. We can drop the suffix ${a}$ here, since ${x_{a}}$ is just a dummy name for the integration variable in ${\left\langle \psi_{1a}\left|x_{a}\right|\psi_{1a}\right\rangle }$.

For identical particles,

$\displaystyle \psi_{\pm}\left(\mathbf{r}_{a},\mathbf{r}_{b}\right)=\frac{1}{\sqrt{2}}\left[\psi_{1}\left(\mathbf{r}_{a}\right)\psi_{2}\left(\mathbf{r}_{b}\right)\pm\psi_{2}\left(\mathbf{r}_{a}\right)\psi_{1}\left(\mathbf{r}_{b}\right)\right] \ \ \ \ \ (5)$

This time, working out ${\left\langle x_{a}\right\rangle }$ is a bit messier but not too bad if we use the orthogonality of the two states.

 $\displaystyle 2\left\langle x_{a}\right\rangle$ $\displaystyle =$ $\displaystyle \left\langle \psi_{1a}\left|x_{a}\right|\psi_{1a}\right\rangle \left\langle \left.\psi_{2b}\right|\psi_{2b}\right\rangle +\left\langle \psi_{2a}\left|x_{a}\right|\psi_{2a}\right\rangle \left\langle \left.\psi_{1b}\right|\psi_{1b}\right\rangle \ \ \ \ \ (6)$ $\displaystyle$ $\displaystyle$ $\displaystyle \pm\left\langle \psi_{1a}\left|x_{a}\right|\psi_{2a}\right\rangle \left\langle \left.\psi_{2b}\right|\psi_{1b}\right\rangle \pm\left\langle \psi_{2a}\left|x_{a}\right|\psi_{1a}\right\rangle \left\langle \left.\psi_{1b}\right|\psi_{2b}\right\rangle \ \ \ \ \ (7)$ $\displaystyle \left\langle x_{a}\right\rangle$ $\displaystyle =$ $\displaystyle \frac{1}{2}\left(\left\langle x\right\rangle _{1}+\left\langle x\right\rangle _{2}\right) \ \ \ \ \ (8)$

Thus the mean position of particle ${a}$ is the average of its positions in the two states, which isn’t all that surprising. We’d get the same result for particle ${b}$ of course, since the two particles are identical. This result is true for both bosons and fermions, since the plus/minus terms work out to zero due to the orthogonality of the states ${\psi_{1}}$ and ${\psi_{2}}$.

What is a bit more interesting is the mean square separation of the two particles, that is ${\left\langle \left(x_{a}-x_{b}\right)^{2}\right\rangle }$. This can be worked out using the same procedure as above, and is done by Griffiths in his section 5.1.2, although his notation is a bit different from mine. (I’ve used a numerical suffix on the wave function, since this is the usual notation used for stationary states. Thus a letter suffix indicates which particle and a number suffix indicates which stationary state.) The results are, in my notation, first for distinguishable particles:

$\displaystyle \left\langle \left(x_{a}-x_{b}\right)^{2}\right\rangle =\left\langle x^{2}\right\rangle _{1}+\left\langle x^{2}\right\rangle _{2}-2\left\langle x\right\rangle _{1}\left\langle x\right\rangle _{2} \ \ \ \ \ (9)$

For identical particles, we get

$\displaystyle \left\langle \left(x_{a}-x_{b}\right)^{2}\right\rangle _{\pm}=\left\langle x^{2}\right\rangle _{1}+\left\langle x^{2}\right\rangle _{2}-2\left\langle x\right\rangle _{1}\left\langle x\right\rangle _{2}\mp2\left|\left\langle x\right\rangle _{12}\right|^{2} \ \ \ \ \ (10)$

where the plus (minus) sign on the left and minus (plus) on the right refer to bosons (fermions), and

$\displaystyle \left\langle x\right\rangle _{12}\equiv\left\langle \psi_{1}\left|x_{\;}\right|\psi_{2}\right\rangle \ \ \ \ \ (11)$

In general, then, since the term ${2\left|\left\langle x\right\rangle _{12}\right|^{2}}$ is always positive, bosons tend to be closer together than distinguishable particles while fermions are further apart. This is a sort of pseudo-force which is an entirely quantum mechanical effect of the symmetries of the wave functions. Although it’s not really a force in the sense that electromagnetism and gravity are, it’s known as the exchange force.

As an example, consider 2 particles in the infinite square well. The wave functions for a single particle are

$\displaystyle \psi\left(x\right)=\sqrt{\frac{2}{a}}\sin\frac{n\pi x}{a} \ \ \ \ \ (12)$

where ${a}$ is the width of the well. If the total wave function is a combination of states ${l}$ and ${n}$, then if the particles are distinguishable

 $\displaystyle \left\langle \left(x_{a}-x_{b}\right)^{2}\right\rangle$ $\displaystyle =$ $\displaystyle \left\langle x^{2}\right\rangle _{1}+\left\langle x^{2}\right\rangle _{2}-2\left\langle x\right\rangle _{1}\left\langle x\right\rangle _{2}\ \ \ \ \ (13)$ $\displaystyle$ $\displaystyle =$ $\displaystyle a^{2}\left(\frac{1}{3}-\frac{1}{2l^{2}\pi^{2}}\right)+a^{2}\left(\frac{1}{3}-\frac{1}{2n^{2}\pi^{2}}\right)-2\left(\frac{a}{2}\right)\left(\frac{a}{2}\right)\ \ \ \ \ (14)$ $\displaystyle$ $\displaystyle =$ $\displaystyle a^{2}\left(\frac{1}{6}-\frac{l^{2}+n^{2}}{2(\pi ln)^{2}}\right) \ \ \ \ \ (15)$

In line 2, we used the results of our earlier calculations for the infinite square well.

If the particles are identical, then

 $\displaystyle \left\langle x\right\rangle _{ln}$ $\displaystyle =$ $\displaystyle \left\langle \psi_{l}\left|x_{\;}\right|\psi_{n}\right\rangle \ \ \ \ \ (16)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{2}{a}\int_{0}^{a}\sin\left(\frac{l\pi x}{a}\right)\sin\left(\frac{n\pi x}{a}\right)xdx\ \ \ \ \ (17)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left(-1+\left(-1\right)^{n+l}\right)\frac{4anl}{\left[\pi\left(n^{2}-l^{2}\right)\right]^{2}} \ \ \ \ \ (18)$

This term is zero if ${n+l}$ is even, so there is a difference in the separation only when ${n+l}$ is odd. In general, we have

$\displaystyle \left\langle \left(x_{a}-x_{b}\right)^{2}\right\rangle _{\pm}=a^{2}\left(\frac{1}{6}-\frac{l^{2}+n^{2}}{2(\pi ln)^{2}}\right)\mp2\left[\left(-1+\left(-1\right)^{n+l}\right)\frac{4anl}{\left[\pi\left(n^{2}-l^{2}\right)\right]^{2}}\right]^{2} \ \ \ \ \ (19)$