Appendix B - Spin States of the Hydrogen Atom

Index

1. Introduction

The ferromagnetic alignments discussed in chapter 2 of Schroeder's Thermal Physics are a simplification on what is actually some fundamental quantum mechanics. While what is written below is not strictly necessary for understanding the sheer classical thermodynamics in section 1, this subject is important enough on its own to deserve further discussion.

Pure, classical paramagnets are systems with particles aligned randomly in the absence of an external field; only becoming alinged after some field is applied to them. We tend to think of the systems as non-self-interacting. In this picture, the alignment of the dipoles are due to a torque being applied to them by the external field — the dipoles lining up parallel to the field lines. The quantum picture of this alignment is closer to the truth, and requires greater detail. To understand, we first need to grasp the concepts of orbital angular momentum and spin in a quantum mechanical setting. We focus on what these terms mean in the setting of the hydrogen atom, but much of what is discussed is adaptable to other settings, including ferromagnetic materials. It is important to distinguish what we mean by each of these terms, if we are to understand them individually we must understand the differences between them.

Sans much detail, what we call "spin" in quantum mechanics is a property of quantum objects. Spin should be understood as an intrinsic angular momentum fundamentally ascribed to a particle, unrelated to any rotation of the object. It can be shown to exist observationally, by recalling the Stern-Gerlach experiment. Crucially, it is also be seen to be a quantized value: the deflection of silver atoms in the inhomogeneous magnetic field causes the particles to cluster in two points. Classically, one would expect a uniform, continuous distribution — instead, two distinct concentrations of particles are seen. The beam being split in two necessitates that this intrinsic ground-state spin could not be quantized as an integer. If it were, the beam would be split into 3 parts: -1, 0, and 1. It is concluded that the particles in the beam had an intrinsic angular momentum of $\frac{1}{2}$, which is where we get the term "spin-$\frac{1}{2}$". All known fermions have half-integer spins, including the proton and electron that make up the hydrogen atom. After discussing in detail the physics of the hydrogen atom, we will proceed with a more detailed overview of spin that justifies our claims thus far. As is so often the case with a seemingly simple question, the best answers are quite complicated.

The orbital angular momentum is a familiar concept, it is in many ways similar to the angular momentum of bodies in classical physics. The difference is that in the case of quantum objects, the orbital angular momentum is "quantized". It takes on certain, discrete values. In order to understand these values, we must look back to the Schrodinger equation.

2. The Schrodinger Equation for Hydrogen

We ask what the Schrodinger equation predicts for the stationary orbits of the electron in the hydrogen atom. We get the following: $$ \frac{-\hbar^2}{2\mu}\nabla^2 \psi - \frac{e^2}{4\pi \epsilon_0 r}\psi = E\psi, $$ remembering that the potential energy for the orbit of the electron is $\frac{-e^2}{4\pi \epsilon_0r}$, and $\mu =\frac{m_em_p}{m_e+m_p}$ is the center of mass between the electron and the proton. We choose to work in spherical coordinates. Now, try not to get discouraged: what I am about to write is the meat and potatoes of what is learned in a second course on quantum theory at the undergraduate level. I will break it down into parts for easier understanding. We take a deep breath and write the Schrodinger equation for the hydrogen atom in spherical coordinates now. $$ \frac{-\hbar^2}{2 \mu} \left[\frac{1}{r^2}\frac{\partial}{\partial r} \left(r^2\frac{\partial \psi}{\partial r}\right) + \frac{1}{r^2\sin \theta}\frac{\partial}{\partial \theta} \left(\sin \theta \frac{\partial \psi}{\partial \theta} \right) + \frac{1}{r^2\sin^2 \theta} \frac{\partial^2 \phi}{\partial \phi^2} \right] - \frac{e^2}{4\pi \epsilon_0 r}\psi = E\psi. $$

From here we use a classic trick in partial differential equations called "seperation of variables". The notion is to seperate the PDE into seperate ODEs with independent variables. In our case, we assume we can break our solution into a radial component, $R(r)$, and an angular component, $Y(\theta,\phi)$. That is, we write (up to a normalization) $$ \psi(r,\theta,\phi) = R(r) \cdot Y(\theta,\phi). $$ This means we can write $$ \frac{\partial \psi}{\partial r} = Y\frac{dR}{dr}, $$ and so on for the angular partial derivatives. With this in mind, we're able to re-arrange the Schrodinger equation: $$ \frac{Y}{r^2}\frac{d}{dr}\left(r^2\frac{dR}{dr}\right)+ \frac{R}{r^2 \sin\theta}\frac{\partial}{\partial \theta} \left(\sin \theta \frac{\partial Y}{\partial \theta}\right) + \frac{R}{r^2 \sin^2 \theta} \frac{\partial^2Y}{\partial \phi^2} + \frac{2\mu}{\hbar^2} \left(-E + \frac{e^2}{4 \pi \epsilon_0 r}\right)RY =0. $$ Simplifying further, we can multiply by $r^2$ and divide by $RY$. This gives us a well-seperated equation: $$ \frac{1}{R}\frac{d}{dr}\left(r^2\frac{dR}{dr}\right)+ \frac{R}{Y \sin\theta}\frac{\partial}{\partial \theta} \left(\sin \theta \frac{\partial Y}{\partial \theta}\right) + \frac{1}{Y \sin^2 \theta} \frac{\partial^2Y}{\partial \phi^2}+\frac{2\mu r^2}{\hbar^2} \left(-E + \frac{e^2}{4 \pi \epsilon_0 r}\right) =0. $$ The first and fourth terms are radial, and the second and third terms are angular. They are only able to balance if they equate to the same constant, with opposite sign. We call that constant $C$, and write a radial and angular equation: $$ \begin{align} & \frac{d}{dr}\left(r^2\frac{dR}{dr}\right) + \frac{2 \mu r^2}{\hbar^2} \left(-E + \frac{e^2}{4 \pi \epsilon_0 r}\right)R - CR = 0,\\ \\ & \frac{1}{\sin \theta}\frac{\partial}{\partial \theta} \left( \sin \theta \frac{\partial Y}{\partial \theta}\right) + \frac{1}{\sin^2 \theta}\frac{\partial^2Y}{\partial \phi^2} +CY = 0. \end{align} $$

