Solution (source code)

= Solution

Expand the square root and quotient in the <logarithmic-height plume equations>:
$$
\widehat M^{1/2}=1+\tfrac12\epsilon M'+O(\epsilon^2),\qquad\frac{\widehat F\widehat Q}{\widehat M}=1+\epsilon(F'+Q'-M')+O(\epsilon^2).
$$
The <linear stability analysis> gives
$$
\boxed{\frac d{d\xi}\begin{pmatrix}Q'\\M'\\F'\end{pmatrix}=J\begin{pmatrix}Q'\\M'\\F'\end{pmatrix},\qquad J=\begin{pmatrix}-q&q/2&0\\m&-2m&m\\f&0&-f\end{pmatrix}.}
$$
In particular, the second row contains $-2mM'$, since the reciprocal of $\widehat M$ contributes another $-M'$.

For the <growing mode of power-law plume similarity>, the <characteristic polynomial> is
$$
P(\lambda)=\lambda^3+(q+2m+f)\lambda^2+\left(\frac32qm+qf+2mf\right)\lambda+qmf.
$$
Within the physical interval $q,m>0$ and $f<0$, $P(0)<0$ while $P(\lambda)\to+\infty$ as $\lambda\to+\infty$. There is a positive real <eigenvalue>, so the similarity equilibrium is \b[unstable with increasing height]. Perturbations behave as $e^{\lambda\xi}\propto z^\lambda$; the result is not stability in physical time. Thus this exact positive-buoyancy similarity trajectory is not a generic attracting far field.

At the upper endpoint,
$$
\boxed{\beta\to-\frac83:\qquad f\to0^-,\quad q\to\frac53,\quad m\to\frac43.}
$$
The limiting polynomial is $\lambda(\lambda+1)(\lambda+10/3)$, reproducing the three stated <eigenvalues>. The positive <eigenvalue> tends to zero, with $\lambda\sim-2f/3$. The two negative modes decay as $z^{-1}$ and $z^{-10/3}$; the remaining mode is neutral at the limiting linear system. Its <eigenvector> is $(1,2,3)$ in relative-flux coordinates.

The <neutral buoyancy-flux mode of a pure plume> has a simple interpretation: a constant surviving <buoyancy flux> fixes the <amplitudes> of a <pure plume>, with $Q\propto F^{1/3}z^{5/3}$ and $M\propto F^{2/3}z^{4/3}$. A small fractional change in $F$ produces changes in $Q,M$ in the ratio $1:2:3$. It is retained rather than damped out.

There is a qualification at the exact endpoint. For fixed nonzero $K$ and $\beta=-8/3$, taking $Q\propto z^{5/3}$ gives $F_z\propto-z^{-1}$, not zero. The flux coefficients above also diverge as $f\to0$. Thus the constant-flux powers are the formal limiting scaling, not an exact constant-buoyancy solution in a nonzero endpoint stratification; cumulative ambient effects can require logarithmic corrections.