Imagine first a disk that is mostly stationary, except for random ("thermal" if you like) "rotational fluctuations" around its axis (which is fixed). Something a bit like what's shown in the figure below, except that the rotational fluctuations should be much smaller than what the arrows suggest:

This post is largely about making this "rotational fluctuations" business a bit more precise, and considerably more general.
Since the space of all rotations is 1-dimensional, one could conceivably model these "rotational fluctuations" as any other 1-dimensional random walk/Brownian motion, at least when the "rotational position $\theta$ is in the vicinity of 0"1, 2.
Now, jettison the disk (it was just a visualization aid), and rephrase the problem as one of defining something like a random walk, or like Brownian motion, in the space of all rotations of the plane about the origin, aka the space of all orthogonal transformations of $\mathbb{R}^2$ with determinant 1.
Finally, what I'm really interested in is this: what is the $n$-dimensional generalization of the $n = 2$ example that I've described (very clumsily) above?
Specifically:
what are the analogues of (discrete) random walks and (continuous) Brownian motion in the space of all orthogonal transformations with determinant 1 of $\mathbb{R}^n$ onto itself?3.
what algorithms can I use to sample from the space of rotations in $\mathbb{R}^n$ such that the resulting sequence of "rotational fluctuations" will simulate these rotational analogues of random walks and Brownian motion?
1Or, if one envisions $\theta$ as a function of time $t$, then maybe one can replace "$\theta$ in the vicinity of 0" with "$\theta(t)$ for small $t$".
2One important question arises even at this most elementary stage: if we remove the "near 0", or "small $t$", condition, is it still valid to use 1-d translational models to model these rotational fluctuations, and if not, what should be used instead?
3 I'm expressly avoiding learned lingo like $SO(n)$, since it would encourage answers that would certainly be well over my head. Somehow, I learn the talk much, much, much faster than I learn the walk.
This is a non-rigorous answer.
If $t \mapsto R_t$ is a differentiable function from $\mathbb R$ to $SO(n)$, then its derivative satisfies $$ \frac{dR_t}{dt} = R_t A_t ,$$ where $A_t$ is an anti-symmetric matrix. I could rewrite this non-rigorously as $$ dR_t = R_t A_t \, dt = R_t (\exp(A_t \, dt) - I) ,$$ where $\exp(A)$ denotes the matrix exponential of a matrix $A$.
Now to get a random walk on $SO(n)$, we would prefer to choose the anti-symmetric matrix to be white noise. Thus we could let $$ A_t^{(j,k)} = -A_t^{(k,j)} = dW_t^{(j,k)} \tag 1$$ for $1 \le j < k \le n$, where $W_t^{(j,k)}$ is an independent collection of Wiener processes. Now by the Ito calculus, we have that $dW_t^{(j_1,k_1)} \cdot dW_t^{(j_2,k_2)} = dt$ if $j_1 = j_2$, $k_1 = k_2$, and zero otherwise. Hence multiplying $A$ by itself, we get $$ A_t^2 = -(n-1)I \, dt .$$ Hence using Taylor's series, we see $\exp(A_t) = I + A_t - \frac{n-1}2 I \, dt$, and hence $$ dR_t = R_t ( A_t - \tfrac{n-1}2 I\, dt) \tag 2$$ defines a random walk on $SO(n)$.