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

In 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)$

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}$