We're not close to the end. In fact, there is much heavy lifting yet to do. In particular, at this stage we need to seperate the angular equation into its constituent parts, and solve those portions. Then we need to solve the radial equation. We start with the angular components. Assuming that $Y(\theta, \phi) = \Theta(\theta) \cdot \Phi(\phi)$, we write $$ \frac{\Phi}{\sin \theta}\frac{d}{d\theta} \left(\sin\theta\frac{d\Theta}{d\theta}\right) + \frac{\Theta}{\sin^2\theta}\frac{d^2\Phi}{d\phi^2} + C\Phi\Theta = 0. $$ In a similar way as we did with the Schrodinger equation, we again multiply by $\sin^2 \theta$ and divide by $\Phi\Theta$. Our new equation is $$ \frac{\sin\theta}{\Theta}\frac{d}{d\theta} \left(\sin\theta\frac{d\Theta}{d\theta}\right) + \frac{1}{\Phi}\frac{d^2\Phi}{d\phi^2} + C\sin^2\theta = 0. $$ Taking a new seperation constant $B$, we have a polar equation and an azimuthal equation: $$ \begin{align} & \frac{\sin\theta}{\Theta}\frac{d}{d\theta} \left(\sin\theta\frac{d\Theta}{d\theta}\right) + C\sin^2\theta - B= 0,\\ \\ & \frac{1}{\Phi}\frac{d^2\Phi}{d\phi^2} + B = 0. \end{align} $$

3. The Azimuthal Equation

Now we can begin solving. I'm going to start with the azimuthal portion, because it's the simplest. As anyone who took first year physics should be able to recognize, this is a 2nd order linear ODE with constant coefficients — it has a solution $$ \Phi(\phi) = a_1e^{im\phi}+a_2e^{-im\phi} $$ with $B = m^2$. We may choose a coordinate system in which $a_2 = 0$, and do so for simplicity. Something else of interest: we see immediately from this that $m$ has to be integer-valued: if it were not, there would be different values for $\theta = 0$ and $\theta = 2\pi$. We dust off our hands and say the solution to the azimuthal equation is $$ \Phi(\phi) = a_1e^{im\phi}. $$ Next up is the much more involved polar equation. Deep breaths, we'll get through it.

4. The Polar Equation

