I would like to find an example of a manifold $(M,g)$ along with a non-constant harmonic function having bounded gradient and bounded Hessian.
What came to my mind was that it is known (by the work of Anderson) that for (complete, simply connected) negatively curved manifolds, there are lots of bounded harmonic functions. The gradient estimates then implies bounded harmonic functions have bounded gradient.
But I am not sure about the Hessian. I would be surprised if every bounded non-constant harmonic function on every negatively curved manifold has unbounded Hessian. On the other hand I do not have an argument for the bound. The Bochner formula might be useful, but I also do not know how to estimate $\Delta |\nabla u|^2$ for harmonic function $u$.
For the construction of a concrete example, I have looked at the unit disk with hyperbolic metric, since being harmonic is conformally invariant, bounded harmonic functions are given by Poisson kernel and any boundary data. However the computation is a mess, there are terms in the integral that blows up on the boundary which I guess can be cancelled by the integral but again I do not have a smart choice of boundary function to make the computation doable.
Ideas, comments, hints are all welcomed, thank you!
Edit: The manifold in question should be open (complete and non-compact without boundary). For closed manifolds this is obvious.