Solution (source code)

= Solution

Put $R=h_t-\nu\nabla^2h$, where the <diffusion coefficient> $\nu$ is constant. Then
$$
\delta S=2\int R(\partial_t\delta h-\nu\nabla^2\delta h)\,d^2x\,dt.
$$
One temporal and two spatial <integrations by parts> give
$$
\delta S=-2\int(\partial_tR+\nu\nabla^2R)\delta h\,d^2x\,dt,
$$
with compactly supported variations. Hence
$$
\boxed{(\partial_t+\nu\nabla^2)(\partial_t-\nu\nabla^2)h=0,
\qquad h_{tt}-\nu^2\nabla^4h=0.}
$$
Here $\nabla^4=(\nabla^2)^2$ is the <biharmonic operator>. The cross terms cancel because $\nu$ is constant and spatial and temporal derivatives commute. Equivalently, $P^0=2R$ and $Q^{ij}=-2\nu R\delta^{ij}$ in part (b), producing $-2\partial_tR-2\nu\nabla^2R=0$.

This <squared diffusion-residual action> does not have only the first-order <diffusion equation> as its stationary equation. The zero-residual paths $h_t=\nu\nabla^2h$ have $S=0$ and are certainly solutions, but the variational equation is the higher-order one above. For a spatial <Fourier mode> with nonzero <wavevector> $\mathbf k$, it gives
$$
\ddot h_{\mathbf k}-\nu^2|\mathbf k|^4h_{\mathbf k}=0,\qquad
h_{\mathbf k}(t)=A_{\mathbf k}e^{-\nu|\mathbf k|^2t}+B_{\mathbf k}e^{+\nu|\mathbf k|^2t}.
$$
Boundary data or an additional physical selection principle are needed to retain only the decaying diffusion branch. The spatially uniform mode instead has $h_0(t)=A_0+B_0t$ in the stationary equation.