mjhjmtu

It is probably the most well-known result in quantum mechanics that the harmonic oscillator can be solved by supersymmetry.

More precisely, the operator

$$-fracd^2dx^2+x^2$$

can be decomposed as

$$-fracd^2dx^2+x^2 = left(-fracddx-xright)left(fracddx-xright)=:a^*a.$$

When I tried to see whether the same holds true in spherical coordinates I noticed a complete suprise

$$left(-fracddr-fracn-1r-rright)left(fracddr-rright) = -fracd^2dr^2 -fracn-1r fracddr +r^2+(n-1).$$

The first three terms comprise the harmonic oscillator in spherical coordinates that is

$$-Delta_r + r^2.$$

For convenience, I discarded the angular part in this calculation.

But for some reasons I seem to get this superficial term $(n-1)$. Why is that? And can I do anything smarter to avoid this term?

Indeed, the supersymmetric operators do not factorize the Hamiltonian of the three-dimensional harmonic oscillator, there is an additional term. See Creation and annihilation operators, symmetry and supersymmetry of the 3D isotropic harmonic oscillator, equation 16.

Supplementing @CarloBeenakker's answer: a different sort of "factorization" of the Laplacian can be done in terms of the Dirac operator, which in some regards better preserves the nature of the situation, but/and has scalar coefficients in an appropriate Clifford algebra, which in general is not commutative. E.g., on $mathbb R$, it is $mathbb C$, which is convenient. But already in two dimensions the relevant Clifford algebra is the Hamiltonian quaternions...

Depending what you want to do, this might be a useful direction to consider, in addition to the "scalarization" of immediate reduction to the spherically symmetric case.