A closer look at how the principal, orbital and magnetic quantum numbers arise in Schrödinger’s theory of the hydrogen atom

SchrodingerIn this note I want to explore some aspects of the solution of Schrödinger’s three-dimensional wave equation in spherical polar coordinates which pertain to the three main quantum numbers characterising the electron in a hydrogen atom: the electron’s principal, orbital and magnetic quantum numbers. In particular, I have noticed that many physics and mathematics sources tend to gloss over some (or all) of the mathematical details of how permissibility conditions on the three quantum numbers emerge naturally from validity constraints on the relevant underlying differential equations. I want to bring out these details clearly in the present note.

In general, four quantum numbers are needed to fully describe atomic electrons in many-electron atoms. These four numbers and their permissible values are:

Principal quantum number \tilde{n} = 1, 2, 3, \ldots

Orbital quantum number l = 0, 1, 2, \ldots, (\tilde{n} - 1)

Magnetic quantum number m_l = 0, \pm 1, \pm 2, \ldots, \pm l

Spin magnetic quantum number m_s = -\frac{1}{2}, +\frac{1}{2}

The principal quantum number determines the electron’s energy, the orbital quantum number its orbital angular-momentum magnitude, the magnetic quantum number its orbital angular-momentum direction, and the spin magnetic quantum number its spin direction.

I have noticed that it is often not explained clearly why, for example, the orbital quantum number cannot exceed the principal quantum number minus one, or why the magnitude of the magnetic quantum number cannot exceed that of the orbital quantum number. I will explore issues like these in depth in the context of the time-independent Schrödinger equation for the hydrogen atom, which only involves the first three quantum numbers. I will not discuss the spin magnetic quantum number in the present note.

Schrödinger’s wave equation for the hydrogen atom

In Cartesian coordinates, Schrödinger’s three-dimensional equation for the electron in the hydrogen atom is

\frac{\partial^2 \psi}{\partial x^2} + \frac{\partial^2 \psi}{\partial y^2} + \frac{\partial^2 \psi}{\partial z^2} + \frac{2m_e}{\hbar^2}(E - U) \psi = 0

where m_e denotes the electron mass. The potential energy U is the electric potential energy of a charge -e given that it is at distance r from another charge +e, namely

U = -\frac{e^2}{4 \pi \epsilon_0 r}

It is necessary to change variables in Schrödinger’s equation since the potential energy is a function of radial distance r rather than the Cartesian coordinate variables x, y and z. Given the spherical symmetry of the atom, it is sensible to proceed by changing the variables in Schrödinger’s equation to those of spherical polar coordinates (rather than changing the r variable in U to Cartesian coordinates using r = \sqrt{x^2 + y^2 + z^2}). Only the variables in the Laplacian part of Schrödinger’s equation need to be changed, so we can use the approach in my previous note on changing variables in Laplace’s equation to get

\frac{1}{r^2} \frac{\partial }{\partial r}\big( r^2 \frac{\partial \psi}{\partial r}\big) + \frac{1}{r^2 \sin \theta}\frac{\partial }{\partial \theta}\big( \sin \theta \frac{\partial \psi}{\partial \theta} \big) + \frac{1}{r^2 \sin^2 \theta}\frac{\partial^2 \psi}{\partial \phi^2} + \frac{2m_e}{\hbar^2}(E - U) \psi = 0

I will now temporarily simplify things by using the representation of the square of the angular momentum operator in spherical polar coordinates which I obtained in another previous note, namely

L^2 = - \hbar^2 \big( \frac{1}{\sin \theta}\frac{\partial}{\partial \theta} \big( \sin \theta \frac{\partial}{\partial \theta}\big) + \frac{1}{\sin^2 \theta}\frac{\partial^2}{\partial \phi^2} \big)

= - \hbar^2 r^2 \big( \frac{1}{r^2 \sin \theta}\frac{\partial}{\partial \theta} \big( \sin \theta \frac{\partial}{\partial \theta}\big) + \frac{1}{r^2 \sin^2 \theta}\frac{\partial^2}{\partial \phi^2} \big)

