Since several presentations on this website make use of the radial Schrödinger equation in spherical coordinates for an arbitrary number of spatial dimensions, it is long past time to provide an explicit derivation. With quantum mechanical linear momentum defined by the gradient $\mathbf{P}=\frac{\hslash}{i}{\nabla}_{\left(n\right)}$ , the Schrödinger equation in general coordinates is

$[\frac{{\mathbf{P}}^{2}}{2m}+V\left({x}_{k}\right)]\psi \left({x}_{k}\right)=[-\frac{{\hslash}^{2}}{2m}{\nabla}_{\left(n\right)}^{2}+V\left({x}_{k}\right)]\psi \left({x}_{k}\right)=E\psi \left({x}_{k}\right)$

where the functional argument implies all coordinates. The Laplacian operator in Cartesian coordinates generalizes easily to higher dimensions,

${\nabla}_{\left(n\right)}^{2}=\frac{{\partial}^{2}}{\partial {x}_{1}^{2}}+\frac{{\partial}^{2}}{\partial {x}_{2}^{2}}+\cdots +\frac{{\partial}^{2}}{\partial {x}_{n}^{2}}=\sum _{k=1}^{n}\frac{{\partial}^{2}}{\partial {x}_{k}^{2}}$

but curvilinear coordinates are not so simple. Consider first two-dimensional coordinates given by

${x}_{1}={r}_{\left(2\right)}sin\phi \phantom{\rule{5em}{0ex}}{x}_{2}={r}_{\left(2\right)}cos\phi $

which are easily inverted to

${r}_{\left(2\right)}=\sqrt{{x}_{1}^{2}+{x}_{2}^{2}}\phantom{\rule{5em}{0ex}}\phi ={tan}^{-1}\left(\frac{{x}_{1}}{{x}_{2}}\right)$

These latter forms must be used in applying the chain rule to partial derivatives in order to keep track of which variables remain constant in the process. With standard derivatives one has

$\begin{array}{l}\frac{\partial}{\partial {x}_{1}}=\frac{\partial {r}_{\left(2\right)}}{\partial {x}_{1}}\frac{\partial}{\partial {r}_{\left(2\right)}}+\frac{\partial \phi}{\partial {x}_{1}}\frac{\partial}{\partial \phi}\\ \phantom{\frac{\partial}{\partial {x}_{1}}}=\frac{{x}_{1}}{{r}_{\left(2\right)}}\frac{\partial}{\partial {r}_{\left(2\right)}}+\frac{\frac{1}{{x}_{2}}}{1+\frac{{x}_{1}^{2}}{{x}_{2}^{2}}}\frac{\partial}{\partial \phi}\\ \frac{\partial}{\partial {x}_{1}}=sin\phi \frac{\partial}{\partial {r}_{\left(2\right)}}+\frac{cos\phi}{{r}_{\left(2\right)}}\frac{\partial}{\partial \phi}\end{array}$

$\begin{array}{l}\frac{\partial}{\partial {x}_{2}}=\frac{\partial {r}_{\left(2\right)}}{\partial {x}_{2}}\frac{\partial}{\partial {r}_{\left(2\right)}}+\frac{\partial \phi}{\partial {x}_{2}}\frac{\partial}{\partial \phi}\\ \phantom{\frac{\partial}{\partial {x}_{2}}}=\frac{{x}_{2}}{{r}_{\left(2\right)}}\frac{\partial}{\partial {r}_{\left(2\right)}}-\frac{\frac{{x}_{1}}{{x}_{2}^{2}}}{1+\frac{{x}_{1}^{2}}{{x}_{2}^{2}}}\frac{\partial}{\partial \phi}\\ \frac{\partial}{\partial {x}_{2}}=cos\phi \frac{\partial}{\partial {r}_{\left(2\right)}}-\frac{sin\phi}{{r}_{\left(2\right)}}\frac{\partial}{\partial \phi}\end{array}$

Squaring and adding together, remembering to include terms from the left partial derivative acting on the right, leads to

