The character $\chi_R : G \to \mathbb{C}$ of a representation $R$ is defined by $\chi_R(U) = Tr_R(U)$, namely by taking the trace of $U$ in the representation $R$. See for example Appendix A of Aharony et al.. Then the equation you write seems reasonable: presumably there are $N_f$ hypermultiplets, each with fundamental and anti-fundamental fields that give $Tr(U^m)$ and $Tr(U^{-m})$
respectively in the sum over $R$. The trace is taken of the matrix $U$ because the fields are in the fundamental.
The factor of 2 should come from the details of which representations exactly are being summed over, the field content of the hypermultiplets, etc., which I didn't try to follow.
Now, regarding the distribution of eigenvalues that appears in (C.5), let me give you a physical explanation of where this comes from. They consider the thermodynamics of a theory on a sphere, so the topology is that of a sphere times $S^1$ -- the Euclidean compact time. In this case there is a zero-mode of the gauge field, which comes from the fact that you can turn on a constant field in the thermal direction. This mode cannot be removed by a gauge transformation (in general), so in a path integral you need to integrate over it. This is explained for example in section 4.1 of Aharony et al. You can use a gauge transformation to diagonalize this zero-mode matrix, so you are left with the discrete eigenvalues of the matrix to integrate over. You can also show that these eigenvalues live on a circle, because you can shift them by $2\pi$ (in some normalization) using a gauge transformation, so you need to integrate them over a circle. In the large $N$ limit you have an infinite number of such eigenvalues, but they are still restricted to live on a circle. So you need to describe them using a density, and that is what they call the distribution $\rho$.
Edit: adding explanation of (C.6)
To derive (C.6), consider first the sum
$$ \sum_i \cos(n \alpha_i) = N \int_0^{2\pi} \rho(\theta) \cos(n\theta) $$
where $\rho(\theta) = \frac{1}{N} \sum_i \delta(\theta - \alpha_i)$. Notice that $\rho$ has the correct normalization.
When taking $N \to \infty$ the delta functions in $\rho$ will become very dense, and we will be able to approximate $\rho$ by a smooth function, whose value at $\theta$ depends on the density of the delta functions across the small interval $[\theta,\theta+\epsilon]$. This will be a good approximation because the function $\cos(n\theta)$ that we are integrating over will not vary much in these intervals, so averaging it instead of sampling it will not change the result by much. When we reach the limit $N=\infty$ this will no longer be an approximation, because the delta functions will become continuous.
So this explains the derivation of the second term in (C.6). As for the first term, the idea is similar except that you have two integrals over $\theta,\theta'$ corresponding to the $i,j$ sums. The contributions from the $i=j$ terms ($\theta=\theta'$ in the integral) are subleading in the large $N$ limit: they scale as $N$ while the rest of it scales as $N^2$, so they are neglected.
Now, replace $\cos(n(\theta-\theta')) = \cos(n\theta) \cos(n\theta') + \sin(n\theta) \sin(n \theta')$. Notice that in the path integral over $\rho$ you may consider only eigenvalue distributions that are symmetric under $\theta \to -\theta$, because the original integrand is invariant under this. This is explicitly mentioned for example in Schnitzer, who does the same computation. So this means that $\int d\theta \rho \sin(n\theta) = 0$. I think it should be clear from this point.
