Solution (source code)

= Solution

Let $b=-g\rho'/\rho_0$ be the buoyancy perturbation and $p$ the pressure perturbation divided by $\rho_0$. For stable stratification, $N>0$. The nonrotating <Linearized Boussinesq equations> are
$$
u_t=-p_x,\qquad w_t=-p_z+b,\qquad b_t+N^2w=0,\qquad u_x+w_z=0.
$$
Eliminating $u,p,b$ gives $(\partial_x^2+\partial_z^2)w_{tt}+N^2w_{xx}=0$. Since $w=\zeta_t$, a nonzero-frequency plane wave obeys the same equation for its displacement. Substituting its phase yields the <dispersion relation> for a <plane internal gravity wave>:
$$
\boxed{\Omega^2=\frac{N^2k^2}{k^2+m^2}.}
$$
Thus the frequency depends on the <wavevector> direction rather than its magnitude.

Advection of the background density gives $\rho'=-\zeta\,\bar\rho_z$ to first order. With constant $\bar\rho_z<0$, the instantaneous density gradient is $\rho_z=\bar\rho_z(1-\zeta_z)$. A region has <unstable density stratification> when this becomes positive, namely when $\zeta_z>1$. The maximum of $\zeta_z$ is $|mA_\zeta|$, so the <monochromatic internal-wave overturning criterion> is
$$
\boxed{|mA_\zeta|>1.}
$$
Equality gives a locally vanishing gradient. This is the prediction of the displacement field extrapolated to overturning; the small-amplitude approximation itself ceases to be reliable there.

For the rising packet, distinguish its conserved <absolute frequency> $\omega$ from its actual <intrinsic frequency> $\sigma=\omega-kU(z)$. The printed terminology calls $\omega$ intrinsic while also assigning it to a stationary observer; the stationary-observer interpretation is the one consistent with the displayed Doppler shift. On the positive-frequency branch, the ray Hamiltonian is
$$
\omega(z,k,m)=kU(z)+\sigma(k,m),\qquad \sigma=\frac{Nk}{\sqrt{k^2+m^2}}.
$$
The <Hamiltonian ray-tracing equations> give
$$
\dot x=\partial_k\omega,\quad \dot z=\partial_m\omega,\quad
\dot k=-\partial_x\omega=0,\quad \dot m=-\partial_z\omega=-ks,\quad
\frac{d\omega}{dt}=\partial_t\omega=0.
$$
The last identity follows also by differentiating the Hamiltonian along its canonical trajectory: the spatial and <wavevector> terms cancel in pairs. Thus <absolute-frequency conservation in steady shear> gives constant $k$, constant $\omega$, and constant stationary-observer horizontal phase speed $c_x=\omega/k$. In contrast, $\sigma=\omega-ksz$ decreases as the packet rises. At its initial height,
$$
\omega=\frac{Nk}{\sqrt{k^2+m_0^2}},\qquad
\boxed{z_c=\frac{\omega}{ks}=\frac{N}{s\sqrt{k^2+m_0^2}}.}
$$
This is the <critical level of an internal gravity wave>. In fact $m(t)=m_0-kst$ and $z(t)=[\omega-\sigma(k,m(t))]/(ks)$, so the inviscid ray approaches $z_c$ as $t\to\infty$, rather than reaching it at a finite time.

Write $\theta=|\Theta|=\arctan(|m|/k)$, so $\sigma=N\cos\theta$ with $0<\theta<\pi/2$. The intrinsic <internal-wave phase and group velocity> calculation gives
$$
c_{gx}=\frac{Nm^2}{(k^2+m^2)^{3/2}}=\frac Nk\sin^2\theta\cos\theta,\qquad
c_{gz}=-\frac{Nkm}{(k^2+m^2)^{3/2}}=\frac Nk\sin\theta\cos^2\theta>0.
$$
The observer-frame horizontal ray velocity is $U+c_{gx}$. Dividing it by $c_{gz}$ proves the <internal-wave ray in uniform vertical shear>:
$$
\boxed{\frac{dx}{dz}=\tan\theta+\frac{ksz}{N\sin\theta\cos^2\theta},\qquad
\tan^2\theta=\frac{N^2}{(\omega-ksz)^2}-1.}
$$
The angle increases toward $\pi/2$ and the vertical group speed tends to zero near the <critical level>.

The <wave-action conservation law> fixes the prescribed upward flux. For a nonzero packet, $B>0$, and the given flux relation implies
$$
A_\zeta^2=\frac{B\sigma}{N^2c_{gz}}=
\frac{Bk\sigma}{N^3\sin\theta\cos^2\theta}.
$$
Apply the <monochromatic internal-wave overturning criterion>, using $|m|=k\tan\theta$. After multiplying by the positive trigonometric factors, the exact instability condition is
$$
\boxed{\cot^4\theta<\frac{Bk^3\sigma}{N^3\sin^3\theta}.}
$$
At marginal overturning near a <critical level>, $\sin\theta\simeq1$, so the <wave-action criterion for critical-level overturning> gives
$$
\cot\theta\sim\left(\frac{Bk^3}{N^3}\right)^{1/4}(\omega-ksz)^{1/4}.
$$
With fixed $B,k,N$, this is the requested quarter-power order estimate; the prefactor supplies the dimensions suppressed in that notation. It is an onset balance, not a replacement for $\cos\theta=\sigma/N$. Combining the two relations instead gives $\cos^3\theta=(Bk^3/N^2)\sin\theta$ at onset. Since $m^2A_\zeta^2$ diverges as $\sigma^{-3}$ toward $z_c$, any nonzero packet flux eventually violates the linear overturning criterion before reaching that level, within this nondissipative ray model.