${\nabla}_{\left(2\right)}^{2}=\frac{{\partial}^{2}}{\partial {r}_{\left(2\right)}}+\frac{1}{{r}_{\left(2\right)}}\frac{\partial}{\partial {r}_{\left(2\right)}}+\frac{1}{{r}_{\left(2\right)}^{2}}\frac{\partial}{\partial \phi}$

For more that two dimensions one can then assume a form

${\nabla}_{\left(n\right)}^{2}=\frac{{\partial}^{2}}{\partial {r}_{\left(n\right)}^{2}}+\frac{m\left(n\right)}{{r}_{\left(n\right)}}\frac{\partial}{\partial {r}_{\left(n\right)}}+\frac{1}{{r}_{\left(n\right)}^{2}}{\nabla}_{{\Omega}_{\left(n\right)}}^{2}$

where the integer $m\left(n\right)$ and the additional angular dependence are to be determined. To this end consider a recursion relation

${\nabla}_{(n+1)}^{2}={\nabla}_{\left(n\right)}^{2}+\frac{{\partial}^{2}}{\partial {x}_{n+1}^{2}}$

The additional Cartesian variable will be the product of a cosine with the higher-dimensional radius, while the original radial variable will be the product of a sine with the new radial variable:

${r}_{\left(n\right)}={r}_{(n+1)}sin\theta \phantom{\rule{5em}{0ex}}{x}_{n+1}={r}_{(n+1)}cos\theta $

These must again be inverted to apply the chain rule,

${r}_{(n+1)}=\sqrt{{r}_{\left(n\right)}^{2}+{x}_{n+1}^{2}}\phantom{\rule{5em}{0ex}}\theta ={tan}^{-1}\left[\frac{{r}_{\left(n\right)}}{{x}_{n+1}}\right]$

and since the evaluation is the same as before with different variable names one has immediately

$\begin{array}{l}\frac{\partial}{\partial {r}_{\left(n\right)}}=sin\theta \frac{\partial}{\partial {r}_{(n+1)}}+\frac{cos\theta}{{r}_{(n+1)}}\frac{\partial}{\partial \theta}\\ \\ \frac{\partial}{\partial {x}_{n+1}}=cos\theta \frac{\partial}{\partial {r}_{(n+1)}}-\frac{sin\theta}{{r}_{(n+1)}}\frac{\partial}{\partial \theta}\end{array}$

Squaring and adding these terms proceeds as before, and the original radial variable dividing terms in the assumed form needs to be replaced with its higher-dimensional version. The recursion relation becomes

$\begin{array}{l}{\nabla}_{(n+1)}^{2}=\frac{{\partial}^{2}}{\partial {r}_{(n+1)}^{2}}+\frac{m\left(n\right)+1}{{r}_{(n+1)}}\frac{\partial}{\partial {r}_{(n+1)}}\\ \phantom{\rule{6em}{0ex}}+\frac{1}{{r}_{(n+1)}^{2}}[\frac{{\partial}^{2}}{\partial {\theta}^{2}}+m\left(n\right)\frac{cos\theta}{sin\theta}\frac{\partial}{\partial \theta}+\frac{1}{{sin}^{2}\theta}{\nabla}_{{\Omega}_{\left(n\right)}}^{2}]\end{array}$

The integer on the second term increases by one with each iteration: beginning with $m\left(2\right)=1$ in general one has $m\left(n\right)=n-1$ . The coefficient of the cotangent is retained from the previous step of the iteration, starting from $m(n-1)=n-2$ and decreasing to zero for the last angular variable. The form of the Laplacian in spherical coordinates in an arbitrary number of dimensions is thus

$\begin{array}{l}{\nabla}_{\left(n\right)}^{2}=\frac{{\partial}^{2}}{\partial {r}_{\left(n\right)}^{2}}+\frac{(n-1)}{{r}_{\left(n\right)}}\frac{\partial}{\partial {r}_{\left(n\right)}}+\frac{1}{{r}_{\left(n\right)}^{2}}{\nabla}_{{\Omega}_{\left(n\right)}}^{2}\\ \\ {\nabla}_{{\Omega}_{\left(n\right)}}^{2}=\sum _{k=2}^{n}\frac{1}{\prod _{l=2}^{k-1}{sin}^{2}{\theta}_{l}}[\frac{{\partial}^{2}}{\partial {\theta}_{k}^{2}}+(n-k)\frac{cos{\theta}_{k}}{sin{\theta}_{k}}\frac{\partial}{\partial {\theta}_{k}}]\end{array}$