Using this to replace the two middle terms in Schrödinger’s equation and rearranging we get

\frac{1}{r^2} \frac{\partial }{\partial r}\big( r^2 \frac{\partial \psi}{\partial r}\big) + \frac{2m_e}{\hbar^2}(E - U) \psi = \frac{L^2}{\hbar^2 r^2} \psi

This equation can now be solved by the usual separation of variables approach. We assume that the \psi function can be expressed as a product

\psi(r, \theta, \phi) = R(r) Y(\theta, \phi)

and then substitute this back into the wave equation to get

\frac{Y}{r^2} \frac{d}{ d r}\big( r^2 \frac{d R}{d r}\big) + \frac{2m_e}{\hbar^2}(E - U) R Y = \frac{R}{\hbar^2 r^2} L^2 Y

Dividing through by \frac{R Y}{r^2} we get

\frac{1}{R} \frac{d}{d r}\big( r^2 \frac{d R}{d r}\big) + \frac{2m_e r^2}{\hbar^2}(E - U) = \frac{1}{\hbar^2 Y} L^2 Y

Since the left-hand side of this equation depends only on r while the right-hand side depends only on \theta and \phi, both sides must be equal to some constant which we can call \lambda. Setting the left and right-hand sides equal to \lambda in turn and rearranging slightly we finally get the radial equation

\frac{1}{r^2} \frac{d}{d r}\big( r^2 \frac{d R}{d r}\big) + \big[ \frac{2m_e}{\hbar^2}(E - U) - \frac{\lambda}{r^2} \big] R = 0

and the angular equation

L^2 Y = \lambda \hbar^2 Y

We can now apply separation of variables again to the angular equation. Rewriting the operator L^2 in full the angular equation becomes

- \hbar^2 \big( \frac{1}{\sin \theta}\frac{\partial}{\partial \theta} \big( \sin \theta \frac{\partial Y}{\partial \theta}\big) + \frac{1}{\sin^2 \theta}\frac{\partial^2 Y}{\partial \phi^2} \big) = \lambda \hbar^2 Y

which simplifies to

\frac{1}{\sin \theta}\frac{\partial}{\partial \theta} \big( \sin \theta \frac{\partial Y}{\partial \theta}\big) + \frac{1}{\sin^2 \theta}\frac{\partial^2 Y}{\partial \phi^2} + \lambda Y = 0

We assume that the Y function can be written as the product

Y(\theta, \phi) = \Theta (\theta) \Phi (\phi)

Substituting this into the angular equation gives

\frac{\Phi}{\sin \theta}\frac{d}{d \theta} \big( \sin \theta \frac{d \Theta}{d \theta}\big) + \frac{\Theta}{\sin^2 \theta}\frac{d^2 \Phi}{d \phi^2} + \lambda Y \Theta \Phi = 0

Multiplying through by \frac{\sin^2 \theta}{\Theta \Phi} and rearranging we get

-\frac{1}{\Phi} \frac{d^2 \Phi}{d \phi^2} = \frac{\sin^2 \theta}{\Theta}\big[ \frac{1}{\sin \theta}\frac{d}{d \theta}\big(\sin \theta \frac{d \Theta}{d \theta}\big) + \lambda \Theta\big]

Since the left-hand side of this equation depends only on \phi while the right-hand side depends only on \theta, both sides must be equal to some constant which we can provisionally call k. Setting the left and right-hand sides equal to k in turn and rearranging we get

\frac{d^2 \Phi}{d \phi^2} + k \Phi = 0

and

\frac{1}{\sin \theta}\frac{d}{d \theta}\big(\sin \theta \frac{d \Theta}{d \theta}\big) + \big( \lambda - \frac{k}{\sin^2 \theta}\big) \Theta = 0

We now have three ordinary differential equations that need to be solved, one for \Phi, one for \Theta and one for R. We will solve each of them in turn.

The equation for \Phi

The equation for \Phi is a straightforward second-order differential equation with auxiliary equation

\zeta^2 + k = 0

implying \zeta = \pm \sqrt{-k} if k < 0 and \zeta  = \pm i \sqrt{k} if k > 0. Therefore it has a general solution of the form

