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 :)