Consider a smooth, connected and complete Riemannian manifold $M$. It is well known that harmonic functions defined on some open subset of $M$ are $C^\infty$.
I'm interested in knowing whether there are quantitative $C^1$ estimates for such harmonic functions which depend only on some bound from below on the curvature (either sectional or Ricci) and an upper bound on the dimension.
To give an example of what I mean, consider the same question with "$C^1$" replaced by "Lipschitz". It is then a classical fact of geometric analysis that Lipschitz estimates follows from a lower bound on the Ricci and an upper bound on the dimension, the proof being an application of the Bochner inequality together with elliptic regularity theory.
Yet, I don't know any result which improves such Lipschitz estimate to $C^1$. I'm interested both in positive statement (e.g. with a lower bound on the sectional there is an estimate) and in negative ones (e.g. a lower bound on the Ricci is not sufficient because along a certain sequence of manifolds the $C^1$ norm of properly defined harmonic functions blows-up).
The question comes from the study of non-smooth spaces arising as (measured)-Gromov-Hausdorff limits of manifolds with prescribed curvature bounds. Basically, what I'm trying to understand is "how smooth" can be functions defined on them, so that it is natural to start looking at harmonic ones. Still, in fact I don't really know how smooth harmonic functions are on Riemannian manifolds, and given that any estimate on limit spaces is (more-or-less) also a quantitative estimate on the smooth ones, before trying to reinvent the wheel I'd like to know if there are known results available in the literature.
Any help is appreciated, thanks.