Alright. Taking $B=m^2$ again, the polar equation is $$ \frac{\sin\theta}{\Theta}\frac{d}{d\theta} \left(\sin\theta\frac{d\Theta}{d\theta}\right) + C\sin^2\theta - m^2= 0. $$ Rearranging, $$ \frac{1}{\sin\theta}\frac{d}{d\theta} \left(\sin\theta\frac{d\Theta}{d\theta}\right) + \left(C-\frac{m^2}{\sin^2\theta}\right)\Theta = 0. $$ Now we do something very clever. We make a substitution: $P(\cos\theta) \equiv \Theta(\theta) , x \equiv \cos\theta$. Our differential changes, now, to $$ \frac{d}{d\theta} = \frac{dx}{d\theta}\frac{d}{dx} = -\sin\theta\frac{d}{dx} $$ and $$ \frac{1}{\sin\theta}(-\sin\theta)\frac{d}{dx} \left(\sin\theta(-\sin\theta\frac{dP}{dx}\right)+ \left(C - \frac{m^2}{\sin^2\theta}\right)P = 0. $$ We can make this cleaner. This rearranges to $$ \frac{d}{dx}\left(\sin^2\theta\frac{dP}{dx}\right)+ \left(C-\frac{m^2}{\sin^2\theta}\right)P =0. $$ Making the further substitution that $\sin^2t = 1- \cos^2t = 1-x^2$, we get the following: $$ \begin{align} & \frac{d}{dx}\left((1-x^2)\frac{dP}{dx}\right)+ \left(C-\frac{m^2}{1-x^2}\right)P =0\\ \\ & = (1-x^2)\frac{d^2P}{dx^2} - 2x\frac{dP}{dx} + \left(C - \frac{m^2}{1-x^2}\right)P = 0. \end{align} $$ Here's where things get a little ugly. This is still a second-order ODE, but it no longer has constant coefficients. The solutions are called "associated Legendre polynomials", which are hairy enough that it makes me want to write an appendix to the appendix. I'm gonna do it here, but you can skip to the solution if you're not up to it. Here goes nothing.

4.1 The Associated Legendre Polynomial

We start with the case that $m=0$. Our ODE becomes $$ (1-x^2)\frac{d^2P}{dx^2} - 2x\frac{dP}{dx} + CP = 0. $$ We notice that the coefficient on the 2nd term is the derivative of the coefficient on the first term. This leads us to suspect that our kernel is probably $(1-x^2)$. The multiple derivatives indicate we should make a substitution along the lines of $w = (1-x^2)^n$ , $w' = -2xn(1-x^2)^{n-1}$ which gives us the following relation: $$ w'(1-x^2) +2nxw =0. $$ However, repeated derivatives of this function don't give us the form of equation we're looking for. We Take our original equation and flip its sign: $$ (x^2-1)\frac{d^2P}{dx^2} + 2x\frac{dP}{dx} - CP = 0, $$ and we can use a similar substitution to find that, letting $\omega = (x^2-1)^l$, $$ \omega'(x^2-1) - 2nx\omega =0. $$ The $l+1$th derivative of this statement gives you $$ (x^2-1)\omega^{(l+2)}+2x\omega^{(l+1)}-l(l+1)\omega^{(l)}, $$ which is in the form we're looking for. This means we may write, noting that $C = l(l+1)$ now, that $$ P = \frac{d^l}{dx^l}(x^2-1), $$ which is the Legendre polynomial. This polynomial may be rewritten as: $$ P_l =\frac{1}{2^l l!} \frac{d^l}{dx^l}(x^2-1)^l. $$ This is the "Rodrigues formula". And it's all well and good, but it's still not what we're after — we don't generally have $m=0$. Let's explore the non-zero case.

Unfortunately, the same solution does not work. Our new formulation would look something like $$ (1-x^2)z''-2x(k+1)z'+(l(l+1)-k(k+1))z=0, $$ if we took additional derivatives to power $k$ and replaced things as we had in the previous paragraph. If, instead, we recast the original equation (with non-zero $m$), we can get something of the form $$ (1-x^2)^2y''-2x(k+1)y'+(l(l+1)(1-x^2)-m^2)z=0. $$ This indicates that the final function should be of the form $$ \xi = (1-x^2)^a\frac{d^b}{dx^b}y_n. $$ Immediately, we see that $b$ has to be a positive integer — no Feynman half-derivative tricks allowed here. Taking derivatives appropriately, we can see that if we force $y=x$ and let $b$ be 1, $$ \begin{align} &(1-x^2)\xi'' -2x\xi' + \left(l(l+1)-\frac{m^2}{1-x^2}\right) \xi\\ & = -oa(1-x^2)^a+a(a-1)(4x^2)(1-x^2)^{a-1}\\ &\ + yx^2a(1-x^2)^{a-1}+l(l+1)(1-x^2)^{a-1}\\ &\ - m^2(1-x^2)^{a-1} = 0. \end{align} $$ For this whole thing to be zero, we need the coefficients to go to zero. This means we can coefficient-match our terms to solve for $a$. Doing so, $$ a(a-1)4x^2(1-x^2)^{a-1}+4x^2a(1-x^2)^{a-1}-m^2(1-x^2)^{a-1} = 0, $$ which yeilds that the relationship between coefficients of x is $$ 4a^2 = m^2 \rightarrow a = \frac{|m|}{2}. $$ Now we can finally put this all together, and the relationship:

$$ P_l^m(x) = (1-x^2)^{\frac{|m|}{2}}\left(\frac{d}{dx} \right)^{|m|}\left(\frac{1}{2^l l!}\right)\left(\frac{d}{dx} \right)^l(x^2-1)^l, $$ or, $$ P_l^m(x) = (1-x^2)^{\frac{|m|}{2}} \left(\frac{d}{dx}\right)^{|m|}P_l(x). $$ The famous "associated Legendre function", our solution to the polar equation. Notice the absolute values on $m$ here, this will have consequences for our allowable quantum numbers that will be discussed below. We often include in our solution a factor of $(-1)^m$, though this is a matter of convention as it is sometimes included in the definition of the spherical harmonics.

The last step (this is truly the most painful part) is for us to find a constant that can normalize the angular portion of our solution. We say that the non-radial solution thus far is $$ Y_l^m(\theta,\phi) = Ae^{im\phi}P_l^m\cos\theta. $$ Using the Born rule, and noting that the radial and angular portions of the wave function may be normalized seperately, we see that $$ \begin{align} & \int_0^{2\pi}\int_0^{\pi}\|Y_l^m\|^2\sin\theta d\theta d\phi =1,\\ \\ & 2\pi A^2\int_{-1}^{1} \left(P_l^m(\cos\theta))\right)^2 \sin\theta d \theta =1. \end{align} $$ We may substitute with our Rodrigues formula and with the definition of the associated Legendre polynomials, as well as setting $(x^2-1)=X$, to recover: $$ 2\pi A^2 \frac{(-1)^m}{2^{p+q}p!q!} \int_{-1}^1X^m \frac{d^m}{dx^m} \left(\frac{d}{dx}\right)^pX^p \frac{d^m}{dx^m} \left(\frac{d^q}{dx^q}X^q\right)dx =1, $$ which looks a little daunting. We can break this down. We use integration by parts $m+q$ times to get $$ 2\pi A^2 \frac{(-1)^{2m+q}}{2^{p+qp!q!}}\int_{-1}^1 \frac{d^{m+q}}{dx^{x+q}}\left(X^m\frac{d^{m+p}}{dx^{x+p}}X^p \right)X^qdx = 1. $$ Then, the Leibniz formula allows us to express our integrand as $$ X^q \frac{d^{m+q}}{dx^{x+q}}\left(X^m\frac{d^{m+p}}{dx^{x+p}}X^p \right) = X^q \sum_{i=0}^{m+q} \frac{(m+q)!}{(m+q-i)!} \frac{d^{m+q-i}}{dx^{m+q-i}}X^m\frac{d^{m+p+i}}{dx^{m+p+i}} X^p $$ First, we can see that $X^m$ has no power of $x$ greater than $x^{2m}$, otherwise the integrand would go to zero. Same goes for $X^p$ having no powers of $x$ greater than $x^{2p}$. Jointly, this means that $m+q-i \leq 2m$ and $m+p+i \leq 2p$, or $i \leq p-m, i \geq q-m$. This implies that $q=p$, becuase otherwise the integral vanishes. We can then write $i=q-m$, and proceed with that substitution: $$ \int_{-1}^1 (P_q^m(x))dx = \frac{(-1)^{2m+q}(q+m)!(2q)!} {2^{2q}(q!)^2(q-m)!(2m)!}\int_{-1}^1X^q\left(\frac{d^{2m}} {dx^{2m}}X^m\right)\left(\frac{d^{2q}}{dx^{2q}}X^q\right)dx. $$ Noting that $\frac{d^{2m}}{dx^{2m}}X^m = (2m)!$, this is the same as writing $$ \int_{-1}^1 (P_q^m(x))dx = \frac{(-1)^{2m+q}(q+m)!(2q)!} {2^{2q}(q!)^2(q-m)!}\int_{-1}^1X^q dq. $$

