Sketch 10

Particle in a Box and QHO: Solving PDEs with ODEs, Ansatz, and Algebraic Approach

Huh...

I've been rereading some things about the particle in a box and QHO setups, but when I was redoing a lot of these integral computations myself, I keep running into tiny errors that keep me from getting the right answer.

Ideally, I feel that my work should be readable and simple enough that I shouldn't feel like I need a solutions manual to ensure its correctness. That's the main reason I'm making this sketch, to determine if there is some way that I can tabularize my computation, or determine a handful of things to memorize that will make organization and verification much simpler.

Even if I can't figure out a better way to organize my work, the making of this sketch should at least make me redo my computations, and on top of this, I want to do some experiments simulating these PDEs.

Date Started: July 6, 2024
Date Finished: In Progress...

Edit: I may have gotten a bit distracted from the original purpose of this sketch, :-), but I think the sketch is still pretty good!

The Premise

A while ago, I was trying to solve Problem 2.5 in Griffiths's Quantum Mechanics book (essentially about computing a bunch of things for the infinite square well). The problem is, everything got extremely messy, and now I feel like not even returning to the book, because if I go back to the book, I feel like I have to complete Problem 2.5, but looking at my past work for Problem 2.5, it just looks miserable:

To me, my work isn't thaaat messy (ehh, maybe it is). But it gets the wrong answer (I think), and I can't immediately see where I went wrong. Do I want to erase it and start from scratch? Maybe I should just skip this problem? Ugh. There has to be a better way of doing all this computation with a low error rate that doesn't involve plugging everything into a symbolic computation tool.

In any case, writing it up in LaTeX should do the trick, so that's why I'm making this sketch.

1D Particle in a Box

