If Delta f=0 , show f_x,f_y,f_(xx),f_(xy),f_(yy) are also…
If $\Delta f=0$, show $f_x,f_y,f_{xx},f_{xy},f_{yy}$ are also harmonic, and generalise to any linear PDE with constant coefficients.
Explanation of Method
Differentiate the equation $\Delta f=0$ itself and use that partial derivatives commute (Clairaut) and that the Laplacian's coefficients are constants, so it commutes with $\partial/\partial x$ and $\partial/\partial y$.
Step 1 · First derivatives
\[ \Delta(f_x) = \partial_x(f_{xx})+\partial_x(f_{yy}) = \partial_x(\Delta f)=\partial_x(0)=0, \] and identically $\Delta(f_y)=\partial_y(\Delta f)=0$.
Step 2 · Second derivatives
Since $f_x,f_y$ are themselves harmonic, applying the same argument to them gives $\Delta(f_{xx})=\partial_x(\Delta f_x)=0$, $\Delta(f_{xy})=\partial_y(\Delta f_x)=0$, $\Delta(f_{yy})=\partial_y(\Delta f_y)=0$.
Step 3 · General linear operator with constant coefficients
For $L=\sum_{|\alpha|\le m}c_\alpha\partial^\alpha$ with constant $c_\alpha$, and any $D_j=\partial/\partial x_j$: since constants commute with differentiation and mixed partials commute, \[ L(D_jf) = \sum c_\alpha\,\partial^\alpha(D_jf) = \sum c_\alpha\,D_j(\partial^\alpha f) = D_j(Lf) = D_j(0)=0, \] i.e. $L\circ D_j=D_j\circ L$. By induction every derivative $\partial^\beta f$ of any order satisfies $L(\partial^\beta f)=0$.
Why constant coefficients matter
If $c_\alpha=c_\alpha(x)$, the product rule adds an extra term $\sum(\partial c_\alpha/\partial x_j)\partial^\alpha f$ that generally doesn't vanish, breaking $L\circ D_j=D_j\circ L$ — so the result is special to constant-coefficient operators (the Laplacian being one example).
Result
Conclusion: