Gaussian distributions can be generalized from Euclidean space to a wide class of Riemannian manifolds. Gaussian distributions on manifolds are harder to make use of in applications since the normalisation factors, which we will refer to as partition functions, are complicated, intractable integrals in general that depend in a highly non-linear way on the mean of the given distribution. Nonetheless, on Riemannian symmetric spaces, the partition functions are independent of the mean and reduce to integrals over finite dimensional vector spaces. These are generally still hard to compute numerically when the dimension (more precisely the rank $N$) of the underlying symmetric space gets large. On the space of positive definite Hermitian matrices, it is possible to compute these integrals exactly using methods from random matrix theory and the so-called Stieltjes-Wigert polynomials. In other cases of interest to applications, such as the space of symmetric positive definite (SPD) matrices or the Siegel domain (related to block-Toeplitz covariance matrices), these methods seem not to work quite as well. Nonetheless, it remains possible to compute leading order terms in a large $N$ limit, which provide increasingly accurate approximations as $N$ grows. This limit is inspired by realizing a given partition function as the partition function of a zero-dimensional quantum field theory or even Chern-Simons theory. From this point of view the large $N$ limit arises naturally and saddle-point methods, Feynman diagrams, and certain universalities that relate different spaces emerge.