Solution (source code)

= Solution

Write total <pressure> and <gravitational potential> as
$$
p=p_0(r)+\pi(r,z,t),\qquad \Phi=\Phi_0(r)+\psi(r,z,t).
$$
Here $\pi$ is the <Eulerian perturbation of a fluid variable> applied to <pressure>, while $\psi$ is the corresponding <gravitational potential> perturbation. They have units $\mathrm{Pa}$ and $\mathrm{m^2\,s^{-2}}$, respectively. The constant interior <mass density> is $\rho$.

The <Euler equations> are
$$
\partial_t\mathbf u+(\mathbf u\cdot\nabla)\mathbf u
=-\frac1\rho\nabla(p_0+\pi)-\nabla(\Phi_0+\psi).
$$
Subtract the <hydrostatic equilibrium> relation $\nabla p_0/\rho+\nabla\Phi_0=0$ and discard terms quadratic in the perturbation. <Mass conservation> gives <incompressible flow>, and the gravitational <Poisson equation> has no bulk density perturbation. Thus, strictly inside the fluid,
$$
\boxed{\partial_t\mathbf u=-\frac{\nabla\pi}{\rho}-\nabla\psi,\qquad
\nabla\cdot\mathbf u=0,\qquad\nabla^2\psi=0.}
$$
The last equation is <Laplace equation>, not a statement that self-gravity has disappeared. The displaced surface carries a first-order <surface mass density> $\rho\eta$, where $\eta$ is its radial displacement, and this supplies the source and matching conditions for $\psi$ at the moving interface. The <axisymmetric gravitational potential of a displaced cylindrical interface>, expressed using <modified Bessel functions>, incorporates that boundary source.

It is useful to define $W=\psi+\pi/\rho$. Taking the <divergence> of the momentum equation gives $\nabla^2W=0$. For an axisymmetric <Fourier mode> with axial <wavenumber> $k>0$, regularity at $r=0$ selects radial dependence $I_0(kr)$, the <Modified Bessel function of the first kind>, since $K_0(kr)$ is singular there. In what follows $k=|k_z|$ can always be chosen nonnegative.