Wsrtjtyk

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

More precisely, the operator

$$-frac{d^2}{dx^2}+x^2$$

can be decomposed as

$$-frac{d^2}{dx^2}+x^2 = left(-frac{d}{dx}-xright)left(frac{d}{dx}-xright)=:a^*a.$$

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

$$left(-frac{d}{dr}-frac{n-1}{r}-rright)left(frac{d}{dr}-rright) = -frac{d^2}{dr^2} -frac{n-1}{r} frac{d}{dr} +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.