An algorithm for simulating Brownian increments on a sphere

TL;DR

The paper tackles the problem of simulating increments of Brownian motion on the sphere $\mathbb{S}^{d-1}(R)$ where a closed-form transition density is unavailable. It introduces an exact density representation obtained by transforming the geodesic-distance process into a Wright–Fisher diffusion with time-change $\tau=2Dt/R^2$ and leveraging moment duality with coalescent processes, yielding a Beta-mixture form for the angular density $\rho(\theta,\tau)$. This leads to an exact, dimension- and scale-agnostic sampling algorithm for spherical Brownian increments, generalizing prior methods to arbitrary $d$, $R$, and $D$, and centered on precomputed coefficients $q^{d-1}_m(\tau)$ and inversion sampling. The work includes a thorough numerical analysis of stability, shows practical applicability for moderate to large time-steps, and discusses fallback approximations for very small $\tau$, positioning the method as a robust complement to tangent-space approximations in manifold diffusion simulations. Overall, the approach provides a principled, exact mechanism for simulating spherical increments with broad applicability in physics, biology, and related disciplines.

Abstract

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ò (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.

Figures (3)

• Figure 1: Absolute errors $\left|\rho^d_1-\rho^d_2\right|$ on the left plot and relative error $\rho^d_1/\rho_2^d-1$ on the right plot for $d=3,4$.
• Figure 2: Plots of values $b^{(\tau,d-1)}_{m+i}(m)$ and of partial sums $\sum_{k=m}^{m+i}(-1)^{k-m}b^{(\tau,d-1)}_k(m)$ for $i=0,\ldots,20$ and parameters for the left plot: $d=3,\tau=0.5,m=2$ and for the right plot: $d=3,\tau=0.1,m=13$.
• Figure 3: Plots of values $q_m^{d-1}(\tau)$ for $m=0,\ldots,40$ and parameters for the left plot: $d=3$ and times $\tau=0.1,0.2,0.3,0.5$ and for the right plot: $\tau=0.1$ and dimensions $d=3,8,15,30.$