Now for something a little less straightforward. We re-substitute $x = \cos \theta$ and write our integral in terms of sine. Then we express the integrand as a beta function, and further as gamma functions that will leave us with a fairly clean factorial term as our result. It should be noted that "seeing this" is not really expected. The way I was taught, there wasn't this much detail — we were supposed to just accept that the normalization factor presented was just correct and move on. I would never expect someone to calculate this result by hand and expect the physicists who came up with this model, most notably Bohr, used an integral table to get the answer that they did. Anyway, for the curious I have done the calculation below: $$ \begin{align} & \int_{-1}^1X^qdx = (-1)^q2\int_0^{\frac{\pi}{2}} \sin^{2q+1} \theta d\theta = (-1)^{q+1}\int_0^{\frac{\pi}{2}} \sin^{2(\frac{1}{2}(2q+2))-1}\theta \cos^{2(\frac{1}{2})-1} \theta d\theta\\ \\ & = (-1)^q\frac{1}{2}\left(B\left( 2q+1,\frac{1}{2}\right)\right) = (-1)^q\frac{1}{2}\frac{\Gamma\left(\frac{2q+2}{2}\right) \Gamma\left(\frac{1}{2}\right)}{\Gamma\left(\frac{2q+1}{2}+1 \right)} = (-1)^q \frac{\sqrt{\pi}}{2}\frac{\Gamma(q+1} {\Gamma \left(q+\frac{3}{2}\right)}. \end{align} $$ For the last term's denominator, we invoke the Legendre duplication formula: $\Gamma(z)\Gamma(z+\frac{1}{2}) = 2^{1-2z}\sqrt{\pi}\Gamma(2z).$ Writing $\Gamma(q+\frac{3}{2} )$ as $\Gamma((q+1)+\frac{1}{2})$, we rearrange and pass to the Legendre formula. Plugging it all in, we get $$ \int_{-1}^1X^qdx = \frac{(-1)^q(q!)^2(2^{2q+1}}{(2q+1)!}. $$ Now we plug this back into our formula above, and we get that $$ \begin{align} & \left(\frac{(-1)^{2m+q}(q+m)!(2q)!}{2^{2q}(q!)^2(q-m)!} \right)\left(\frac{(-1)^q(q!)^2(2^{2q+1})}{(2q+1)!}\right)\\ \\ & =(-1)^{2m+2q}\frac{2(q+m)!}{(2q+1)(q-m)!}= 1. \end{align} $$ So letting $q=l$ and $\alpha = (-1)^{l+m}$, we write $$ \begin{align} & A = \sqrt{\frac{(2l+1)(l-m)!}{2(l+m)!}},\\ \\ & Y_l^m(\theta , \phi) = \alpha \sqrt{\frac{(2l+1)(l-|m|)!}{4 \pi(l+|m|)!}}e^{im\phi} P_l^m(\cos \theta). \end{align} $$ At this point, I encourage the reader who has worked through this with me to take a sigh of relief and maybe a break, because that's what I did.

5. The Radial Equation

Ok, so with the angular part squared away we can finish our description of the hydrogen atom by examining the radial equation, which is $$ \frac{d}{dr}\left(r^2\frac{dR}{dr}\right) +\frac{2\mu r^2}{\hbar^2}\left(-E + \frac{e^2}{4\pi \epsilon_0 r} - \frac{\hbar^2l(l+1)}{2 \mu r^2}\right)R =0. $$ While it is possible to proceed with this, we are also able to rewrite it in the followisng way, which I believe is more useful: setting $4\pi \epsilon_0 =1$ and realising that it is possible to rewrite the first term as $r \frac{d^2}{dr^2}(rR)$, $$ r\frac{d^2}{dr^2}(rR) +\frac{2\mu r^2}{\hbar^2}\left(-|E| + \frac{e^2}{r} - \frac{\hbar^2l(l+1)}{2 \mu r^2}\right)R =0. $$ We're writing an absolute value around $E$ because we're actually interested in the negative energy eigenstates, as they correspond to bounded (rather than scattering) solutions. We make a load of substitutions here to simplify our equation. I'll introduce them a few at a time rather than all at once. If we take $\rho = 2\kappa r$, $d\rho = 2\kappa dr$, $r^2 = \frac{\rho^2}{4\kappa^2}$, and $dr^2 = \frac{1}{4\kappa^2}d\rho^2$, the equation above becomes $$ \rho \frac{d^2}{d\rho^2}(\rho R) +\frac{\mu \rho}{2\kappa^2 \hbar^2}\left(-|E| + \frac{e^2}{\rho} - \frac{2\kappa^2 \hbar^2l(l+1)}{2 \mu \rho^2}\right)\rho R =0. $$ We divide by a factor of $\rho$ and make the replacements $|E| = \frac{\hbar^2 \kappa^2}{2\mu}$, $a_0 = \frac{\hbar^2}{\mu e^2}$, $\lambda = \frac{1}{\kappa a_0}$, $u(\rho) = \rho R$ to get $$ \left(\frac{d^2}{d\rho^2} - \frac{l(l+1)}{\rho^2} + \frac{\lambda}{\rho} -\frac{1}{4} \right)u =0. $$ The general solution to this kind of equation comes in the form of a wave function, which we will call $j$, which is described by its singular values, and a polynomial, $P$, which we will discuss more later. That is, we write $u = \left(\prod_ij_i\{s\}\right)P(r)$, where the set $\{s\}$ are the singular points particular to $j_i$. We need singular points at the center of the atom, and at $n \rightarrow \infty$. We cover the form of $j$ in those cases. First, the infinite distance case — as $r$ gets very large, so does $\rho$ and our middle two terms drop out. That leaves us with $$ \begin{align} & \frac{d^2u}{d\rho^2}u -\frac{1}{4}u = 0i,\\ & \rightarrow u = e^{\frac{\pm \rho}{2}},\\ & \frac{d^2u}{d\rho^2}u -\frac{1}{4}e^{\frac{\pm \rho}{2}} =0. \end{align} $$ We in fact discard the positive solution, as it diverges for large $\rho$.

Now in the limit as $r \rightarrow 0$, we see that the terms of import are only the first and second — the other two don't grow quickly enough in the limit. We write $$ \frac{d^2u}{d\rho^2} -\frac{l(l+1)}{\rho^2}u = 0. $$ This is a second order ODE with variable coefficients, and to solve it we need Frobenius's method of assuming a power series solution i.e. that $u(\rho)$ is representable as a power series. We can find part of our $$ \begin{align} & u(\rho) = \sum_{i=0}^{\infty} \alpha_i \rho^{i+b},\\ & u'(\rho) = \sum_{i=0}^{\infty}(i+b) \alpha_i \rho^{i+b},\\ & u''(\rho) = \sum_{i=0}^{\infty}(i+b)(i+b-1) \alpha_i \rho^{i+b-2}.\\ \end{align} $$ Then, $$ \begin{align} & \sum_{i=0}^{\infty}(i+b)(i+b-1) \alpha_i \rho^{i+b-2} - (l(l+1)\rho^{-2} \sum_{i=0}^{\infty} \alpha_i \rho^{i+b}=0\\ & \rightarrow \sum_i \alpha_i \rho^{i+b-2} ((i+b)(i+b-1) - l(l+1)) = 0. \end{align} $$ As our sum $\sum_i \alpha_i \rho^{i+b-2} \neq 0$, we can extract our $i=0$ term and find the allowable solutions for $b$. We get that $$ b(b-1) - l(l+1) \rightarrow b = l+1, -l. $$ Our second asymptotic solution for $u$ then is $\lim_{\rho \rightarrow 0}u(\rho) = A\rho^{l+1}+B\rho^{-l}$, and we see that $B=0$ as otherwise the term blows up to infinity.

At this point, we write a cursory soluton: that $u(\rho) = \rho^{l+1}e^{\frac{-\rho}{2}}L(\rho)$. Plugging this into our original equation and performing a large amount of tedious but ultimately straightforward algebra, we get that $$ \rho \frac{d^2L}{d\rho^2}+ (2l+2-\rho)\frac{dL}{d\rho} + (\lambda -l -1)L = 0. \label{eq:laguerre} \tag{*} $$ To get a recurrence relation for the coefficients of the polynomial $L$, we suppose that $L$ is expressible as a power series in terms of $\rho$ similar to how we found in the Frobenius method above: we say $L= \sum_i a_i\rho^i$, and so $$ \begin{align} & \sum_i(i)(i-1)a_i\rho{i-1}+(2l+2)\sum_ia_i\rho^{i-1} -\sum_ia_i\rho^i +(\lambda -l -1)\sum_ia_i\rho^i = 0\\ & = \sum_i[(i)(i-1) +(2l+2)(i)]a_i\rho^{i-1} - \sum_i[(i)+\lambda-l-1](a_i)\rho^i =0. \end{align} $$ If we adjust the index of the first term to be $i' = i-1$, we can rewrite our formula as $$ \sum_{i}[(i'+1)(i') +(2l+2)(i'+1)]a_{i'+1}\rho^{i'} - \sum_i[(i)-\lambda+l+1](a_i)\rho^i =0. $$ Now everything is of the same power, and we can replace $i'$ with $i$ and establish our recursion relationship: $$ \begin{align} & \sum_{i}[(i+1)(i) +(2l+2)(i+1)]a_{i+1}\rho^{i} - \sum_i[(i)-\lambda+l+1](a_i)\rho^i =0\\ & \rightarrow a_{i+1} = \frac{a_i(i-\lambda+l+1)}{(i+1)(i+2l+2)}. \end{align} $$

