This fact is immediate if $M = \mathbb R^n$ since we can write $\partial_i f = \lim \frac{f(x+te_i) - f(x)}{t}$. But I didn't know whether it holds for complete Riemannian manifolds.
Recently, I'm reading Geometric Analysis by Peter Li. The proof of the Cheeger-Gromoll splitting theorem in his book(Theorem 4.4) is more simple than the one on Peterson's book, which avoid the application of smooth support functions, but I don't think it is rigorous enough, so I want some help.
In the last few lines on Page 35, he said that ''By regularity theory, $\beta^+$ is a smooth harmonic function with $|\nabla \beta^+| \leq 1$'', where $\beta^+$ denotes the Buseman function. Obviously, $\beta^+$ is Lipschitz, but I'm wondering how to derive the property that $|\nabla \beta^+| \leq 1$.
For your reference, the book can be found at https://www.cambridge.org/core/books/geometric-analysis/D0A2375D56122B91A0BA370530978248
After viewing the comments, it is true and can be proved as follows.
We choose the normal coordinates centered at $x$, then $|grad \beta(x)| = \sum_j (\partial_j \beta)^2$. Suppose by contradiction that $|grad \beta(x)| > 1$, then there exists a small neighborhood $U$ near $x$ such that $\sum_j (\partial_i \beta)^2 > 1$ for all $y\in U$. Choose $y\in U$ such that the unique minimizing normal geodesic $\gamma$ with $\gamma(0) = x$ and $\gamma(1) = y$ satisfies $\dot\gamma_k(t) = \partial_k\beta(\gamma(t))$ for all $t\in [0,1]$, $k = 1,\ldots, m$. Then by the Lipschitz condition, $$ 1\geq \frac{\beta(\gamma(t)) - \beta(x)}{t} = \frac1{t}\int_0^t \langle{\nabla \beta, \dot\gamma\rangle}\, ds = \frac1{t}\int_0^t \sum_{k=1}^m \partial_k\beta \dot\gamma_k\, ds = \frac1{t}\int_0^t \sum_k (\partial_k \beta)^2 > 1, $$ which is a contradiction.
The reverse implication is also true, as I'll show bellow. I'll present another approach (which is essentially the same as the ones suggested in the comments) based on some other posts I've read on the site. Fix $p \in M$ arbitrarily and let $\gamma: I \to M$ be a geodesic starting at $p$ with initial tangent vector $v$ of unit length (which, from basic Riemannian geometry, is given by $t \mapsto \exp_p(tv)$ for small enough $t$). Suppose that $f$ is $K$-Lipschitz (i.e, $|f(p) - f(q)| \leq K \mathrm{dist}(p, q)$ for all $p, q$ in $M$). Then we have:
$$ \begin{aligned} \langle (\nabla f)_{p}, v \rangle &= \left.\frac{\mathrm{d}}{\mathrm{d} t} (f \circ \gamma) \right|_{t = 0} \\ &= \lim_{t \to 0} \frac{f(\gamma(t)) - f(p)}{t} \\ & \leq \lim_{t \to 0} K \frac{\mathrm{dist}(\exp_p(tv), p)}{t } = K \end{aligned}$$ where in the last inequality we used the hypothesis that $f$ is $K$-Lipschitz and also the basic fact that $\mathrm{dist}(\exp_p(tv), p) = t$ for small enough $t$. Taking $v = \frac{\| (\nabla f)_p \|}{\| \nabla f \|_p}$ we see that $\| (\nabla f)_p \| \leq K$, and since $p$ was taken arbitrarily, it follows that $\| \nabla f \| \leq K$.
Conversely, suppose now that $\| \nabla f \| \leq K$. Fix $\varepsilon > 0$ and $p, q \in M$ arbitrarily. By the definition of the Riemannian distance and the definition of infimum, there exists a smooth curve $\gamma: [0, 1] \to M$ such that $\gamma(0) = p$, $\gamma(1) = q$ and $$\ell(\gamma) = \int_{0}^1 \| \gamma'(t) \| \ \mathrm{d} t \leq \mathrm{dist}(p, q) + \varepsilon$$
Notice also that $f \circ \gamma$ is a smooth curve (with image in $\mathbb{R}$) which joins $f(p)$ and $f(q)$. Again, by the definition of distance and infimum, we have
$$ \begin{aligned} |f(p) - f(q)| &\leq \int_{0}^{1} \vert (f \circ \gamma)'(t) \vert \ \mathrm{d} t \\ &= \int_{0}^{1} \vert \langle (\nabla f)_{\gamma(t)}, \gamma'(t) \rangle \vert \ \mathrm{d} t \leq K \cdot \ell(\gamma) \leq K \cdot( \mathrm{dist}(p, q) + \varepsilon) \end{aligned}$$
where in the first inequality from the second line we used Cauchy-Schwarz. Therefore
$$|f(p) - f(q) | \leq K \cdot \inf_{\varepsilon > 0} (\mathrm{dist}(p, q) + \varepsilon) = K \cdot \mathrm{dist}(p, q) $$
i.e, $f$ is $K$-Lipschitz.