\Phi(\phi) = Ae^{\sqrt{-k} \phi} + Be^{- \sqrt{-k} \phi}

if k < 0 and \Phi(\phi) = Ae^{i \sqrt{k} \phi} + Be^{- i \sqrt{k} \phi} if k > 0, where A and B are arbitrary constants. Now, the azimuth angle \phi can take any value in (-\infty, \infty) but the function \Phi must take a single value at each point in space (since this is a required property of the quantum wave function which \Phi is a constituent of). It follows that the function \Phi must be periodic since it must take the same value at \phi and \phi + 2\pi for any given \phi. This imposes two constraints on the form of the general solution: (1) it cannot consist only of exponential functions with real arguments since these are not periodic (thus ruling out the first general solution above and thereby implying that the separation constant k must be nonnegative); (2) \sqrt{k} must be an integer. Given these constraints, it is customary in quantum mechanics to denote \pm \sqrt{k} by the letter m (it is called the magnetic quantum number) and to specify the separation constant in the angular equations as m^2, which guarantees its nonnegativity. We then state the general solution of the equation for \Phi as

\Phi(\phi) = Ae^{i m \phi} + Be^{- i m \phi}

In practice, a particular electron wave function will involve only a single value of m so only the first of the two terms in the general solution will be necessary. We can therefore state the general solution of the equation for \Phi for a given magnetic quantum number m as

\Phi(\phi) \propto e^{i m \phi}

The equation for \Theta

Given that we now know the separation constant for the angular equations is either zero or a positive square number k = m^2, we can write the equation for \Theta as

\frac{1}{\sin \theta}\frac{d}{d \theta}\big(\sin \theta \frac{d \Theta}{d \theta}\big) + \big( \lambda - \frac{m^2}{\sin^2 \theta}\big) \Theta = 0

Expanding the first term we get

\frac{1}{\sin \theta} \cos \theta \frac{d \Theta}{d \theta} + \frac{d^2 \Theta}{d \theta^2} + \big( \lambda - \frac{m^2}{\sin^2 \theta}\big) \Theta = 0

I am now going to multiply and divide the first two terms by \sin^2 \theta to get

\sin^2 \theta \big(\frac{\cos \theta}{\sin^3 \theta} \frac{d \Theta}{d \theta} + \frac{1}{\sin^2 \theta} \frac{d^2 \Theta}{d \theta^2} \big) + \big( \lambda - \frac{m^2}{\sin^2 \theta}\big) \Theta = 0

\iff

\sin^2 \theta \big(- \frac{\cos \theta}{\sin^3 \theta} \frac{d \Theta}{d \theta} + \frac{1}{\sin^2 \theta} \frac{d^2 \Theta}{d \theta^2} + \frac{2 \cos \theta}{\sin^3 \theta} \frac{d \Theta}{d \theta}\big) + \big( \lambda - \frac{m^2}{\sin^2 \theta}\big) \Theta = 0

Now we can make the change of variable x = \cos \theta which implies dx = - \sin \theta d \theta and therefore

\frac{d \theta}{dx} = - \frac{1}{\sin \theta}

\frac{d \Theta}{d x} = \frac{d \Theta}{d \theta} \frac{d \theta}{d x} = - \frac{1}{\sin \theta} \frac{d \Theta}{d \theta}

\frac{d^2 \Theta}{d x^2} = \frac{d}{d \theta} \big[ - \frac{1}{\sin \theta} \frac{d \Theta}{d \theta} \big] \frac{d \theta}{d x} = - \frac{\cos \theta}{\sin^3 \theta} \frac{d \Theta}{d \theta} + \frac{d^2 \Theta}{d \theta^2}

Using these in the amended form of the \Theta equation together with the fact that \sin^2 \theta = 1 - x^2, the \Theta equation becomes

(1 - x^2) \big(\frac{d^2 \Theta}{d x^2} - \frac{2x}{1 - x^2} \frac{d \Theta}{d x} \big) + \big(\lambda - \frac{m^2}{1 - x^2} \big) \Theta = 0

\iff