We need to make sure these polynomials don't blow up in the limit of large $i$. Taking the ratio between $a_{i+1}$ and $a_i$, we get that this approacjes $\frac{1}{i}$ for large $i$, the same power series ratio between terms of a positive exponential — we need the coefficients to tend to zero if this is going to converge. The recursion relationship is a strictly mathematical result, to have a physical meaning we need to allow for a maximal value of $i$. To find this maximal value, we set $a_{i+1} =0$, and find $0 = i_{\max} - \lambda +l+1$, or $$ \lambda = l+1+i_{\max}. $$ From this relationship, we see that $\lambda$ takes an integer value and that $\lambda \geq l+1$. If $l = 0$, as it is allowed, we have that $\lambda =1$. This is the lowest possible energy-value, showing that there is no zero energy level for the bound energy states of the hydrogen atom. That's pretty cool. We make the notational adjustment $\lambda \rightarrow n$ and call this integer the "principal quantum number". It is possible to show by direct substitution that the energy corresponding to different energy levels is given by the formula $$ E_n = - \frac{\mu e^4}{32 \pi^2 \epsilon_0^2 \hbar^2 n^2}. $$

5.1 The Associated Laguerre Polynomial

Final stretch. A polynomial that satisfies the equation $\eqref{eq:laguerre}$ is called an associtated Laguerre polynomial. These are well-studied, and have the form, again due to Rodrigues, $$ L_n^k(x)= \frac{e^xx^{-k}}{n!} \frac{d^n}{dx^n}(e^{-x}x^{n+k}). $$ They have the orthogonality and recursion relationships $$ \begin{align} & \int_0^{\infty} e^{-x}x^kL_n^k(x)L_m^k(x) = \frac{(n+k)!}{n!} \delta_{nm},\\ & xL_n^k = (2n+k+1)L_m^k - (n+k)L_{n-1}^k - (n+1)L_{n+1}^ka. \end{align} $$ We need to normalize this radial solution, and these relationships will help us do that. We remember that the form of our wavefunction is $\psi_{nlm}(r,\theta, \phi) = N_{n,l,m} R_{n,l}(r)*Y_{l,m}$, where N is our normalization constant.

