  This paper presents a novel formula for the transition density of the
Brownian motion on a sphere of any dimension and discusses an algorithm for the
simulation of the increments of the spherical Brownian motion based on this
formula. The formula for the density is derived from an observation that a
suitably transformed radial process (with respect to the geodesic distance) can
be identified as a Wright-Fisher diffusion process. Such processes satisfy a
duality (a kind of symmetry) with a certain coalescent processes and this in
turn yields a spectral representation of the transition density, which can be
used for exact simulation of their increments using the results of Jenkins and
Span\`o (2017). The symmetry then yields the algorithm for the simulation of
the increments of the Brownian motion on a sphere. We analyse the algorithm
numerically and show that it remains stable when the time-step parameter is not
too small.
