I'll make an attempt at providing the steps you are seeking to go "from Euler equation to Bessel function".
You start from the Euler equation, describing conservation of momentum,
$$\rho\frac{\partial \vec{u}}{\partial t}+\rho\vec{u}\cdot\nabla\vec{u}=-\nabla p$$
and the continuity equation, describing conservation of mass,
$$\frac{\partial\rho}{\partial t}=-\nabla\cdot(\rho \vec{u})$$
These are nonlinear equations, to make them tractable you'll want to linearize them, both in the velocity $\vec{u}$ and in the deviations $\delta\rho=\rho-\rho_0$ of the density from the uniform density $\rho_0$. This approximation throws away lots of interesting physics (shock waves, turbulence,...), but without it no simple solution exists.
The linearized equations read
$$\rho_0\frac{\partial \vec{u}}{\partial t}=-\nabla p$$
$$\frac{\partial\delta\rho}{\partial t}=-\rho_0\nabla\cdot\vec{u}$$
We may also assume a linear relation $p=C^2\delta\rho$ between the pressure $p$ and the density variations. (This is a socalled adiabatic equation of state, the coefficient $C^2$ must be positive for mechanical stability.) We define $\xi=\delta\rho/\rho_0$, take the divergence of the first equation and the time derivative of the second equation,
$$\nabla\cdot\frac{\partial \vec{u}}{\partial t}=-C^2\nabla^2 \xi$$
$$\frac{\partial^2\xi}{\partial t^2}=-\frac{\partial}{\partial t}\nabla\cdot\vec{u}$$
Finally, we substitute the first equation into the second one, exchanging the order of differentiation with respect to time and space, to arrive at a wave equation for $\xi$,
$$\frac{\partial^2\xi}{\partial t^2}=C^2\nabla^2 \xi$$
The quantity $C>0$ represents the speed of sound.
We seek a solution of this equation that is a harmonic function of time, so it oscillates with frequency $\omega$. Rather than working with sines or cosines, it is more convenient to use a complex notation, writing
$$\xi(\vec{r},t)={\rm Re}\;e^{-i\omega t}f(\vec{r})$$
The complex function $f$ satisfies the Poisson equation,
$$C^2\nabla^2 f=-\omega^2 f$$
Let's seek a solution with cylindrical symmetry, so $f(R)$ depends only on the radial coordinate $R=\sqrt{x^2+y^2}$. The Poisson equation in cylindrical coordinates takes the form
$$\frac{d^2}{dR^2}f(R)+\frac{1}{R}\frac{d}{dR}f(R)=-(\omega/C)^2f$$
The solution is a Bessel function
$$f(R)={\rm constant}\times J_0(\omega R/C)$$
The full solution thus becomes
$$\delta\rho/\rho_0=A\cos(\omega t+B)J_0(\omega R/C)$$
where $A$ and $B$ are arbitrary coeffients.
And we're done :)