At this point, we write again our cursory solution — $u(\rho) = e^{\frac{-\rho}{2}}\rho^{l+1}L_{n-l-1}^{2n+1}$. Realizing that $u(\rho) = \rho R(\rho)$, we can say that $$ R(\rho) = e^{\frac{-\rho}{2}}\rho^{l}L_{n-l-1}^{2l+1}. $$ Here's the grand finale. Our entire wavefunction must follow the the relationship $$ \int_0^{2\pi}\int_0^{\pi}\int_0^{\infty} |\psi_{nlm}|^2 r^2\sin^2\theta dr d\theta d\phi =1. $$ Having already normalized the angular portion, we make the substitution $\rho = 2\kappa r$, $d\rho = 2\kappa dr$, and proceed with the following: $$ \begin{align} & \frac{N^2}{(2\kappa)^3}\int_0^{\infty} e^{-\rho}\rho^{2l+2} \left(L_{n-l-1}^{2l+1}\right)^2\\ & = \frac{N^2}{(2\kappa)^3}\int_0^{\infty} e^{-\rho}\rho^{2l+1} \left(\rho L_{n-l-1}^{2l+1}\right) \left(L_{n-l-1}^{2l+1}\right)\\ & = \frac{N^2}{(2\kappa)^3}\frac{2n(n+l)!}{(n-l-1)!} = 1\\ &\rightarrow N= \sqrt{\frac{(2\kappa)^3(n-l-1)!}{2n(n+l)!}}. \end{align} $$ We may now write the wavefunction, explicitly, and in full: $$ \psi_{n,l,m}(r,\theta,\phi) = \alpha \left(\frac{2}{na_0}\right)^{\frac{3}{2}} \sqrt{\frac{(n-l-1)!}{(2n)(n+l)!}}e^{\frac{-r}{na_0}} \left(\frac{2r}{na_0}\right)^lL_{n-l-1}^{2l+1}\left( \frac{2r}{na_0}\right)\sqrt{\frac{(2l+1)(l-|m|)!} {4\pi(l+|m|)!}}e^{im\phi}P_l^m(\cos\theta). $$ This, in fact, may be generalized slightly for "hydrogen-like atoms", those with a single electron in the valence. We send $\kappa \rightarrow Z\kappa$, $Z$ being the atomic number — the number of protons in the nucleus.

6. Orbital Angular Momentum

From here, it's a procedural matter to find the orbital angular momentum: all one needs to do is apply the angular momentum operator, $\vec{\textbf{L}}^2$, to the wave function for hydrogen in order to determine its azimuthal quantum number and determine the orbital angular momentum, $$ \vec{\textbf{L}}^2\ket{\psi_{nlm}} = \hbar l(l+1)\ket{\psi_{nlm}} $$ or you might apply the operator $\vec{\textbf{L}}_z$ to find the magnetic quantum number: $$ \vec{\textbf{L}}_z\ket{\psi_{nlm}} = m_l\hbar \ket{\psi_{nlm}}. $$

But this is not the total angular momentum of an electron orbiting a proton, and it is at this point we need to revisit spin. The initial work to describe the total angular momentum was done by Pauli in the late 1920's, but the equation he came up with is not derivable from the strict limitations of the Schrodinger equation — instead, we must ask a relativistic question, come up with a relativistic answer, and take a non-relativistic limit to come to the answer we know today. We will need Dirac's physics, but to get to it we first look at the too-often overlooked Klein-Gordon equation.

7. The Klein-Gordon Equation

We will be working in natural units for this section: $\hbar = c = 1$. Taking a second time derivative of the Schrodinger equation, multiplying again by $i$, we recover $$ -\frac{\partial^2}{\partial t^2}\ket{\psi(t)} = H^2 \ket{\psi(t)}. $$ Using our knowledge, from Einstein, that $H^2 = p^2 + m^2$, defining $\Psi(\textbf{x},t) \equiv \braket{\textbf{x}|\psi(t)}$ and the fact that $\bra{\textbf{x}}p^2\ket{\psi(t)} = - \nabla^2 \psi(\textbf{x},t)$, we can express this result as $$ \left[ \frac{\partial^2}{\partial t^2} - \nabla^2 +m \right] \Psi(\textbf{x},t) = 0. $$ This is the Klein-Gordon equation. This is a notable equation for many reasons, the first of which is that it it provides a covariant wavefunction. This is more easily realised by observing that Lorentz transformation, which leaves the square of the spacetime interval invariant implies that the derivatives in this equation operate the same way on $(t,\textbf{x})$ as $(t',\textbf{x}')$.