(1 - x^2) \frac{d^2 \Theta}{d x^2} - 2x \frac{d \Theta}{d x} + \big(\lambda - \frac{m^2}{1 - x^2} \big) \Theta = 0

We will solve this equation first for the case m = 0 (the solutions will be Legendre polynomials) and use these results to construct solutions for the case m \neq 0 (the solutions here will be the associated Legendre functions). Setting m = 0 we get

(1 - x^2) \frac{d^2 \Theta}{d x^2} - 2x \frac{d \Theta}{d x} + \lambda \Theta = 0

which has the form of a well known differential equation known as Legendre’s equation. It can be solved by assuming a series solution of the form

\Theta = a_0 + a_1 x + a_2 x^2 + a_3 x^3 + a_4 x^4 + \cdots + a_n x^n + \cdots

and then differentiating it term by term twice to get

\Theta^{\prime} = a_1  + 2 a_2 x + 3 a_3 x^2 + 4 a_4 x^3 + \cdots + n a_n x^{n-1} + \cdots

and

\Theta^{\prime \prime} = 2 a_2 + 6 a_3 x + 12 a_4 x^2 + 20 a_5 x^3 + \cdots + n (n - 1) a_n x^{n-2} + \cdots

We now substitute these into Legendre’s equation and set the coefficient of each power of x equal to zero (because \Theta must satisfy Legendre’s equation identically). We find that the coefficient of the x^n term satisfies

(n + 2)(n + 1) a_{n + 2} + (\lambda - n(n + 1)) a_n = 0

which implies

a_{n + 2} = - \frac{(\lambda - n(n + 1))}{(n + 2)(n + 1)} a_n

This formula makes it possible to find any even coefficient as a multiple of a_0 and any odd coefficient as a multiple of a_1. The general solution of our Legendre equation is then a sum of two series involving two arbitrary constants a_0 and a_1:

\Theta = a_0 \bigg \{1 - \frac{\lambda}{2!} x^2 + \frac{\lambda (\lambda - 6)}{4!} x^4 - \frac{\lambda(\lambda - 6)(\lambda - 20)}{6!} x^6 + \cdots \bigg \}

+ a_1 \bigg \{x - \frac{(\lambda - 2)}{3!} x^3 + \frac{(\lambda - 2)(\lambda - 12)}{5!} x^5 - \frac{(\lambda - 2)(\lambda - 12)(\lambda - 30)}{7!} x^7 + \cdots \bigg \}

Both of the series in this sum converge for x^2 < 1 but in general they do not converge for x^2 = 1. This is a problem for us because in our change of variables we set x = \cos \theta and we want solutions that converge for all possible values of \theta including those that result in x^2 = 1. It turns out that the only way to get such solutions is to choose integer values of \lambda that make either the a_0 or the a_1 series in the above sum terminate (the other series will generally be divergent so we remove it by setting the corresponding arbitrary constant equal to zero). This requires \lambda to take values in the quadratic sequence 0, 2, 6, 12, 20, 30, 42, 56 \ldots The l-th term of this sequence is l(l + 1), so the separation constant \lambda must be of this form, i.e., \lambda = l(l + 1) for some l = 0, 1, 2, 3, \dots. When l takes an even value the a_0 series will terminate and we can set a_1 = 0 to make the other series vanish. Conversely, when l takes an odd value the a_1 series will terminate and we can set a_0 = 0 to make the other series vanish.

From the eigenvalue equation for L^2 given earlier (L^2 Y = \lambda \hbar^2 Y) it is clear that the magnitude of the orbital angular momentum is L = \sqrt{l(l + 1)} \hbar. It is interesting to see how the form of this arises mathematically from considering series solutions to Legendre’s equation above. The parameter l is called the orbital angular momentum quantum number.

Note that negative integral values of l are allowed but they simply give solutions already obtained for positive values. For example, l = -2 gives \lambda = 2 and this makes the a_1 series terminate, yielding the polynomial solution

\Theta = a_1 x

