Why do we treat molecular vibrations as linear harmonic oscillations and not 3D isotropic? Born-Oppenheimer approximation leads to this equation for the eigenfunctions and eigenvalues describing molecular vibrations: $$\left[-\frac{\hbar^2}{2\mu}\frac{1}{R^2}\frac{\partial}{\partial R}\left(R^2\frac{\partial}{\partial R}\right)+\frac{\langle\Phi_s|\mathbf{N^2}|\Phi_s\rangle}{2\mu R^2}+E_s(R)-E\right]F_s(R)=0$$
(faithfully copied from page 484, Physics of Atoms and Molecules, Bransden & Joachain 2003). Assuming that we are dealing with a molecule which is not rotating (i.e. $\mathbf{N}=0)$, and approximating the curve $E_s(R)$ with the harmonic potential $\frac{1}{2}k(R-R_0)^2$ in the vicinity of its point of minimum, the B.O. equation above is totally equal to that of a 3D isotropic harmonic oscillator with null angular momentum. Therefore, its solutions should be given by associated Laguerre polynomials times some stuff. Turns out that Hermite polinomials are used instead. Why? The idea that a 3d harmonic oscillator which does not rotate should boil down to a linear harmonic oscillator makes total sense to me. But the equations give these two different polynomials and there's no way to say that $L_n^{1/2}(x^2)\propto H_n(x)$ (namely that associated Laguerre polynomials for null angular momentum case are proportional to Hermite p. of the same degree to make sure that eigenfunctions for 3d and 1d are the same). What am I doing wrong?
 A: For the radial part of the Schrödinger equation for 3D harmonic oscillator with $l=0$,
$$-\frac1{r^2}\frac{\mathrm d}{\mathrm d r}\left(r^2f'(r)\right)+Ar^2 f(r)=Ef(r),$$
one can use the substitution
$$f(r)\to \frac{g(r)}r,$$
which will lead to the equation
$$-g''(r)+Ar^2g(r)=Eg(r),$$
which is isomorphic to the Schrödinger equation of 1D harmonic oscillator. Thus, the solutions of Schrödinger equations for both 1D and (the radial part of) 3D harmonic oscillators are expressible in terms of the same kind of functions.
Interestingly, this works only in number of dimensions $n=1$ or $n=3$ (with generalized substitution $f(r)\to g(r)r^{\frac{1-n}2}$). For any other $n$, you'll get additional terms in the ODE for $g$ that will require Laguerre polynomials instead of Hermite ones.
A: I think Ruslan gave an excellent answer, and I thought I could add a bit of intuition from another perspective. Now, maybe I'm completely wrong about this because currently I don't deal with molecules, but in the field of optics, one encounters the equations that deal with diffraction of a laser beam (Fresnel diffraction) which lead to eigen modes who obey the 2D quantum harmonic oscillator equation.
These can be described with Hermite-Gauss modes when working in cartesian coordinates:
$$I(x,y,z)\propto H_n \left(\frac{x}{\frac{w}{\sqrt{2}}}\right ) H_m \left(\frac{y}{\frac{w}{\sqrt{2}}}\right ) \exp(-\frac{x^2+y^2}{2q^2}) $$
where q,w are parameters related to the size and divergence of the beam and $n,m$ are the quantum numbers of the mode (and z is the direction the beam is propagating in). See image taken from Wikipedia:

and in cylindrical coordinates the modes can be described with Laguerre-Gauss modes, where instead of Hermite polynomials as functions of x,y one encounters a generalized Laguerre polynomial as function of r and another $\exp (i m \varphi)$ term giving the beam "orbital angular momentum". (Image taken from this site)

As you can see, due to "degeneracy", the modes with $l>0$ can be represented both by LG modes and both by HG modes (their both valid bases to describe the propagation of the beam). However, $l=0$ is the same at both situations, because Hermite and Laguerre 0's polynomial are both constants.
What I'm trying to say is: you're doing nothing wrong, just working in a different, equivalent base, which for $l=0$ also give exactly the same modes.
