How to solve the $C^\alpha$ Poisson equation on closed Riemannian manifolds? To be specific, suppose $M$ is a closed oriented manifold, $g$ is a Riemannian metric of $M$. 
Let $\Delta_g$ be the Laplace-Beltrami operator w.r.t. $g$. 
Prove: Suppose $f\in C^\alpha(M)$ satisfies $\int_M f\, dVol_g=0$, then there exists a function $u\in C^{2,\alpha}(M)$ such that $\Delta_g u=f$ in $M$, and $u$ is unique up to plus a constant, here $0<\alpha<1$.
My attempt is that, firstly use $D(u):=\frac{1}{2}\int_M( |\nabla u|^2+fu)dVol_g$ is a convex functional with a lower bound on $W_0^{1,2}(M)$ to show that there exists a weak solution $u\in W^{1,2}(M)$, next use the $L^2$-regularity theory to show that $u\in W^{2,2}(M)$, but I don't know how to improve the regularity of $u$ further. (Actually, I can use the method to prove that if $f$ is $C^\infty$, then $u$ is also $C^\infty$, but I cannot extend this result to $C^\alpha$ case.)
Another attempt is Schauder estimate. However, in Gilbarg and Trudinger's book they assume that $u\in C^{2,\alpha}(M)$ already to get some interior derivative norm bound of $u$, while I don't know how to establish $u\in C^{2,\alpha}(M)$. They give a continuity method to ensure that, but it seems their discussion works for domains in Euclidean space, not for manifolds. Therefore, I want to split the question into coordinate charts, but I failed, because I don't know how to use the condition $\int_M f\, dVol_g=0$ and how to give boundary conditions in every coordinate charts.
Since I'm a novice in PDE, my presentation of the problem might have some errors. Please correct them by comments or answers. Also, any comments or answers are welcome.
Thanks for your help.
Remark: I've already asked this question on math.stackexchange.com, but nobody replied. Maybe this question is not so suitable for MO, but I really want to get an answer.
 A: Careful. C^\infty is not dense in C^{\alpha}. Nonetheless you can first produce an L^2 solution then use Schauder estimates to show that the solution is C^{2,\alpha}. This last step is localizable -- you can multiply by a smooth cutoff function (and use some preliminary
estimates to say that the solution is C^2, for example) and transfer to a coordinate chart, then apply the estimates and continuity method in Gilbarg&Trudinger. There are other cleaner ways to proceed: the solution u = G f is given by a pseudodifferential operator of order -2 and it is not hard to prove (see Taylor, Vol. III) that such a psido is bounded from
C^\alpha to C^{2,\alpha}.
A: Elliptic regularity in Hölder spaces is precisely what you need here.
If you want to prove it along the lines you described, using $L^2$-theory, then you can first establish that smooth $f$ leads to smooth solution $u$. Then you approximate $f$ in the $L^2$-norm by smooth functions with uniformly bounded $C^{\alpha}$-norms (this is important as smooth functions are not dense in Hölder spaces), and use the Schauder estimate. It suffices to do everything locally.