This is exactly the same solution as the one that would be obtained if l = 1. It is therefore customary to restrict l to nonnegative values. Each possible value of l gives a polynomial solution to Legendre’s equation. For l = 0 we get \Theta = a_0, for l = 1 we get \Theta = a_1 x, for l = 2 we get \Theta = a_0 - 3 a_0 x^2, and so on. If the value of a_0 or a_1 in each polynomial equation is selected so that \Theta = 1 when x = 1 the resulting polynomials are called Legendre polynomials, denoted by P_l(x). Given that for each l we have P_l(1) = 1 the first few Legendre polynomials are

P_0(x) = 1 P_1(x) = x

P_2(x) = \frac{1}{2}(3 x^2 - 1)

P_3(x) = \frac{1}{2}\big(5 x^3 - 3x \big)

These are the physically acceptable solutions to Legendre’s equation for \Theta above.

We now consider the solutions for

m \neq 0 of the equation (1 - x^2) \frac{d^2 \Theta}{d x^2} - 2x \frac{d \Theta}{d x} + \big(\lambda - \frac{m^2}{1 - x^2} \big) \Theta = 0

We now know that \lambda = l(l + 1) so we can write this in and we can also add the subscript l to m as the solutions to this equation will involve a link the between the values of the orbital angular momentum and magnetic quantum numbers. The equation we need to solve becomes

(1 - x^2) \frac{d^2 \Theta}{d x^2} - 2x \frac{d \Theta}{d x} + \big[l(l + 1) - \frac{m_l^2}{1 - x^2} \big] \Theta = 0

The link between l and m_l arises from the fact that we are constrained in trying to solve this equation: it encompasses the case m_l = 0 for which the physically acceptable solutions are the Legendre polynomials P_l(x). Therefore the physically allowable solutions for the above equation must include the Legendre polynomials as a special case. We can find these by using the series approach again and it turns out that the physically acceptable solutions are the so-called associated Legendre functions which take the form

P_l^{m_l}(x) = (1 - x^2)^{m_l/2}\frac{d^{m_l}}{d x^{m_l}}P_l(x)

Now, each Legendre polynomial P_l(x) is a polynomial of degree l. Therefore the m_l-th order derivative in P_l^{m_l} will equal zero if |m_l| > l, so for physically acceptable solutions we must impose the constraint |m_l| \leq l in the differential equation for \Theta. This is where the link between the quantum numbers l and m_l comes from in the quantum theory of the hydrogen atom: given a value of l the acceptable values of m_l are integers in the range -l \leq m_l \leq l.

Finally, note two things: (1) P_l^{m_l}(x) reduces to the Legendre polynomial P_l(x) when m_l = 0, which is what we needed. (2) A negative value for m_l does not change m_l^2 in the original differential equation so a solution for positive m_l is also a solution for the corresponding negative m_l. Thus many references define the associated Legendre function P_l^{m_l}(x) for -l \leq m_l \leq l as P_l^{|m_l|}(x).

To conclude, given values for the quantum numbers l and m_l, the general solution of the equation for \Theta can be written as

\Theta(\theta) \propto P_l^{m_l}(\cos \theta)

The radial equation for R

To clarify where the principal quantum number comes from, the final equation we need to deal with is the radial equation

\frac{1}{r^2} \frac{d}{d r}\big( r^2 \frac{d R}{d r}\big) + \big[ \frac{2m_e}{\hbar^2}(E - U) - \frac{\lambda}{r^2} \big] R = 0

Writing \lambda = l(l + 1) and replacing U with the formula for the potential energy we get

\frac{1}{r^2} \frac{d}{d r}\big( r^2 \frac{d R}{d r}\big) + \big[ \frac{2m_e}{\hbar^2}\big(\frac{e^2}{4 \pi \epsilon_0 r} + E\big) - \frac{l(l + 1)}{r^2} \big] R = 0

\iff

\frac{d^2 R}{d r^2} + \frac{2}{r} \frac{d R}{d r} + \frac{2m_e}{\hbar^2} \big[ E + \frac{e^2}{4 \pi \epsilon_0 r} - \frac{l(l + 1) \hbar^2}{2 m_e r^2} \big] R = 0