where an empty product is taken to be unity. With the exception of ${\theta}_{n}$ , the domain of the angular variables is $[\mathrm{0,\pi ]}$ . The domain of ${\theta}_{n}$ is $[\mathrm{0,2\pi ]}$ , since along with the radial variable it parametrizes two Cartesian variables.

This Laplacian can now be inserted in the equation at the outset for a general answer. To determine the radial equation, one needs eigenfunctions of the angular part of the Laplacian. Write the entire wavefunction as the separable product

$\psi \left({x}_{k}\right)=R\left({r}_{\left(n\right)}\right)F\left({\theta}_{k}\right)=R\left({r}_{\left(n\right)}\right)\prod _{k=2}^{n}{f}_{k}\left({\theta}_{k}\right)$

with the functional dependence made explicit in the product. Normalization constants will be omitted in what follows for simplicity, but are understood to be necessary for application.

First look at the overall structure of the angular part of the Laplacian. Subscripted partials here represent the two bracketed terms for almost every angular variable except the last:

${\nabla}_{{\Omega}_{\left(n\right)}}^{2}=[{\partial}_{2}+\frac{1}{{sin}^{2}{\theta}_{2}}[{\partial}_{3}+\frac{1}{{sin}^{2}{\theta}_{3}}[{\partial}_{4}+\frac{1}{{sin}^{2}{\theta}_{4}}[\cdots \left]\right]\left]\right]$

At each stage of the expression the quantities in brackets depend on only one angular variable, as long as all enclosed brackets have constant values. This can be accomplished by setting each bracketed quantity in turn proportional to the relevant function to form an eigenvalue equation. Working downwards from the highest-numbered variable, all enclosed brackets will then be constant.

The highest angular eigenfunction has a simple eigenvalue equation,

$\frac{{\partial}^{2}}{\partial {\theta}_{n}^{2}}{f}_{n}={c}_{n}{f}_{n}\phantom{\rule{2em}{0ex}}\to \phantom{\rule{2em}{0ex}}{f}_{n}={e}^{i{l}_{n}{\theta}_{n}}$

where the eigenvalue is taken to be negative to produce a normalizable eigenfunction. In three dimensions this angular momentum quantum number is normally written with an *m*, but will be denoted here with an angular index for generality.

The second highest angular eigenfunction has a more complicated eigenvalue equation:

$[\frac{{\partial}^{2}}{\partial {\theta}_{n-1}^{2}}+\frac{cos{\theta}_{n-1}}{sin{\theta}_{n-1}}\frac{\partial}{\partial {\theta}_{n-1}}-\frac{{l}_{n}^{2}}{{sin}^{2}{\theta}_{n-1}}]{f}_{n-1}={c}_{n-1}{f}_{n-1}$

If the eigenvalue is taken to be negative and chosen appropriately, then this is the differential equation for associated Gegenbauer polynomials in trigonometric form for $\alpha =\frac{1}{2}$ :

$\begin{array}{c}[\frac{{\partial}^{2}}{\partial {\theta}_{n-1}^{2}}+\frac{cos{\theta}_{n-1}}{sin{\theta}_{n-1}}\frac{\partial}{\partial {\theta}_{n-1}}+{l}_{n-1}({l}_{n-1}+1)-\frac{{l}_{n}^{2}}{{sin}^{2}{\theta}_{n-1}}]{f}_{n-1}=0\\ \\ \to \phantom{\rule{2em}{0ex}}{f}_{n-1}={C}_{{l}_{n-1},{l}_{n}}^{(1/2)}\left({\theta}_{n-1}\right)\end{array}$

This equation is better known as that of the associated Legendre polynomials, since for $\alpha =\frac{1}{2}$ the two types of polynomials are equivalent. For generality only Gegenbauer polynomials will be notated.

The third highest angular eigenfunction has an even somewhat more complicated eigenvalue equation,

