We consider the quasi-static problem governing the localized surface plasmon modes and permittivity eigenvalues ǫ of smooth, arbitrarily shaped, axisymmetric inclusions. We develop an asymptotic theory for the dense part of the spectrum, i.e., close to the accumulation value ǫ = −1 at which a flat interface supports surface plasmons; in this regime, the field oscillates rapidly along the surface and decays exponentially away from it on a comparable scale. With τ = −(ǫ + 1) as the small parameter, we develop a surface-ray description of the eigenfunctions in a narrow boundary layer about the interface; the fast phase variation, as well as the slowly varying amplitude and geometric phase, along the rays are determined as functions of the local geometry. We focus on modes varying at most moderately in the azimuthal direction, in which case the surface rays are meridian arcs that focus at the two poles. Asymptotically matching the diverging ray solutions with expansions valid in inner regions in the vicinities of the poles yields the quantization rulewhere n ≫ 1 is an integer and Θ a geometric parameter given by the product of the inclusion length and the reciprocal average of its cross-sectional radius along its symmetry axis. For a sphere, Θ = π, whereby the formula returns the exact eigenvalues ǫ n = −1 − 1/n. We also demonstrate good agreement with exact solutions in the case of prolate spheroids.