We are only interested in solutions for which the electron is bound within the atom, so we take E < 0 (the negative energy of the electron is the amount of energy that must be supplied to it to free it from the atom). In order to solve the above equation it is customary to make the change of variable

\rho = \big(-\frac{8 m_e E}{\hbar^2}\big)^{1/2} r

and define the dimensionless constant

\tau = \frac{e^2}{4 \pi \epsilon_0 \hbar} \big(-\frac{m_e}{2 E} \big)^{1/2}

If we then specify R = R(\rho) we have

\frac{d R}{d r} = \frac{d R}{d \rho} \frac{d \rho}{d r} = \big(-\frac{8 m_e E}{\hbar^2}\big)^{1/2} \frac{d R}{d \rho}

\frac{2}{r} \frac{d R}{d r} = \big(-\frac{8 m_e E}{\hbar^2}\big) \frac{2}{\rho} \frac{d R}{d \rho}

\frac{d^2 R}{d r^2} = \big(-\frac{8 m_e E}{\hbar^2}\big)^{1/2} \frac{d^2 R}{d \rho^2} \frac{d \rho}{d r} = \big(-\frac{8 m_e E}{\hbar^2}\big) \frac{d^2 R}{d \rho^2}

\frac{2 m_e}{\hbar^2} E + \frac{2 m_e}{\hbar^2} \frac{e^2}{4 \pi \epsilon_0 r} = \big(-\frac{8 m_e E}{\hbar^2}\big) \big \{\frac{1}{4}\frac{e^2}{4 \pi \epsilon_0 r} \big(-\frac{1}{E}\big) - \frac{1}{4}\big \} = \big(-\frac{8 m_e E}{\hbar^2}\big) \big(\frac{\tau}{\rho} - \frac{1}{4}\big)

\frac{l(l + 1)}{r^2} = \big(-\frac{8 m_e E}{\hbar^2}\big) \frac{l(l + 1)}{\rho^2}

Using these results we can rewrite the differential equation as

\big(-\frac{8 m_e E}{\hbar^2}\big) \bigg \{\frac{d^2 R}{d \rho^2} + \frac{2}{\rho} \frac{d R}{d \rho} + \big[\frac{\tau}{\rho} - \frac{1}{4} - \frac{l(l + 1)}{\rho^2}\big] R(\rho) \bigg \} = 0

\iff

\frac{d^2 R}{d \rho^2} + \frac{2}{\rho} \frac{d R}{d \rho} + \big[\frac{\tau}{\rho} - \frac{1}{4} - \frac{l(l + 1)}{\rho^2}\big] R(\rho) = 0

To make further progress we consider the behaviour of this differential equation as \rho \rightarrow \infty. It reduces to

\frac{d^2 R}{d \rho^2} - \frac{1}{4} R = 0

which is a straightforward second-order differential equation with auxiliary equation

\zeta^2 - \frac{1}{4} = 0

\implies \zeta = \pm \frac{1}{2}

The positive solution to the auxiliary equation implies a term in the general solution of the form e^{\rho/2} which is unacceptable since it explodes as \rho \rightarrow \infty. Therefore we only accept the negative solution to the auxiliary equation and the general solution for R as \rho \rightarrow \infty must be of the form

R \propto e^{-\rho/2}

This suggests we can try an exact solution of the full differential equation of the form

R = e^{-\rho/2} F(\rho)

Differentiating this twice we get

\frac{d R}{d \rho} = -\frac{1}{2} e^{-\rho/2} F(\rho) + e^{-\rho/2} F^{\prime} (\rho)

\frac{d^2 R}{d \rho^2} = \frac{1}{4} e^{-\rho/2} F(\rho) - \frac{1}{2} e^{-\rho/2} F^{\prime}(\rho) - \frac{1}{2} e^{-\rho/2} F^{\prime}(\rho) + e^{-\rho/2} F^{\prime \prime} (\rho)

Substituting these into the differential equation

\frac{d^2 R}{d \rho^2} + \frac{2}{\rho} \frac{d R}{d \rho} + \big[\frac{\tau}{\rho} - \frac{1}{4} - \frac{l(l + 1)}{\rho^2}\big] R(\rho) = 0

