The paper studies spectral theory of Schrödinger operators H=ℏ2Δ+V on the sphere from the standpoint of integrability and separation. Our goal is to uncover the fine structure of spec H, i.e., asymptotics of eigenvalues and spectral clusters, determine their relation to the underlying geometry and classical dynamics and apply this data to the inverse spectral problem on the sphere. The prototype model is the celebrated Neumann Hamiltonian p2+V with quadratic potential V on Sn. We show that the quantum Neumann Hamiltonian (Schrödinger operator H) remains an integrable and find an explicit set of commuting integrals. We also exhibit large classes of separable potentials {V} based on ellipsoidal coordinates on Sn. Several approaches to spectral theory of such Hamiltonians are outlined. The semiclassical problem (small ℏ) involves the EKB(M)-quantization of the classical Neumann flow along with its invariant tori, Maslov indices, etc., all made explicit via separation of variables. Another approach exploits Stäckel–Robertson separation of the quantum Hamiltonian and reduction to certain ODE problems: the Hill’s and the generalized Lamè equations. The detailed analysis is carried out for S2, where the ODE becomes the perturbed classical Lamè equation and the Schrödinger eigenvalues are expressed through the Lamè eigendata.