My approach is this: let's look for a family of special functions that carry a nice representation of $E(2)$, the Euclidean motions. It's easy to handle the rotations. Any family of functions of the form $f_n(x,y) = b_n(r)e^{in\theta}$ behaves nicely under rotations. So now we look for a suitable choice of the $b_n$ that behaves nicely under translations too. As $E(2)$ isn't commutative we can't expect anything as simple as we did with rotations. We can look at the action of the generators of translations in $E(2)$, ie. the derivatives of the $f_n$ w.r.t. $x$ and $y$.

We get

\begin{align} \frac{\partial f_n}{\partial x} &= \frac{e^{in\theta}}{r}\Big(-\frac{inyb_n(r)}{r}+xb_n'(r)\Big) \label{partialx}\\ \frac{\partial f_n}{\partial y} &= \frac{e^{in\theta}}{r}\Big(\frac{inxb_n(r)}{r}+yb_n'(r)\Big) \label{partialy} \end{align}

With $z=x+iy$ these give

\begin{align*} 2\frac{\partial f_n}{\partial\bar{z}} & = e^{i(n+1)\theta} \Big(-\frac{nb_n(r)}{r}+b_n'(r)\Big)\\ 2\frac{\partial f_n}{\partial z} & = e^{i(n-1)\theta} \Big(\frac{nb_n(r)}{r}+b_n'(r)\Big) \end{align*}

If we pick $$ -\frac{nb_n(r)}{r}+b_n'(r)=0 $$ our $f_n$ can be chosen to be just the polynomials in $z$. We know these carry a nice enough representation of translations and rotations. The effect of translations is given by the binomial theorem.

But looking at the way $\partial/\partial z$ and $\partial/\partial \bar{z}$ raise and lower the $n$ in the complex exponent, there's another path begging to be taken. Choose \begin{align*} b_{n+1}(r) & = \frac{nb_n(r)}{r}-b_n'(r) \\ b_{n-1}(r) & = \frac{nb_n(r)}{r}+b_n'(r) \end{align*} and we get the pretty \begin{align} \frac{2nb_n(r)}{r} & = b_{n-1}(r)+b_{n+1}(r) \label{sum} \\ 2b_n'(r) & = b_{n-1}(r)-b_{n+1}(r) \label{deriv} \end{align} If we view the second relation as a family of couple ODEs we can choose any initial conditions we want. So pick a "unit" vector $b_n(0)=\delta_{n0}$. The resulting $b_n$ are the Bessel functions $B_n$.

So by looking for functions that behave nicely under translations and rotations we're led to the Bessel functions as a particularly nice choice. (The actual behaviour under finite translations is a little bit more complicated than the derivative and is given by Graf's addition theorem.)

(This is part of the content to the references you link to but I've tried to unpack some details and make it more elementary.)

My approach is this: let's look for a family of special functions that carry a nice representation of $E(2)$, the Euclidean motions. It's easy to handle the rotations. Any family of functions of the form $f_n(x,y) = b_n(r)e^{in\theta}$ behaves nicely under rotations. So now we look for a suitable choice of the $b_n$ that behaves nicely under translations too. As $E(2)$ isn't commutative we can't expect anything as simple as we did with rotations. We can look at the action of the generators of translations in $E(2)$, ie. the derivatives of the $f_n$ w.r.t. $x$ and $y$.

We get

\begin{align} \frac{\partial f_n}{\partial x} &= \frac{e^{in\theta}}{r}\Big(-\frac{inyb_n(r)}{r}+xb_n'(r)\Big) \label{partialx}\\ \frac{\partial f_n}{\partial y} &= \frac{e^{in\theta}}{r}\Big(\frac{inxb_n(r)}{r}+yb_n'(r)\Big) \label{partialy} \end{align}

With $z=x+iy$ these give

\begin{align*} 2\frac{\partial f_n}{\partial\bar{z}} & = e^{i(n+1)\theta} \Big(-\frac{nb_n(r)}{r}+b_n'(r)\Big)\\ 2\frac{\partial f_n}{\partial z} & = e^{i(n-1)\theta} \Big(\frac{nb_n(r)}{r}+b_n'(r)\Big) \end{align*}

If we pick $$ -\frac{nb_n(r)}{r}+b_n'(r)=0 $$ our $f_n$ can be chosen to be just the monomials in $z$. We know these carry a nice enough representation of translations and rotations. The effect of translations is given by the binomial theorem.

But looking at the way $\partial/\partial z$ and $\partial/\partial \bar{z}$ raise and lower the $n$ in the complex exponent, there's another path begging to be taken. Choose \begin{align*} b_{n+1}(r) & = \frac{nb_n(r)}{r}-b_n'(r) \\ b_{n-1}(r) & = \frac{nb_n(r)}{r}+b_n'(r) \end{align*} and we get the pretty \begin{align} \frac{2nb_n(r)}{r} & = b_{n-1}(r)+b_{n+1}(r) \label{sum} \\ 2b_n'(r) & = b_{n-1}(r)-b_{n+1}(r) \label{deriv} \end{align} If we view the second relation as a family of couple ODEs we can choose any initial conditions we want. So pick a "unit" vector $b_n(0)=\delta_{n0}$. The resulting $b_n$ are the Bessel functions $B_n$.

So by looking for functions that behave nicely under translations and rotations we're led to the Bessel functions as a particularly nice choice. (The actual behaviour under finite translations is a little bit more complicated than the derivative and is given by Graf's addition theorem.)

(This material is part of the content in the references you link to but I've tried to unpack some details and make it more elementary.)