gives

F^{\prime \prime}(\rho) + \frac{(2 - \rho)}{\rho} F^{\prime}(\rho) + \big[\frac{(\tau - 1)}{\rho} - \frac{l(l + 1)}{\rho^2}\big] F(\rho) = 0

\iff

\rho^2 F^{\prime \prime}(\rho) + \rho (2 - \rho) F^{\prime}(\rho) + \big[\rho (\tau - 1) - l(l + 1)\big] F(\rho) = 0

We can now try to solve this latest version of the differential equation by the method of Frobenius, which involves assuming a generalised power series solution of the form

F(\rho) = a_0 \rho^s + a_1 \rho^{s + 1} + a_2 \rho^{s + 2} + \cdots

Differentiating twice we get

F^{\prime}(\rho) = s a_0 \rho^{s-1} + (s + 1) a_1 \rho^s + (s + 2) a_2 \rho^{s + 1} + \cdots

F^{\prime \prime}(\rho) = (s - 1) s a_0 \rho^{s-2} + s (s + 1) a_1 \rho^{s-1} + (s + 1) (s + 2) a_2 \rho^s + \cdots

Then the terms appearing in the differential equation have the generalised power series forms

\rho^2 F^{\prime \prime}(\rho) = (s - 1) s a_0 \rho^s + s (s + 1) a_1 \rho^{s+1} + (s + 1) (s + 2) a_2 \rho^{s+2} + \cdots

2 \rho F^{\prime}(\rho) = 2 s a_0 \rho^s + 2 (s + 1) a_1 \rho^{s+1} + 2 (s + 2) a_2 \rho^{s + 2} + \cdots

-\rho^2 F^{\prime}(\rho) = - s a_0 \rho^{s+1} - (s + 1) a_1 \rho^{s+2} - (s + 2) a_2 \rho^{s + 3} - \cdots

(\tau - 1) \rho F(\rho) = (\tau - 1) a_0 \rho^{s+1} + (\tau - 1) a_1 \rho^{s + 2} + (\tau - 1) a_2 \rho^{s + 3} + \cdots

-l(l + 1) F(\rho) = -l(l + 1) a_0 \rho^s - l(l + 1) a_1 \rho^{s + 1} - l(l + 1) a_2 \rho^{s + 2} - \cdots

Summing these terms (remembering that the sum must be identically equal to zero) we find the coefficient of \rho^s to be

[s(s - 1) + 2s - l(l+1)]a_0 = 0

\implies

s(s + 1) - l(l+1) = 0

\implies

s = l or s = -l - 1

Now, when s = -l - 1 the first term of the power series for F(\rho) is a_0/\rho^{l+1} which explodes as \rho \rightarrow 0. This is unacceptable so we discard this solution and set s = l.

For the coefficient of \rho^{s + n} we get

[(s+n)(s+n-1) + 2(s+n) - l(l+1)]a_n + [(\tau - 1) - (s+n-1)]a_{n-1} = 0

Setting s = l and rearranging gives us the recurrence equation

a_n = \frac{(l + n - \tau)}{(l+n+1)(l+n) - l(l+1)} a_{n-1}

From this recurrence equation we observe that

a_n \rightarrow \frac{1}{n}a_{n-1} = \frac{1}{n!}a_0

as n \rightarrow \infty. We deduce from this that the series for F(\rho) becomes like a_0 \rho^l \sum \frac{\rho^n}{n!} as n \rightarrow \infty and therefore R = e^{-\rho/2} F(\rho) becomes like a_0 \rho^l e^{\rho/2}. However, this diverges as \rho \rightarrow \infty which is unacceptable, so we conclude that the series for F(p) must terminate at some value of n which we will call N. In this case we have a_{N+1} = 0 which the recurrence equation tells us can only happen if

\tau = l + N + 1 \equiv \tilde{n}

This is how the principal quantum number \tilde{n} first appears. Now, we have

\tau = \frac{e^2}{4 \pi \epsilon_0 \hbar} \big(-\frac{m_e}{2 E} \big)^{1/2} = \tilde{n}

