We find necessary and sufficient conditions for a finite $K$-bi-invariant measure on a compact Gelfand pair $(G, K)$ to have a square-integrable density. For convolution semigroups, this is equivalent to having a continuous density in positive time. When $(G,K)$ is a compact Riemannian symmetric pair, we study the induced transition density for $G$-invariant Feller processes on the symmetric space $X = G/K$. These are obtained as projections of $K$-bi-invariant Lévy processes on $G$, whose laws form a convolution semigroup. We obtain a Fourier series expansion for the density, in terms of spherical functions, where the spectrum is described by Gangolli's Lévy-Khintchine formula. The density of returns to any given point on $X$ is given by the trace of the transition semigroup, and for subordinated Brownian motion, we can calculate the short time asymptotics of this quantity using recent work of Bañuelos and Baudoin. In the case of the sphere, there is an interesting connection with the Funk-Hecke theorem.