As a part of my bachelor thesis on cosmic structure formation I have been dealing with a sum of many randomly distributed phase factors, so in principle with a pearson random walk. If there are $N$ particles with random locations $\vec{q_j}$ within a sphere of Radius $R$ the sum I have been talking about would be $$s = \sum_{j=1}^N \mathrm{e}^{i \vec{k}\cdot \vec{q_j}}$$
$\phi_s$ denotes the phase of $s$ and it is a quantity which fluctuates between $-\pi$ and $\pi$ on a scale of $1/R$ (can I state this more precisely?). Essentially I want to evaluate the asymptotic proportionalities for $R\rightarrow \infty$ of the integral
$$\int_{\Omega} \mathrm{d}^3k f(\vec{k}) \phi_s(\vec{k})$$
with $\Omega \subset \mathbb{R}^3$. The result should be some distribution with mean zero and a standard deviation which depends on $R$.
What I did so far:
If the problem was in a onedimensional setting with $f = 1$ I would start like this. In order to integrate $\phi_s$ over a finite interval I would "discretize" it: Instead of integrating $\phi_s$ itself I integrate a step function with step size $\frac{1}{R}$ and a random value between $-\pi$ and $\pi$. Now the integral has in principle become a product of the step size $\frac{1}{R}$ and the result of a 1D-random walk: $$\int_a^b \phi_s(k)\mathrm{d}k \approx \sum_{j=1}^{\left \lceil{(b-a)R}\right \rceil} (\pm \pi) \frac{1}{R}$$ with $\pm \pi (j)$ yielding $\pi$ and $-\pi$ with equal probability (actually it should be some value between $-\pi$ and $\pi$ but this should only be wrong by a constant prefactor). For $R\rightarrow \infty$ the average value of this sum is 0 while its standard deviation should be $\frac{\pi}{R}\sqrt{(b-a)R} = \pi\sqrt{\frac{(b-a)}{R}}$.
Now this method can be generalized in various directions. In order to integrate the product of $\phi_s$ with some function $f$ over a possibly unbounded domain in one dimension one can use the same line of thought but now one has to deal with a random walk with varying step size according to $f$. $$\int_a^b f(x) \phi_s(x) \mathrm{d}x \approx \sum_{j=1}^{(b-a)R} f(j)\phi_s(j) \frac{1}{R}$$ Using $\phi_s = \pm \pi$ and the additivity of variances one can write the standard deviation of this term as follows: $$\sigma = \frac{\pi}{R}\sqrt{\sum_{j=1}^{(b-a)R} (f(j))^2} = \frac{\pi}{\sqrt{R}}\sqrt{\sum_{j=1}^{(b-a)R} (f(j))^2 \frac{1}{R}} \xrightarrow{R\rightarrow \infty} \frac{\pi}{\sqrt{R}} \sqrt{\int_a^b (f(x))^2 \mathrm{d}x}$$ In the last step I used the definition of the (possibly improper) Riemann integral. The $D$-dimensional case is completely analogous and yields $$\sigma\left(\int_{\Omega} f(\vec{x}) \phi_s(\vec{x}) \mathrm{d}^D x\right) \approx \frac{\pi}{\sqrt{R}^D} \sqrt{\int_{\Omega} (f(\vec{x}))^2 \mathrm{d}^Dx}$$ for the case where $\phi_s$ fluctuates randomly on a scale of $R$ in every direction and $$\sigma\left(\int_{\Omega} f(\vec{x} \phi_s(\vec{x}) \mathrm{d}^D x\right) \approx \frac{\pi}{\Pi_{i=1}^D \sqrt{R_i}} \sqrt{\int_{\Omega} (f(\vec{x}))^2 \mathrm{d}^Dx}$$ for the case where $\phi_s$ fluctuates randomly on a scale of $R_i$ in the direction $i$.