In Euclidean vector space, we can parameterize multivariate normal distribution by its mean vector and covariance matrix.
On a general sphere $S^n$, starting from a dirac delta and simulate the heat equation, one can also obtain an gaussian-like distribution, as answered in this question.
Now that I wonder if its possible to define a "spherical multivariate normal distribution". Obviously, the "mean vector" on a sphere would be an unit vector. However, how should one properly define "covariance matrix" on $S^n$?
Take a loot at 5-parameter Fisher–Bingham distribution/Kent distribution.