I'm a physicist, somewhat familiar with stochastic processes, but I'm a little unsure of what follows. What I basically have is a complicated quantity involving a vector that is equivalent to diffusion on a unit sphere. I need to average this quantity, and I am able to do so if I have the characteristic functional of the process. Here are the details:

Suppose we have standard Brownian motion in $d$ dimensions, represented by the usual Langevin equation

$\frac{d}{dt} \mathbf{x}(t) = \boldsymbol{\xi}(t)$,

where $\boldsymbol{\xi}(t)$ is Gaussian noise. Then the characteristic functional (not function) is a Gaussian functional,

$P[\mathbf{k}] = \langle e^{i \int_{0}^{\infty} \mathbf{k}(t') \cdot \mathbf{x}(t') \, dt'} \rangle = e^{i \int_{0}^{\infty} k_{i}(t') \langle x_{i}(t') \rangle \, dt' - \frac{1}{2} \int_{0}^{\infty} \int_{0}^{\infty} k_{i}(t') \langle x_{i}(t') x_{j}(t'') \rangle k_{j}(t'') \, dt' dt''}$.

In the above equation, repeated Cartesian indices are implicitly summed.

I'd like the equivalent of $P[\mathbf{k}]$ for Brownian motion on a unit sphere in $d$ dimensions (or even just $3$ dimensions). From what I've read in the literature, this process has the Langevin equation

$\frac{d}{dt} n_{i}(t) = \sqrt{2 D} (\delta_{i j} - n_{i}(t) n_{j}(t)) \xi_{j}(t) - 2 D n_{i}(t)$,

where $\boldsymbol{\xi}(t)$ is again Gaussian noise which we can take to be uncorrelated.

What is the characteristic functional for this process? Is it Gaussian? If it is not Gaussian and the formula is unobtainable, can I justify approximating it by a Gaussian?

Edit: For full information, the quantity that I need to compute is $\bigg\langle e^{i \mathbf{a} \cdot \int_{0}^{t} \mathbf{n}(t') \,dt'} \frac{d}{d t} \mathbf{n}(t) \bigg\rangle$ for arbitrary $\mathbf{a}$. If there is an easier way to do this, that would be great. I am particularly interested in the long time limit.

Edit 2: Okay, I believe the characteristic functional is definitely not Gaussian. For rotational diffusion, the probability distribution is a sum of spherical harmonics. That should be equivalent if we set the $\mathbf{k}(t) = \mathbf{k} \delta(t)$ and Fourier invert. Anyway, is there an easy way to derive this functional from a higher dimensional Brownian motion or justify the Gaussian approximation?