\iff

\big(\frac{e^2}{4 \pi \epsilon_0}\big)^2 \big(-\frac{m_e}{2 \hbar^2}\big) \frac{1}{E} = \tilde{n}^2

\iff E_{\tilde{n}} = \big(-\frac{m_e}{2 \hbar^2}\big) \big(\frac{e^2}{4 \pi \epsilon_0}\big)^2 \frac{1}{\tilde{n}^2}

These are the famous bound-state energy eigenvalues for \tilde{n} = 1, 2, \ldots. This is the same formula for the energy levels of the hydrogen atom that Niels Bohr obtained by intuitive means in his 1913 ‘solar system’ model of atomic structure.

As stated above, the integer \tilde{n} is called the principal quantum number. Recall that \tilde{n} = l + N + 1 and N cannot be smaller than zero. It follows that

\tilde{n} - l - 1 \geq 0

\iff

l \leq \tilde{n} - 1

This explains why for a given value of \tilde{n} the allowable values of l are l = 0, 1, 2, \dots, (\tilde{n}-1).

Returning to the solution of

\rho^2 F^{\prime \prime}(\rho) + \rho (2 - \rho) F^{\prime}(\rho) + \big[\rho (\tau - 1) - l(l + 1)\big] F(\rho) = 0

the above discussion suggests that we should look for a solution of the form

F(\rho) = a_0 \rho^l L(\rho)

where L(\rho) is a polynomial (rather than an infinite series). Differentiating this twice gives

F^{\prime}(\rho) = a_0 l \rho^{l-1}L(\rho) + a_0 \rho^l L^{\prime}(\rho)

F^{\prime \prime} (\rho) = a_0 (l-1) l \rho^{l-2}L(\rho) + 2 a_0 l \rho^{l-1} L^{\prime}(\rho) + a_0 \rho^l L^{\prime \prime}(\rho)

Substituting these into the differential equation and setting \tau = \tilde{n} we get

\rho^{l+2} L^{\prime \prime}(\rho) + (2l + 2 - \rho) \rho^{l+1}L^{\prime}(\rho) + (\tilde{n} - 1 - l)\rho^{l+1}L(\rho) = 0

\iff

\rho L^{\prime \prime}(\rho) + (2l + 2 - \rho) L^{\prime}(\rho) + (\tilde{n} - 1 - l) L(\rho) = 0

\iff

\rho L^{\prime \prime}(\rho) + (\alpha + 1 - \rho) L^{\prime}(\rho) + \tilde{n}^{*} L(\rho) = 0

where \alpha \equiv 2l + 1 and \tilde{n}^{*} \equiv \tilde{n} - 1 - l. This last form is a well known differential equation whose physically acceptable solutions in the present context are associated Laguerre polynomials given by the formula

L_{\tilde{n}^{*}}^{(\alpha)} = \sum_{j = 0}^{\tilde{n}^{*}} (-1)^j \frac{(\tilde{n}^{*} + \alpha)!}{(\tilde{n}^{*}-j)!(\alpha + j)!}\frac{\rho^j}{j!}

For given quantum numbers \tilde{n} and l, the solution of the radial equation for R is then

R_{\tilde{n}l}(r) \propto e^{-\rho/2}p^l L_{\tilde{n} - l - 1}^{(2l + 1)}

Final form of the electron wave function \psi

Putting everything together, for given principal quantum number \tilde{n}, orbital quantum number l and magnetic quantum number m_l, the wave function of the electron in the hydrogen atom is

\psi_{\tilde{n} l m_l}(r, \theta, \phi) \propto e^{-\rho/2}p^l L_{\tilde{n} - l - 1}^{(2l + 1)} P_l^{m_l}(\cos \theta) e^{i m_l \phi}

where

\rho = \big(-\frac{8 m_e E_{\tilde{n}}}{\hbar^2}\big)^{1/2} r

and

E_{\tilde{n}} = \big(-\frac{m_e}{2 \hbar^2}\big) \big(\frac{e^2}{4 \pi \epsilon_0}\big)^2 \frac{1}{\tilde{n}^2}