I'll start with the 1D Particle in the Box. The goal? It's in the name: we want to simulate a free particle in a box. To frame this mathematically, we can create a "box" by creating an appropriate potential energy function: $$V(x)=\{\begin{array}{cc}0& -a\le x\le a,\\ \mathrm{\infty }& \text{else.}\end{array}$$ This is a box because particles have finite energy to begin with, so it could not leave the domain $(-a,a)$ else it would have infinite energy.

The quantum mechanics perspective allows us to determine the "trajectory" of the particle via its wave function $\psi :\mathbb{R}\times [0,\mathrm{\infty })\to \mathbb{C}$, which is a solution to the Schrödinger equation: $$-i\hslash {\partial }_t\psi =-\frac{{\hslash }^2}{2m}{\partial }_{xx}\psi +V\psi$$ Generally, we solve this equation via separation by variables, which gives us so called "separable solutions". These solutions are nice for a few reasons, for which one can check out Griffiths. Most importantly, we can find general solutions by taking linear combinations of these separable solutions. So our aim is to first find the separable solutions.

Separation of Variables

For good measure, let's do the separation of variables. Let $\psi (x,t)=f(x)g(t)$, for $f:\mathbb{R}\to \mathbb{C}$ and $g:\mathbb{R}\to \mathbb{C}$. This notation is a bit different than what I've seen elsewhere, but the notation elsewhere that I've seen tends to be very notationally-abusive, so I prefer to just keep the variables $f$ and $g$ even though it isn't canonical. Then $$i\hslash f{g}^{\prime }=-\frac{{\hslash }^2}{2m}{f}^{″}g+Vfg$$ To make some progress here, we pray that $\psi \ne 0$ and divide both sides by $fg$. Then $$i\hslash \frac{g^{\prime }(t)}{g(t)}=-\frac{{\hslash }^2}{2m}\frac{f^{″}(x)}{f(x)}+V(x)$$ But now notice, the LHS is a function of time, and the RHS of space. Since this holds for all $x$ and $t$, we must have that $$i\hslash \frac{g^{\prime }(t)}{g(t)}=-\frac{{\hslash }^2}{2m}\frac{f^{″}(x)}{f(x)}+V(x)=E$$ for some $E\in \mathbb{C}$. Which gives us two "mini-Schrödinger equations": $$i\hslash {\partial }_{t}g=Eg$$ $$-\frac{{\hslash }^2}{2m}{\partial }_{xx}f=(E-V)f$$ which are affectionately called the time-dependent and time-independent Schrödinger equations respectively (TDSE and TISE). The solution to TDSE can be seen quickly to be: $$g=e^{-\frac{iE}{\hslash }t}$$ so finding these so called "separable solutions" reduces to solving TISE. You may be saying: hold on, what about boundary conditions? Well, in our case, any Dirichlet or Neumann boundary conditions would result in a constant being multiplied out front, so a solution such as $g=C{e}^{-\frac{iE}{\hslash }t}$. But we can always push that constant over to $f$, the solution to TISE, so it sufficies to just take $g=e^{-\frac{iE}{\hslash }t}$ here.

TISE for 1D Particle in a Box

We now return to the particle in a box situation. First, we need to find the separable solutions, which from above just means solving TISE. For $x\in [-a,a]$, we know that $V=0$, so we wish to solve ${\partial }_{xx}f=-\frac{2mE}{{\hslash }^2}f$, whose solutions are well known: sinusoids. $$f=A\cos (kx)+B\sin (kx)$$ where $k=\frac{\sqrt{2mE}}{\hslash }$. For $x\notin [-a,a]$, $V=\mathrm{\infty }$, so the only even slightly reasonable solution is where $f=0$ (technically this is a bit awkward to be multiplying $0$ and $\mathrm{\infty }$, but perhaps we could define it in that way in our extended real number line, or instead work with $V$ as a limit of a sequence of functions).

Physically, it makes the most sense if $f$ is reasonably regular, so we enforce that $f$ should be continuous. Then we must have: $$0=A\cos (-ka)+B\sin (-ka)$$ $$0=A\cos (ka)+B\sin (ka)$$ so $A\cos (-ka)=0$ and $B\sin (ka)=0$. If we further want our solution to be "normalizable" (AKA, we want functions with finite but nonzero $L^2$ norm), we force either $A=0$ and $ka=n\pi$ for $n\in \mathbb{Z}$ or $B=0$ and $ka=\frac{(2n-1)\pi }{2}$ for $n\in \mathbb{Z}$. For a given $n$, we may take $k_n$ to be the corresponding value of $k$ as per these relations. In any case, $2k_{n}a=n\pi$ for $n\in \mathbb{Z}$, and for even $n$, we obtain a solution of the form $B\sin (k_{n}x)$, and for odd $n$, we obtain a solution of the form $A\cos (k_{n}x)$. Up to sign (which doesn't matter because we normalize anyway), then the solutions become $$f_n(x)=\{\begin{array}{cc}C\sin (k_{n}x+k_{n}a)& x\in [-a,a],\\ 0& \text{else.}\end{array}$$ (we add in a $k_{n}a$ to turn the cosine into a sine term by shifting). To summarize, our solution to TISE for particle in a box satisfies a few properties:

• It is continuous, but not continuously differentiable.
• It has finite but nonzero $L^2$ norm for $C\ne 0$ (AKA normalizable).
• It has compact support, and $$\text{supp}(f_n)=\{x:V(x)<\mathrm{\infty }\}$$ (this is why I made the interval on which $V$ was 0 to be $[-a,a]$ as opposed to $(-a,a)$).

Integration

I thought of a few ways to try and be more organized. I couldn't come up with a good way to tabularize all the computation, but what could be an improvement (while not as spacially efficient as tabularization) is indentation. Here's what it might look like for the computations to find $C$ above:

We would like to know several things about the particle in a box. What are the expectations of the position and momentum? Perhaps, can we compute the acceleration of the particle? Additionally, what are the variances? It would be nice to determine a few moments of the probability distributions for positions, momentums, etc. of these particles so that we can get a good feel for them distributionally.

Each computation comes with its own integral:

• $E[x]={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\psi }^{\ast }x\psi dx$
• $E[p]=m{\partial }_t(E[x])$
• $\text{Var}[x]=E[x^2]-(E[x]{)}^2={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\psi }^{\ast }{x}^2\psi dx-(E[x]{)}^2$
• $\text{Var}[p]=E[p^2]-(E[p]{)}^2={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\psi }^{\ast }(\frac{h}{i}\frac{\partial }{\partial x}{)}^2\psi dx-(E[p]{)}^2$

There are, in particular, three main integrals which need to be computed: ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\psi }^{\ast }x\psi dx$, ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\psi }^{\ast }{x}^2\psi dx$, and ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\psi }^{\ast }(\frac{h}{i}\frac{\partial }{\partial x}{)}^2\psi dx$.