$[\frac{{\partial}^{2}}{\partial {\theta}_{n-2}^{2}}+2\frac{cos{\theta}_{n-2}}{sin{\theta}_{n-2}}\frac{\partial}{\partial {\theta}_{n-2}}-\frac{{l}_{n-1}({l}_{n-1}+1)}{{sin}^{2}{\theta}_{n-2}}]{f}_{n-2}={c}_{n-2}{f}_{n-2}$

but if the eigenvalue is again taken to be negative and chosen appropriately, this is still the differential equation for associated Gegenbauer polynomials in trigonometric form, but with $\alpha =1$ :

$\begin{array}{c}[\frac{{\partial}^{2}}{\partial {\theta}_{n-2}^{2}}+2\frac{cos{\theta}_{n-2}}{sin{\theta}_{n-2}}\frac{\partial}{\partial {\theta}_{n-2}}+{l}_{n-2}({l}_{n-2}+2)-\frac{{l}_{n-1}({l}_{n-1}+1)}{{sin}^{2}{\theta}_{n-2}}]{f}_{n-2}=0\\ \\ \to \phantom{\rule{2em}{0ex}}{f}_{n-2}={C}_{{l}_{n-2},{l}_{n-1}}^{\left(1\right)}\left({\theta}_{n-1}\right)\end{array}$

And now the general rule becomes apparent: for the eigenvalue equation

$[\frac{{\partial}^{2}}{\partial {\theta}_{n-m}^{2}}+m\frac{cos{\theta}_{n-m}}{sin{\theta}_{n-m}}\frac{\partial}{\partial {\theta}_{n-m}}-\frac{{l}_{n-m+1}({l}_{n-m+1}+m-1)}{{sin}^{2}{\theta}_{n-k}}]{f}_{n-m}={c}_{n-m}{f}_{n-m}$

if the eigenvalue is always taken to be negative and chosen to produce the equation

$\begin{array}{l}[\frac{{\partial}^{2}}{\partial {\theta}_{n-m}^{2}}+m\frac{cos{\theta}_{n-m}}{sin{\theta}_{n-m}}\frac{\partial}{\partial {\theta}_{n-m}}\\ \phantom{\rule{4em}{0ex}}+{l}_{n-m}({l}_{n-m}+m)-\frac{{l}_{n-m+1}({l}_{n-m+1}+m-1)}{{sin}^{2}{\theta}_{n-k}}]{f}_{n-m}=0\end{array}$

then the solution is an associated Gegenbauer polynomial for $\alpha =\frac{m}{2}$ :

${f}_{n-m}={C}_{{l}_{n-m},{l}_{n-m+1}}^{(m/2)}\left({\theta}_{n-m}\right)$

Working downwards from the highest angular variable, values for the index $\alpha $ increase by one half, leading to a linear increase by one in the second factor of the eigenvalue. This is precisely the behavior of associated Gegenbauer polynomials.

Setting $k=n-m$ in this general solution, the complete unnormalized separable wavefunction is

$\psi \left({x}_{k}\right)=R\left({r}_{\left(n\right)}\right){e}^{i{l}_{n}{\theta}_{n}}\prod _{k=2}^{n-1}{C}_{{l}_{k},{l}_{k+1}}^{\left(\right(n-k)/2)}\left({\theta}_{k}\right)$

The eigenvalue of the first angular variable, $-{l}_{2}({l}_{2}+n-2)$ , is the value that will appear in the radial Schrödinger equation. The quantum number in this eigenvalue limits all other quantum numbers for higher angular variables due to the definition of associated Gegenbauer polynomials as derivatives: for integral indices, the second lower index cannot exceed the first if one is to produce a nonzero function.

Dropping the index two in the eigenvalue as no longer needed, as well the the dimension index on the radial variable, the equation sought is

$[\frac{{d}^{2}}{d{r}^{2}}+\frac{(n-1)}{r}\frac{d}{dr}-\frac{l(l+n-2)}{{r}^{2}}+\frac{2m}{{\hslash}^{2}}[E-V(r\left)\right]\phantom{\rule{.2em}{0ex}}]R\left(r\right)=0$

as has appeared in previous presentations.