I find it strange that Besse does not discuss this! Here's my understanding of the issue:

Sun and Wang care about analyticity because they want to apply the Łojasiewicz--Simon inequality. This inequality was discovered in the finite dimensional setting by Łojasiewicz. Leon Simon later figured out how to use it in an amazing way to treat uniqueness of asymptotic limits of parabolic and elliptic problems (the original application was to the uniqueness of tangent cones at certain singularities of minimal surfaces and harmonic maps). Lately it has been applied to Ricci flow in a few places, including Sun--Wang.

The Łojasiewicz--Simon inequality was originally proven in L. Simon "Asymptotics for a class of nonlinear evolution equations, with applications to geometric problems" Ann. Math. 1983 http://www.ams.org/mathscinet-getitem?mr=727703.

The key fact is that the Łojasiewicz inequality for finite dimensional analytic functions $f:\mathbb{R}^n\to \mathbb{R}$ says that there if $\nabla f (0)= 0$ and $f(0) = 0$ then is some neighborhood of $0$, say $U\subset \mathbb{R}^n$, some constants $C>0$ and $\theta\in(0,\frac 12)$ so that
$$
|\nabla f(x)| \geq C|f(x)|^{1-\theta}
$$
for $x \in U$. However, this is FALSE for $f$ merely smooth. It's fun to try to find a counterexample.

@Deane Yang, the basic idea behind this is as follows: suppose that $x(t)$ is the (negative) gradient flow of $f$ and $f \geq 0$ on $U$. Then, we claim that $x(t)$ converges to some $x_\infty$ in $U$ as long as $x(0) \in V$ (some smaller neighborhood of $0$ $U$).

To see this, we first claim that the Łojasiewicz inequality will be valid along the entire flow $x(t)$, as long as we start very close to $0$
$$
\frac{d}{dt} f(x(t))^{\theta} = \theta f^{\theta - 1} \nabla f(x(t)) \cdot \dot x(t) = -\theta f^{\theta - 1} |\nabla f(x(t))|^2 \leq -\theta C|\nabla f(x(t)) = -\theta C |\dot x(t)|
$$
Integrate this from $0$ to $T$, so
$$
f(x(T))^\theta - f(x(0))^\theta \geq - \theta C \int_0^T |\dot x(t)|
$$
This shows that (as long as $f(x(0))$ was close enough to $0$) $x(t)$ cannot have too long of a length. Thus, shrinking $U$ to $V$ if necessary, we can control the length so that $x(t)$ cannot leave $U$. Thus, the Łojasiewicz inequality will be valid along all of $x(t)$.

Now, we can do a similar argument to get the convergence:
$$
\frac{d}{dt} f(x(t))^{2\theta -1} = (1-2\theta) f(x(t))^{2\theta -2} |\nabla f(x(t))|^2 \geq C
$$
Integrating this gives
$$
f(x(t)) \leq Ct^{\frac{1}{2\theta -1}}
$$
Now, you can combine the above two inequalities to get a bound on $\int_{T_1}^{T_2} |\dot x(t)|$ which is strong enough to prove convergence to some critical point of $f$.

Simon's argument basically reduces an infinite dimensional "analytic functional" to a finite dimensional piece (i.e. the kernel of the linearization) and then applies Łojasiewicz's result. The key is that this "reduction" is analytic. (In the Sun--Wang paper, there are two additional complications as compared to Simon's paper: (1) There is a gauge group and (2) The $\mu$ functional is an inf over $f$'s)

I think that the simplest definition of analyticity which is sufficient for your purposes is that a map $G:M\to N$ between Banach manifolds is real-analytic if for $x \in M$, there are coordinate neighborhoods $U$ of $x$ and $V$ of $G(x)$ so that $G:U\to V$ is a real analytic map between subsets of a Banach space. What does this mean? For $y\in U$ and $v \in T_yU$ small enough, we require that there are $G_k(y,v)\in V$ (depending in an appropriate way on $y$ and $v$, I'll say more in a minute) so that
$$
G(y+\lambda v) = \sum_{k\geq 0\} G_k(y,v) \lambda^k
$$
as a convergence sum in $V$. Furthermore we require uniform estimates on the $G_k$. See Simon's paper, p 529, for a version of this (He's only considering real analytic maps from a Banach space to $\mathbb{R}$, because in the end thats all we care about: we'd like to apply the Łojasiewicz inequality to some *functional*)

Then, as you say, it is obvious that e.g. scalar curvature is a real analytic map $Met^k(M)\to H^{k-2}(M)$, where $Met^k$ is the Sobolev space of "$k$-times differentiable" metrics (for $k$ large). A proof of this would go something like this: scalar curvature is built from the metric from a few pieces, e.g. covariant differentiating and tracing, and you might check that each of these pieces gives a real analytic map between the respective Banach manifolds, and then you might prove that the composition of real analytic maps is real analytic.

I can't seem to find a good reference.. You are right that it does not seem to be in Besse (although they implicitly use this fact in the section on "Moduli Spaces of Einstein Structures" when they say e.g. that the moduli space is an analytic variety. If you find anything, please let me know here, as it would be nice to know where its written down!