This answer on physics stack exchange says
For a compact domain without boundary (such as the surface of a sphere), you don't need any boundary conditions: there are no non-constant harmonic functions on such domains. (Slick proof of this: you can prove that harmonic functions never have local maxima or minima, but a nonconstant function on such a domain must have them -- in particular, it must have a global maximum and a global minimum somewhere.)
A harmonic function obeys $$\nabla^2 f = \frac{\partial^2}{\partial x^2}f + \frac{\partial^2}{\partial y^2}f = 0$$ And indeed, if $f$ had a local extremum the sign of $\frac{\partial^2}{\partial x^2}f$ and of $\frac{\partial^2}{\partial y^2}f$ would have to be the same, so the sum cannot be $0$, so $f$ cannot have a local extremum, so on a compact domain without boundary it must be constant.
I have two questions about this line of reasoning:
What about local extrema of the form $x^4 y^4$? The second derivative in both directions could be $0$, so I cannot exclude the possibility that a harmonic function could have an extremum like this.
Suppose we are not considering harmonic functions, but rather solutions of $$\nabla \cdot \left[\begin{matrix} m_{11} & m_{12} \\ m_{21} & m_{22}\end{matrix}\right]\nabla f = 0 \label{a}\tag{1}$$ where the $m_{ij}$ can depend on $x,y$. Under what conditions is it still true that there are no nonconstant solutions of $(\ref{a})$ on compact domains without boundary? (presumably the $m_{ij}$ matrix needs to be invertible. In the case of interest to me, it is Hermitian. I could try to change coordinates to a coordinate system whose axes are the eigenvectors of the the $m_{ij}$ matrix, but since those eigenvectors are not guaranteed to be real that would involve introducing complex coordinates...).
Why can't a nonconstant harmonic function have a local minimum? The averaging property. The value of a plane harmonic function at a point $x$ is equal to its average value on any disk centered at $x$.
So, then, if $f(x)\le f(y)$ for all $y$ with $\|y-x\|_2 \le r$, then $$0=f(x)-f(x)=\frac1{\pi r^2}\iint_{\|y-x\|_2\le r}f(y)-f(x)\,dy \ge 0$$ with equality only if the (continuous) nonnegative function we're integrating is identically zero in the disk $\|y-x\|_2\le r$.
Now, I should note that it's nontrivial to define the Laplacian operator $\Delta = \nabla^2$ on the sphere or something like it; the partial derivative operators we built it from are based on a specific set of coordinates, and there's no single set of coordinates that defines the sphere. We can't fall back on the (coordinate-independent) exterior derivative, either; while the gradient map from $0$-forms to $1$-forms and the divergence map from $(n-1)$-forms to $n$-forms are both cases of the exterior derivative, the map we need to bridge the gap from $1$-forms to $n$-forms is strongly coordinate-dependent. Even in the $n=2$ case.
How do we fix this and get something we can work with? To get an idea of what's going on, note that the Hessian matrix of second derivatives transforms by a congruence relation: if $g(x)=f(Ax)$, then $H_g=A^TH_fA^T$. The Laplacian is the trace of this matrix. The Laplacian is the trace of this matrix, so we would like to restrict ourselves to transformations that preserve the trace. What does that? Orthogonal matrices; if $A^T=A^{-1}$, then that's a similarity relation, and the trace is fixed.
Orthogonal matrices preserve distances and inner products. If we want this generalized Laplacian to work out, we need a notion of inner products on our manifold. That's a Riemannian metric - a system of inner products on the tangent spaces. So, then, to find the Laplacian, choose local coordinates $x_1,x_2,\dots,x_n$ so that the inner product at the point we're interested in is the Euclidean inner product. The Laplacian at that point is $\frac{\partial^2}{\partial x_1^2}+\frac{\partial^2}{\partial x_2^2}+\cdots+\frac{\partial^2}{\partial x_n^2}$, and this is invariant under any isometry.
Now, looking back at our argument - the averaging property depends on working in a flat space. If the metric varies point-to-point in those local coordinates we found, so will the formula for the Laplacian - and the averaging property is no longer something we can explicitly write down. What do we get instead? The Laplacian becomes $$\Delta f(x,y) = a(x,y)\frac{\partial^2 f}{\partial x^2}+2b(x,y)\frac{\partial f}{\partial x\partial y}+c(x,y)\frac{\partial^2 f}{\partial y^2} = \text{tr}(MH)$$ in our local coordinates, where the matrix $M=\begin{bmatrix}a&b\\b&c\end{bmatrix}$ is always positive definite and $H$ is the Hessian matrix of second derivatives.
And there's your second part. We need it in order to deal with harmonic functions on any manifold that isn't flat. The extra condition we need falls right out, too - that matrix has to be positive definite everywhere. If $M$ is indefinite, we can move to a basis of eigenvectors and then set $H$ to a positive diagonal matrix while still having that operator zero, giving us a local minimum.
As for the difference between real symmetric and Hermitian? We're adding the off-diagonal elements of $M$ with their transposes anyway, so the antisymmetric imaginary part zeros out. There's no reason to go to complex coordinates unless we started with them.
This more general result follows from the Hopf maximum principle. On a closed disk in our local coordinate patch, the maximum and minimum must occur on the boundary of the disk, so the center can't be a local extremum unless the function is constant. Every point on the manifold has coordinates of this form, so there aren't any local extrema unless the function is constant.