For the first two, we are actually very lucky, because things commute, and for separable solutions $\psi =fg$, we know from above that $g{g}^{\ast }=1$. The integrals are thus the same as ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{f}_n^{\ast }{f}_{n}xdx$ and ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{f}_n^{\ast }{f}_n{x}^{2}dx$. Also, it turns out that for our particular problem, we may choose $C$ to be real (in the end it doesn't matter for our general solutions since they will be linear combinations), so that $f_n$ is real-valued, and the integrals become ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{f}_n^{2}xdx$ and ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{f}_n^2{x}^{2}dx$.

Since $f_n^2$ is an even function, $f_n^{2}x$ is odd, so ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{f}_n^{2}xdx=0$.

On the other hand, $f_n^2{x}^2$ is even, so ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{f}_n^2{x}^{2}dx=2{\int }_0^{\mathrm{\infty }}{f}_n^2{x}^{2}dx$. Now we have: $$\begin{array}{rl}\text{(Compact Support)}& {\int }_0^{\mathrm{\infty }}{f}_n^2{x}^{2}dx& ={\int }_0^a{f}_n^2{x}^{2}dx& ={\int }_0^a{C}^2{\sin }^2(k_{n}x+k_{n}a)\cdot x^{2}dx\\ & =C^2{\int }_0^a{x}^2{\sin }^2(k_{n}x+k_{n}a)dx\\ \text{(Double Angle Formula)}& & =C^2{\int }_0^a{x}^2(1-2\cos (2k_{n}x+2k_{n}a))dx& =C^2{\int }_0^a{x}^{2}dx-2C^2{\int }_0^a{x}^2\cos (2k_{n}x+2k_{n}a)dx\end{array}$$ Now the first expression becomes $C^2{a}^3/3$. As for the second expression, we do integration by parts tabularly.

Quantum Harmonic Oscillators (QHO)

The importance of the QHO is that in nature, whenever a system is subject to a conservative force, and is around a stable equilibrium, the fluctuation around that equilibrium is approximately a harmonic oscillator. Mathematically, it is also distinct from our previous example because now $V$ isn't constant with respect to $x$. For a harmonic motion, our angular frequency satisfies $\omega =\sqrt{\frac{k}{m}}$, so $k=m{\omega }^2$. This means our Shrodinger equation becomes: $$-\frac{{\hslash }^2}{2m}\frac{d^2{\psi }_E}{d{x}^2}+\frac{1}{2}m{\omega }^2{x}^2{\psi }_E=E{\psi }_E$$ Note that this section of this sketch closely follows the methods from two of my Quantum Mechanics professor's lectures.

Analytic Method

First we will use an analytic method, shown by my quantum mechanics professor. There are four steps:

1. Introduce dimensionless parameters (we write $x=b\xi$, where $b$ has units of lengths and $\xi$ is a dimensionless parameter).
2. Try to guess the ansatz by considering some limits ($\xi \to \pm \mathrm{\infty }$ and $\xi \to 0$).
3. Substitute the ansatz into the TISE and consider the power series.
4. Impose boundary condition to find the allowed energy values.

From the first step, our equation becomes: $$-\frac{{\hslash }^2}{2m{b}^2}\frac{d^2{\psi }_E}{d{\xi }^2}+\frac{1}{2}m{\omega }^2{b}^2{\xi }^2{\psi }_E=E{\psi }_E$$ which rewriting becomes: $$\frac{d^2{\psi }_E}{d{\xi }^2}-\frac{m^2{\omega }^2{b}^4}{{\hslash }^2}{\xi }^2{\psi }_E+\frac{2mE{b}^2}{{\hslash }^2}{\psi }_E=0$$ We now notice that each term should have the same units as ${\psi }_E$, so the coefficients of the second and third terms should be dimensionless. Thus, we can take $b=\sqrt{\frac{\hslash }{m\omega }}$ which defines a length scale based on the fact that $\frac{m^2{\omega }^2{b}^4}{{\hslash }^2}$ should be dimensionless. Then, $\frac{2mE{b}^2}{{\hslash }^2}=\frac{2E}{\hslash \omega }$, which is idneed dimensionless. It is natural here to define $\mathcal{E}=\frac{2E}{\hslash \omega }$. With these definitions, TISE becomes: $$\frac{d^2{\psi }_E}{d{\xi }^2}-{\xi }^2{\psi }_E+\mathcal{E}{\psi }_E=0$$

Now it's time for the second step, guessing an ansatz.
Consider $\xi \to \mathrm{\infty }$, or $\xi >>\mathcal{E}$. Then TISE becomes $\frac{d^2{\psi }_E}{d{\xi }^2}-{\xi }^2{\psi }_E=0$. This is sorta like $\frac{d^2{\psi }_E}{d{\xi }^2}-k^2{\psi }_E=0$, which makes us guess ${\psi }_E={\xi }^m{e}^{\pm {\xi }^2/2}$. Then $${\psi }_E^{\prime }=m{\xi }^{m-1}{e}^{\pm {\xi }^2/2}\pm {\xi }^{m+1}{e}^{\pm {\xi }^2/2}$$ $${\psi }_E^{″}={\xi }^{m+2}{e}^{\pm {\xi }^2/2}[1\pm (2m+1)\frac{1}{{\xi }^2}+m(m-1)\frac{1}{{\xi }^2}]$$ Now in the limit, this indeed does satisfy $\frac{d^2{\psi }_E}{d{\xi }^2}-{\xi }^2{\psi }_E=0$, so this ansatz works. Since TISE is linear, this would mean $\sum _m{A}_m{\xi }^m{e}^{\pm {\xi }^2/2}$ is also a solution.
Consider $\xi \to 0$. This becomes $\frac{d^2\psi }{d{\xi }^2}+\mathcal{E}{\psi }_E=0$. The general solution is then ${\psi }_E\sim A\sin (\sqrt{\mathcal{E}}\xi )+B\cos (\sqrt{\mathcal{E}}\xi )\approx A\sqrt{\mathcal{E}}\xi +B+\frac{B}{2}\mathcal{E}{\xi }^2$ via Taylor expansion. Considering our ansatz from the previous limit, it seems that under a taylor expansion, that would still work here! The $\sum _m{A}_m{\xi }^m$ part is somewhat like the $A\sqrt{\mathcal{E}}\xi$ term and the exponential term is somewhat like the other term.
Based on these two limits, we guess a form ${\psi }_E=h(\xi )e^{-{\xi }^2/2}$ for the solution where $h$ has a power series (we want the minus sign so that ${\psi }_E$ is finite as $\xi$ blows up). You might be saying: if we were going to use a power series, why do the ansatz step? If that's what you're thinking, please see Appendix A for a more thorough discussion of that.

So far so good, for the third step, we will now substitute the ansatz into the TISE. Notice that ${\psi }_E^{\prime }=(h^{\prime }-\xi h)e^{-{\xi }^2/2}$ and ${\psi }_E^{″}=(h^{″}-h-\xi h^{\prime }-\xi (h^{\prime }-\xi h))e^{-{\xi }^2/2}=(h^{″}+({\xi }^2-1)h-2\xi h^{\prime })e^{-{\xi }^2/2}$. The Shrodinger equation then becomes $$h^{″}-2\xi h^{\prime }+({\xi }^2-1+\mathcal{E}-{\xi }^2)h=0$$ Now using that $h$ is analytic, we write $h=\sum _{j=0}^{\mathrm{\infty }}{a}_j{\xi }^j$, so $$\sum _{j=0}^{\mathrm{\infty }}{a}_{j+2}(j+2)(j+1){\xi }^j-2\xi \sum _{j=1}^{\mathrm{\infty }}{a}_{j}j{\xi }^{j-1}+(\mathcal{E}-1)\sum _{j=0}^{\mathrm{\infty }}{a}_j{\xi }^j=0$$ $$2a_{j+2}+a_0(\mathcal{E}-1)+\sum _{j=1}^{\mathrm{\infty }}(a_{j+2}(j+2)(j+1)-2j{a}_j+(\mathcal{E}-1)a_j){\xi }^j=0$$ From this we get the recursive relation: $$a_{j+2}=\frac{(2j+1-\mathcal{E})}{(j+2)(j+1)}{a}_j$$ Basically, once we choose $a_0$ and $a_1$, we can find the result. In order for ${\psi }_E\to 0$ to go to zero as $|\xi |\to \mathrm{\infty }$, we need the series to terminate at a finite term, so we need $2j+1-\mathcal{E}=0$ for some $j=n$. For a demonstration of this result, see Appendix B. Since $\mathcal{E}=\frac{2E}{\hslash \omega }$, $$\boxed{{\displaystyle E=(n+\frac{1}{2})\hslash \omega }}$$ which is called energy quantization.

We can only have bound states because $V(x)$ tends to infinity as $|x|\to \mathrm{\infty }$. Since $V(-x)=V(x)$ so $V$ is even, the bound states are either symmetric or antisymmetric. In the symmetric case, we have $a_1=a_3=\cdots =0$. Now for each value of $n$, we obtain a polynomial $h$, and $a_0$ is an appropriate normalization factor. Notice that $n=0$ and $n=1$ lead to the same polynomial $h$, and similarly with $n=2k$ and $n=2k+1$ in general. In the antisymmetric case, we construct $a_3$, $a_5$, and so on from $a_1$. The $h$ polynomials here are called Hermite polynomials. We then get ${\psi }_n$ for each $n$ being the product of a Hermite polynomial and a Gaussian.

Algebraic Method

Notice that clasically, $H=\frac{p^2}{2m}+\frac{1}{2}m{\omega }^2{x}^2=\frac{1}{2m}(p+im\omega x)(p-im\omega x)$. This motivates trying a similar factorization of the operator in quantum mechanics. We can try to write our Hamiltonian operator in the form: $$\hat{H}\sim {\hat{a}}^{†}\hat{a}$$ (notice that the RHS is always Hermitian for all $\hat{a}$). We can squint at this and write: $$\hat{a}\equiv \frac{1}{\sqrt{2m\hslash \omega }}(m\omega \hat{X}+i\hat{P}).$$ The reason for the $\hslash \omega$ in the denominator is because of the "energy scale" (recall the stuff we did before with dimensionless quantiites). $${\hat{a}}^{†}\equiv \frac{1}{\sqrt{2m\hslash \omega }}(m\omega \hat{X}-i\hat{P}).$$ These are called the annihilation and creation operators respectively. You can compute quickly via linearity and basic properties that the commutation relation between them is $[\hat{a},{\hat{a}}^{†}]=\mathbb{1}$, the identity! So they don't commute. One can compute fairly easily that: $$({\hat{a}}^{†}\hat{a}+\frac{1}{2})\hslash \omega =\hat{H}.$$ from which it becomes not too difficult to compute that $[\hat{a},\hat{H}]=\hslash \omega \hat{a}$ and $[{\hat{a}}^{†},\hat{H}]=([\hat{H},\hat{a}]{)}^{†}=(-[\hat{a},\hat{H}]{)}^{†}=-\hslash \omega {\hat{a}}^{†}$. From this we get two nice results: (assuming they aren't the zero vector) $\hat{a}|E⟩$ is an energy eigenstate of energy $E-\hslash \omega$, and ${\hat{a}}^{†}|E⟩$ is an energy eigenstate of energy $E+\hslash \omega$. First let's examine the usefulness of these results, and then we will prove them later in Appendix C. Once we know these relations, we can define a a ground state $|0⟩$ as an eigenvalue of $\hat{a}$, and notice that it must have energy $\frac{1}{2}\hslash \omega$. Then repeatedly applying the creation operator gets us more and more states of increasing energies. Let $|n⟩$ denote $|E_n⟩$, where $E_n=(n+\frac{1}{2})\hslash \omega$ and $n\ge 0$. Then $$\begin{array}{rl}|c_n{|}^2& =|c_n{|}^2⟨n-1|n-1⟩\\ & =⟨n|{\hat{a}}^{†}\hat{a}|n⟩\\ & =⟨n|(\frac{\hat{H}}{\hslash \omega }-\frac{1}{2})|n⟩\\ & =⟨n|\frac{(n+\frac{1}{2})\hslash \omega }{\hslash \omega }|n⟩-\frac{1}{2}\\ & =n+\frac{1}{2}-\frac{1}{2}\\ & =n\end{array}$$ so $|c_n{|}^2=n$. In other words, $\hat{a}|n⟩=\sqrt{n}|n-1⟩$. Substituting this into $⟨n|{\hat{a}}^{†}\hat{a}|n⟩$, we see that ${\hat{a}}^{†}|n⟩=\sqrt{n+1}|n+1⟩$. Sometimes we define the number operator $\hat{N}={\hat{a}}^{†}\hat{a}$, so that $\hat{H}=\hslash \omega (\hat{N}+\frac{1}{2})$. Aside from being mathematically nice, these are very useful operators when solving problems. There are two more important ideas: $$\hat{X}=\frac{\sqrt{2m\hslash \omega }}{2m\omega }({\hat{a}}^{†}+\hat{a})$$ $$\hat{P}=\frac{i\sqrt{2m\hslash \omega }}{2}({\hat{a}}^{†}-\hat{a})$$ which allows us to much more easily do computations relating to $\hat{X}$ and $\hat{P}$ when working with energy eigenstates.

Now we can actually solve this via an algebraic method. By definition, $\hat{a}|0⟩=0$. But this means, by acting by $⟨x|$ on the left, $$\frac{1}{\sqrt{2m\hslash \omega }}[m\omega x+i(-i\hslash \frac{d}{dx})]{\psi }_0=0$$ Simplifying, this becomes $$m\omega x{\psi }_0+\hslash \frac{d{\psi }_0}{dx}=0$$ Separating the variables, $$\int -\frac{m\omega x}{\hslash }dx=\int \frac{1}{{\psi }_0}d{\psi }_0$$ so ${\psi }_0=A{e}^{-m\omega x^2/2\hslash }$. Since $x=b\xi$ with $b=\sqrt{\frac{\hslash }{m\omega }}$, ${\psi }_0=A{e}^{-{\xi }^2/2}$. So this is arguably a lot easier than doing all the weird analytic stuff. :-)

Appendix A: What if we skip the ansatz step?

This time, let's entirely skip the ansatz step and try to solve the differential equation by Taylor expanding ${\psi }_E$.

Appendix B: Energy Quantization of QHO

Appendix C: PRoving Important Algebraic Relations

References