Solution (source code)

= Solution

For the positive-frequency <internal gravity wave> branch, $\omega=NK_h/K$. Differentiating with respect to the components of the <wave vector> gives the <group velocity>,
$$
\boxed{\mathbf c_g=\left(\frac{Nkm^2}{K_hK^3},\frac{Nlm^2}{K_hK^3},-\frac{NK_hm}{K^3}\right).}
$$
The wavefront-normal <phase velocity> is $\mathbf c_p=\omega\mathbf k/K^2$. Directly,
$$
\mathbf k\cdot\mathbf c_g
=\frac{Nm^2(k^2+l^2)}{K_hK^3}-\frac{NK_hm^2}{K^3}=0,
\qquad\boxed{\mathbf c_p\cdot\mathbf c_g=0.}
$$
The negative-frequency branch changes both propagation signs but not their orthogonality. Another proof is that $\omega$ is homogeneous of degree zero in the <wave vector>, so Euler's homogeneous-function identity gives $\mathbf k\cdot\nabla_{\mathbf k}\omega=0$. \b[Energy propagates along constant-phase surfaces, rather than normal to them.] When $m=0$, $\omega=N$ and $\mathbf c_g=0$; the scalar-product result remains true, but there is no nonzero group-velocity direction to describe geometrically.