Recall that the Ornstein–Uhlenbeck (OU) process in $\mathbb{R}^d$ is defined by the following SDE, $$ d Z_t=\frac{-1}{2} Z_t d t+d W_t, \quad t \geqslant 0 $$ where $\left(W_t\right)_{t \geqslant 0}$ is a standard $d$ dimensional Brownian motion. As it is well known, the process $Z$ admits the following explicit representation $$ Z_t=Z_0 e^{-t / 2}+e^{-t / 2} \int_0^t e^{s / 2} d W_s, \quad t \geqslant 0. $$ Furthermore, the infinitesimal generator of the above process is the following operator $$L:=\Delta -x\cdot \nabla.$$

I was wondering what would be a generalization of this operator to Riemannian Manifolds? I know that Brownian motion can be constructed on manifolds and its generator is simply the Laplace Beltrami operator, however, the OU process has a drift and so I am not sure what is the appropriate way to generalize this to a sphere, for instance.