Adopting relativistic notation, the Klein-Gordgon equation can be rewritten as $$ \left[\partial_{\mu}\partial^{\mu} +m^2\right] \Psi(\textbf{x},t) = 0. $$ The Klein-Gordon equation also fulfills our expectations for the solutions of a free relativistic particle of mass $m$. It has time dependence like $e^{-iEt}$, with $E$ an eigenvalue of the Hamiltonian, and spacial dependence like a plane wave: i.e. $e^{i\textbf{p}\cdot\textbf{x}}$. This is to say, our solution is of the form $$ \Psi(\textbf{x},t) = Ne^{-ip^{\mu}x_{\mu}},\\ p^{\mu} = (E,\textbf{p}). $$ This indeed does solve the equation as long as $-p^{\mu} p_{\mu}+m^2 = -E^2+\textbf{p}^2+m^2=0.$ That's a plus to be sure. However, the Klein-Gordon equation is not without its problems. In further investigation, one realises that this equation does not yeild a postive definite probability density. It is not a problem without solution and a physical interpretation — but unfortunately, the solutions were ill-timed. The difficulties all stem from the fact that the equation is second-order in space and time. It took Dirac's stubborn assertion that a solution existed linear in time and space for relativistic quantum mechanics to really take off: and a byproduct of this is, finally, an explanation for spin.

8. The Dirac Equation

What Dirac asserted was an equation like the following: $$ (i\gamma^{\mu}\partial_{\mu}-m)\Psi(\textbf{x},t) = 0. $$ We must determine this factor $\gamma^{\mu}$, and we insist that our results manage to preserve the the correct energy eigenvalues for a free particle, those we found in our solution for the Klein-Gordon equation. But we can turn this equation into a form resembling the Klein-Gordon equation by operating on it with its conjugate! In so doing, we come up with $$ (\gamma^{\nu}\partial_{\nu}\gamma^{\mu}\partial_{\mu}+m^2) \Psi(\textbf{x},t)=0,\\ \gamma^{\nu}\partial_{\nu}\gamma^{\mu}\partial_{\mu} = \gamma^{\nu}\gamma^{\mu}\partial_{\nu}\partial_{\mu} = \partial^{\mu}\partial_{\mu}=\eta^{\mu \nu}\partial_{\nu} \partial_{\mu}. $$ We can express this condition by symmetrizing: $$ \frac{1}{2} (\gamma^{\mu}\gamma^{\nu}+\gamma^{\nu}\gamma^{\mu})\equiv \frac{1}{2}\{\gamma^{\mu},\gamma^{\nu}\} = \eta^{\mu \nu}. $$ These entities $\gamma^{\mu}$, then, are not simply the complex numbers we are familiar with but rather elements of a Clifford algebra. Given our values for $\eta^{\mu \nu}$, we are can infer that $(\gamma^0)^2 = 1$, $(\gamma^i)^2 = -1$, and $\gamma^{\mu}\gamma^{\nu} = -\gamma^{\nu}\gamma^{\mu}$ if $\mu \neq \nu$.

If we substitute the free particle solution from the Klein-Gordon equation into the Dirac equation, we find $$ \gamma^{\mu}p_{\mu}-m = 0. $$ If we write this in its timelike and spacelike components and multiply by $\gamma^0$, we can rearrange this to be $$ E = \gamma^0\boldsymbol{\gamma}\cdot\textbf{p}+ \gamma^0m. $$ This leads to what has become known as the Dirac Hamiltonian: $H = \boldsymbol{\alpha}\cdot\textbf{p}+ \beta m$. The Clifford algebra can be fulfilled by $\boldsymbol{\alpha}$ and $\beta$ provided they are at least 4-dimensional matrices. We also insist they are Hermitian. We choose to construct these with the Pauli spin matrices, $$ \boldsymbol{\alpha}= \begin{pmatrix} \boldsymbol{\sigma}&0\\ 0&\boldsymbol{\sigma} \end{pmatrix}, \beta = \begin{pmatrix} \mathbb{1}&0\\ 0&-\mathbb{1} \end{pmatrix}. $$ This makes $\Psi(\textbf{x},t)$ a 4-dimensional column matrix. Examining the Dirac equation in the case of a stationary free particle, we can write $$ i\partial_t \Psi = \beta m \Psi. $$ We arrive at 4 independent solutions for the wavefunction, the upper half corresponding to positive energies and the lower half to negative. It will turn out that the two allowable values for each energy state will correspond to spin-up and spin-down states. Our solutions are given by: $$ \Psi_1 = e^{-imt} \begin{bmatrix} 1\\ 0\\ 0\\ 0 \end{bmatrix},\ \Psi_2 = e^{-imt} \begin{bmatrix} 0\\ 1\\ 0\\ 0 \end{bmatrix},\ \Psi_3 = e^{imt} \begin{bmatrix} 0\\ 0\\ 1\\ 0 \end{bmatrix},\ \Psi_4 = e^{imt} \begin{bmatrix} 0\\ 0\\ 0\\ 1 \end{bmatrix}. $$ Pressing forward, for a momentum in the $z$-direction, our Hamiltionian become $$ \begin{bmatrix} m & 0 & p & 0\\ 0 & m & 0 & -p\\ p & 0 & -m & 0\\ 0 & -p & 0 & -m \end{bmatrix} \begin{bmatrix} u_1\\ u_2\\ u_3\\ u_4 \end{bmatrix} = E \begin{bmatrix} u_1\\ u_2\\ u_3\\ u_4 \end{bmatrix}. $$ We notice the coupling of the 1-3 states and the 2-4 states, which turn out to be analagous spins between particle and antiparticle pairs. Going through the motions, it is possible to determine energy eigenvalues of $E = \pm \sqrt{m^2+p^2}$, the positive energies corresponding to matter and the negative energies to antimatter. For the positive energy case, where $E_p = +\sqrt{p^2+m^2}$: $$ u^{+}_{\uparrow}(p) = \begin{bmatrix} 1\\ 0\\ \frac{p}{E_p+m}\\ 0 \end{bmatrix},\ u^{+}_{\downarrow}(p) = \begin{bmatrix} 0\\ 1\\ 0\\ \frac{-p}{E_p+m} \end{bmatrix}. $$ The upper components in both of these vectors are dominating in the nonrelativistic case. If we consider the behaviour of the operator $\boldsymbol{\Sigma} \cdot \textbf{p} = \Sigma_z$, where $$ \boldsymbol{\Sigma} = \begin{pmatrix} \boldsymbol{\sigma} & 0\\ 0 & \boldsymbol{\sigma} \end{pmatrix}. $$ We can see that $\textbf{S}=\frac{\hbar}{2}\boldsymbol {\Sigma}$ does indeed project the expected spin-up result where $u_1 \neq 0$ and the expected spin-down result for $u_2 \neq 0$. There are analagous results for the spin states in the negative energy regime.

