pr.probability – Divergence-free Gaussian vector field with given mean magnitude and correlation function

My general question is how one might construct an isotropic random vector field $vec f: mathbb{R}^3 to mathbb{R}^3$ which has a given mean magnitude $mathbb{E}(||vec f(vec x)||)=mu$ such that vector magnitudes and direction are correlated up to some length scales $l$ (beyond which the correlation goes to zero). Furthermore, we wish that $nabla cdot vec f=0$, although a construction which does not satisfy this condition is already very interesting.

More precisely, we are given a vector $vec{mu}$ in $mathbb{R^3}$ and matrix-valued correlation function $C(vec x_1,vec x_2): mathbb{R^3} times mathbb{R^3} to M_3(mathbb{R})$ which is isotropic, i.e $C(vec x_1,vec x_2)=C(||vec x_1-vec x_2||)$. One may define a Gaussian Process $f(vec x) sim GP(vec{mu},C(vec x_1,vec x_2))$ such that $mathbb{E}(vec f(vec x))=vecmu$ and $Cov(vec f(vec x_1),vec f(vec x_1))=C(vec x_1,vec x_2)$. I believe it is known how to generate such a random field $f$. It seems the equivalent problem for a scalar field $f$ is well-studied, where one can for example draw the field first in Fourier space by normalizing a white noise field with the appropriate power spectrum $P(k)$, and then transform back to real space. However it would be useful if someone could describe here a simple procedure for the case of $mathbb{R}^3$, which I think is known but I haven’t found a clear description anywhere.

Question: how can one generate such a random field $f$ if one imposes a zero mean vector $vec{mu}=vec 0$, and in addition mean magntitude $mathbb{E}(||vec f(vec x)||)=mu$? And now if we impose also $nabla cdot vec f=0$?