Consider a particle of mass \(M\) moving in the \(xy\) plane in the presence of the Coulomb potential \[V = -\frac{\alpha}{\sqrt{x^2 + y^2}} \qquad (\alpha > 0) \,.\] Bound states \(\psi\) and energy lelvels \(E < 0\) are determined by the time-independent Schrödinger equation \[-\frac{\hbar^2}{2 M} \nabla^2 \psi + V \psi = E \psi \,.\] In polar coordinates (\(x = r \cos \theta\), \(y = r \sin \theta\)), the equation reads \[-\frac{\hbar^2}{2 M} \left( \frac{\partial^2}{\partial r^2} + \frac{1}{r} \frac{\partial}{\partial r} + \frac{1}{r^2} \frac{\partial^2}{\partial \theta^2} \right) \psi - \frac{\alpha}{r} \psi = E \psi \,.\] Using separation of variables, along with the condition that \(\psi\) is single-valued, we get \[\psi(r,\theta) = R(r) e^{i m \theta} \qquad (m \in \mathbb{Z}) \,,\] where the radial wave function \(R(r)\) satisfies \[\frac{d^2 R}{d r^2} + \frac{1}{r} \frac{d R}{d r} + \left( \frac{2 M E}{\hbar^2} + \frac{2 M \alpha}{\hbar^2 r} - \frac{m^2}{r^2} \right) R = 0 \,.\] Introducing dimensionless (positive) energy \(\mathcal{E}\) and radius \(\rho\) according to \[E = -\frac{M \alpha^2}{2 \hbar^2} \mathcal{E} \qquad (\mathcal{E} > 0)\,,\] \[r = \frac{\hbar^2}{2 M \alpha} \rho \,,\] we obtain \[\frac{d^2 R}{d \rho^2} + \frac{1}{\rho} \frac{d R}{d \rho} + \left( -\frac{\mathcal{E}}{4} + \frac{1}{\rho} - \frac{m^2}{\rho^2} \right) R = 0 \,.\]

Let us first consider the asymptotic behavior of \(R\) at small and large \(\rho\):

Limit \(\rho \to 0\). Suppose \(R \sim \rho^{\nu}\) with \(\nu \ge 0\). Then, the most singular terms in the equation balance each other out if \(\nu\) is such that \(\nu (\nu-1) + \nu - m^2 = 0\), or \(\nu = |m|\). Hence, \[R \sim \rho^{|m|} \,.\]

Limit \(\rho \to \infty\). Suppose \(R \sim e^{-\lambda \rho}\) with \(\lambda > 0\). Since the left-hand side of the equation is dominated by \(R'' - (\mathcal{E}/4) R\), we get \(\lambda^2 - \mathcal{E}/4 = 0\), or \(\lambda = \sqrt{\mathcal{E}}/2\). Hence, \[R \sim e^{-\sqrt{\mathcal{E}} \rho / 2} \,.\]

In view of the asymptotic behavior of \(R\), we make the following substitution: \[R = \rho^{|m|} e^{-\sqrt{\mathcal{E}} \rho / 2} S(\rho) \,,\] where \(S(\rho)\) is a yet unknown function. A straightforward calculation shows that \(S\) must satisfy \[\rho \frac{d^2 S}{d \rho^2} + \Big( 2 |m| + 1 - \sqrt{\mathcal{E}} \rho \Big) \frac{d S}{d \rho} - \left( \frac{2 |m| + 1}{2} \sqrt{\mathcal{E}} - 1 \right) S = 0 \,.\] Introducing \[z = \sqrt{\mathcal{E}} \rho \,,\] we find that \(S\) is described by Kummer's equation \[z \frac{d^2 S}{d z^2} + (b - z) \frac{d S}{d z} - a S = 0\] with \[a = |m| + \frac{1}{2} - \frac{1}{\sqrt{\mathcal{E}}} \,,\] \[b = 2 |m| + 1 \,.\] The solution (nonsingular at the origin) is given by Kummer's confluent hypergeometric function \[M(a, b, z) = \sum_{k=0}^{\infty} \frac{a^{(k)}}{b^{(k)} k!} z^k = 1 + \frac{a}{b} z + \frac{a (a+1)}{b (b+1) 2!} z^2 + \ldots \,,\] where \(a^{(k)}\) and \(b^{(k)}\) are the rising factorials defined as \[r^{(0)} = 1 \,, \qquad r^{(1)} = r \,, \qquad r^{(2)} = r (r + 1) \,, \qquad \ldots\] \[r^{(k)} = r (r + 1) \ldots (r + k - 1) = \frac{\Gamma(r+k)}{\Gamma(r)} \,.\] Except when \(a = 0, -1, -2, \ldots \,,\) \[M(a,b,z) \sim \frac{z^{a-b} e^z}{\Gamma(a)} \qquad \text{as} \qquad z \to \infty \,,\] implying \(R \to \infty\), which is physically impossible. When \(a = -n\), with \(n = 0, 1, 2, \ldots \,,\) function \(M(-n,b,z)\) is a polynomial in \(z\) of degree not exceeding \(n\). In this case, \(R\) decays to zero at infinity sufficiently fast and is normalizable. This yields the following quantization codition: \[|m| + \frac{1}{2} - \frac{1}{\sqrt{\mathcal{E}}} = -n \,,\] or \[\mathcal{E} = \frac{1}{\left( n + |m| + \frac{1}{2} \right)^2} \qquad (n = 0, 1, 2, \ldots; \; m = 0, \pm 1, \pm 2, \ldots) \,.\] The corresponding (unnormalized) radial eigenfunction is given by \[R = \rho^{|m|} e^{-\sqrt{\mathcal{E}} \rho / 2} M(-n, \, 2 |m| + 1, \, \sqrt{\mathcal{E}} \rho) \,,\] or \[R = \rho^{|m|} \exp \left[ -\frac{\rho}{2 n + 2 |m| + 1} \right] M \left( -n, \, 2 |m| + 1, \, \frac{2 \rho}{2 n + 2 |m| + 1} \right) \,.\]

An alternative (and more conventional) way to label the bound states is to use quantum numbers \[N \equiv n + |m| + 1 = 1, 2, 3, \ldots\] and \[m = -(N-1), \, \ldots \,, N-1 \,.\] Then, \[\mathcal{E} = \frac{1}{\left( N - \frac{1}{2} \right)^2} \qquad (N = 1, 2, 3, \ldots)\] and \[R = \rho^{|m|} \exp \left[ -\frac{\rho}{2 N - 1} \right] M \left( -N + |m| + 1, \, 2 |m| + 1, \, \frac{2 \rho}{2 N - 1} \right) \,.\]

**Reference:** Analytic solution of a two-dimensional hydrogen atom. I. Nonrelativistic theory