I look to solve analytically the Vlasov-Maxwell equations for a magnetized hot plasma $$\frac{\partial f_{s1}}{\partial t}+\vec{v}\cdot\frac{\partial f_{s1}}{\partial \vec{r}}+\frac{q_s}{m_s}\Big(\vec{v}\wedge\vec{B}_0 \Big)\cdot\frac{\partial f_{s1}}{\partial \vec{v}} =-\frac{q_s}{m_s}\Big(\vec{E}_1+\vec{v}\wedge\vec{B}_1 \Big)\cdot \frac{\partial f_{s0}}{\partial \vec{v}}$$ $$\vec{\nabla}\cdot\vec{E}_1=\dfrac{1}{\varepsilon_0}\displaystyle\sum_s q_s \int f_{s1}\,\mathrm{d^3}v$$ $$\vec{\nabla}\cdot\vec{B}_1=0$$ $$\vec{\nabla}\wedge\vec{E}_1=-\frac{\partial \vec{B}_1}{\partial t}$$ $$ \vec{\nabla}\wedge\vec{B}_1=\mu_0 \displaystyle\sum_s q_s \int \vec{v}\; f_{s1}\,\mathrm{d^3}v+\frac{1}{c^2}\frac{\partial \vec{E}_1}{\partial t}$$
With $f_{s0}(\vec{v})\equiv f_{s0}(v_\perp,v_z)=$ the zero order distribution function in cylindrical geometry.
$(f_{s1}, \vec{E}_1, \vec{B}_1)=$ the perturbed parameters of order 1 , $B_{0}=$ the constant magnetic field.
In order to calculate the perturbed distribution function $f_{s1}$, resolution of Vlasov equation is done by the method of characteristics with a cylindrical geometry which involves very difficult calculations. (see Swanson p.91-95)
Instead of this method, other works use the development of $f_{s1}$ on the spherical harmonics $$f_{s1}= \displaystyle\sum_l\sum_m f_l^m(v)\;Y_l^m (\theta,\varphi)$$ then by projection of the $Y_l^m$ on Vlasov equation and the relation of orthogonality of $Y_l^m$ they form an infinit system of coupled equations in $f_l^m$.
So I wanted to apply this method of spherical harmonics in this problem but I couldn't solve the system of coupled equations in $f_l^m$. I therefore concluded that spherical harmonics cannot be applied.
Please, do you know of any other alternative methods to calculate the perturbed distribution function $f_{s1}$ in the same cylindrical geometry?