That's been a serious digression, but an important one if we are to understand the relationship between angular momentum, spin, and the magnetic moment of fermions. We may introduce electromagnetic interactions via gauge transformation: $\tilde{\textbf{p}} \equiv \textbf{p}-e \textbf{A}$, and the dirac equation may be expressed in the form of a $2 \times 2$ matix: $$ \begin{bmatrix} m & \boldsymbol{\sigma} \cdot \tilde{\textbf{p}}\\ \boldsymbol{\sigma} \cdot \tilde{\textbf{p}} & -m\\ \end{bmatrix} \begin{bmatrix} u\\ v \end{bmatrix} = E \begin{bmatrix} u\\ v \end{bmatrix}. $$ For nonrelativistic energy $E = T + m, T \ll m$ and the bottom equation is given by $$ \boldsymbol{\sigma} \cdot \tilde{\textbf{p}}u =(E+m)v \approx 2mv. $$ This lets us write the upper equation as $$ \frac{(\boldsymbol{\sigma} \cdot \tilde{\textbf{p}}) (\boldsymbol[{\sigma} \cdot \tilde{\textbf{p}})}{2m}u =\left[\frac{\tilde{\textbf{p}}^2}{2m} +\frac{i\boldsymbol{\sigma}}{2m}\cdot (\tilde{\textbf{p}} \times \tilde{\textbf{p}})\right]u =Tu. $$ The cross product $\tilde{\textbf{p}} \times \tilde{\textbf{p}}u$ in coordinate representation can be written as $$ \begin{gather} \tilde{\textbf{p}} \times \tilde{\textbf{p}}u =(i \nabla +e\textbf{A}) \times (i \nabla u +e\textbf{A}u)\\ ie[\nabla \times (\textbf{A}u)+ \textbf{A} \times \nabla u] \\ =ie(\nabla \times \textbf{A})u = ie\textbf{B}u. \end{gather} $$ So our upper equation is $$ \left[\frac{\tilde{\textbf{p}}^2}{2m} - \boldsymbol{\mu} \cdot \textbf{B}\right]u = Tu, $$ where $\boldsymbol{\mu} = g\frac{e}{2m}\textbf{S}$, and $g \approx 2$ for an electron's intrinsic spin.

We have seen that to answer the question about the relation between spin and magnetic moment, we have needed to develop much of quantum mechanical theory, and finally invoke relativistic physics. It is a spectacularly deep question!

9. Total Angular Momentum and Quantum Numbers

We will close with some final remarks about the measurement of total angular momentum, relating to our exploration of the hydrogen atom. We choose a frame with zero vector potential, for simplicity. We consider the commutation relation between the Dirac Hamiltonian and the angular momentum operator. The $\beta$ portion clearly commutes, so we are left with $$ \begin{gather} [\boldsymbol{\alpha}\cdot\textbf{p},L_i] = [\alpha_l p_l, \epsilon_{ijk}x_jp_k]\\ = \epsilon_{ijk}\alpha_l[p_l,x_j]p_k \\ = -i \epsilon_{ijk}\alpha_j p_k. \end{gather} $$ The Dirac Hamiltonian does not commute with the orbital angular momentum, and thus orbital angular momentum is not a conserved quantity for spin-$\frac{1}{2}$ particles! But what about the operator, $\Sigma$? $$ [\boldsymbol{\alpha}\cdot\textbf{p},\Sigma_j] = [\alpha_i,\Sigma_j]p_i = 2i \epsilon_{ijk}\alpha_k p_i. $$ So spin, by itself, is not conserved either! However, if we modify the operator $\boldsymbol{\Sigma}$ by a factor of $\frac{\hbar}{2}$ and add it to the orbital angular momentum operator $\textbf{L}$, we can see that $$ \begin{gather} \textbf{J} = \textbf{L}+ \frac{\hbar}{2}\boldsymbol{\Sigma} = \textbf{L} + \textbf{S},\\ [H,\textbf{J}] = 0. \end{gather} $$ Thus the total angular momentum is a conserved quantity. That the total is conserved but the components are not implies that the components can and do trade off with each other: this is a phenomenon known as the "Spin-Orbit Interaction". To go into further elaboration on this phenomenon would digress too far even for me. However, the spin-orbit interaction does allow us to describe total angular momentum with fewer quantum numbers used below.

At last, the entire quantum angular momentum of an electron may, now, be characterized by the quantum numbers $l$ and $s$, the "azimuthal" and "spin" quantum numbers, respectively. For the total angular momentum $\textbf{J}$, we get a total angular momentum number $j = l+s$ which is bounded by $$ |l-s| \leq j \leq l+s. $$ The azimuthal quantum number is determined by the orbital in which it finds itself, and is itself bound by the principal quantum number $n$ which determines the energy level: $$ n = 1,2,3, \ldots, $$ The azimuthal quantum number may range from $$ 0,1,\ldots, n-1 $$ and describe the shape of the atomic orbitals. We also have the magnetic quantum number $m$, which dictate the individual orbital structure from the class in which it belongs. Accordingly, its range is $$ m = -l, -l+1, \ldots, l-1, l. $$ And lastly, there is the spin quantum number, which can take only one of two values: $$ s = \pm \frac{1}{2}. $$