CBE has been a challenging problem in various disciplines of physical sciences. Due to existing difficulties of the general problem, finding a solution to the linearized CBE became the center of attraction in the twentieth century when @L46 and van Kampen (1955) discovered the normal modes of collisionless ensembles. Later in 1970’s, @K71 [@K77] developed a matrix theory that was capable of computing normal modes of stellar systems through solving a nonlinear eigenvalue problem. His theory remained as the only analytical perturbation theory used by the community of galactic dynamicists over the past three decades.

@K77 assumed an exponential form $\exp(-{\rm i}\omega t)$ for the time-dependent part of physical quantities where ${\rm i}=\sqrt{-1}$, and solved the linearized CBE for the perturbed distribution function (DF) $f_1$ in terms of the perturbed potential $V_1$. After expanding the potential and density functions in terms of bi-orthogonal basis sets in the configuration space, he used the weighted residual form of the fundamental equation $$f_1 d \textbf{\textit{v}} d \textbf{\textit{x}}=
\Sigma_1 d \textbf{\textit{x}},$$ to obtain a nonlinear eigenvalue problem for the complex eigenfrequency $\omega$. For self-consistent perturbations the surface density $\Sigma_1=\int f_1 d \textbf{\textit{v}}$ is related to $V_1$ through Poisson’s integral, and the symbols $ d 