The Bessel and Neumann (or equivalently the Hankel) functions provide a basis for the radial equation, when separating the 2D Laplace equation in polar coordinates. For the angular part, which I'd identify with the analog of spherical harmonics in 2D are simply the orthonormal set of exponential functions
[tex]u_m(\varphi)=\frac{1}{\sqrt{2 \pi}} \exp(\mathrm{i} m \varphi), \quad m \in \mathbb{Z},[/tex]
or, if you prefer real basis,
[tex]u_0(\varphi)=\frac{1}{\sqrt{2 \pi}}, \quad u_m^{(1)}(\varphi)=\cos(m \varphi), \quad u_m^{(2)}(\varphi) \sin(\varphi), \quad m \in \mathbb{N}_{>0}.[/tex]
The general solution of the Laplace equation in terms of the corresponding series reads
[tex]\phi(r,\varphi)=\sum_{m=-\infty}^{\infty} [\phi_m^{(1)} J_m(r/r_0) + \phi_m^{(2)} N_m(r/r_0)] u_m(\varphi).[/tex]
This form with the Bessel and Neumann functions is convenient since [itex]J_m[/itex] is the solution of the radial equation which is analytic in [itex]r=0[/itex], while [itex]N_m[/itex] is singular at the origin.