I have no statistics/probability background, so excuse me if this is a dumb question. Consider the following function: \begin{align*} f_N(\vec{x}) = \frac{1}{\sqrt{N}} \sum_{n=1}^N \cos(2\pi \vec{k}_n\cdot \vec{x} + \varphi_n) \end{align*} where $\varphi_n$ is uniformly distributed in $[0,2 \pi]$ and $\vec{k}_n$ is uniformly distributed on the 2-sphere $S_2$. Is there any way that I can make a statement about how that function will look in the limit of $N \to \infty$?
Edit: I encountered this function in my research while trying to simulate the surface of a bicontinuous structure. In something called Cahn's scheme, level sets of this function represent the surface of a randomly assembled bicontinuous structure (e.g. a sponge). The scheme does not mention too much about N's size, so I was wondering how that would look in the limit of large N.
In this answer, I will replace $\varphi_n$ by $\theta_n$ to emphasize that this is a random angle and save the symbol $\varphi$ for other use.
Let $\varphi$ be a compactly supported distribution on $\mathbb{R}^3$ and write $\hat{\varphi}(k) = \int_{\mathbb{R}^3} \varphi(x) e^{-2\pi i k\cdot x} \, dx$ for its Fourier transform. Then
\begin{align*} \mathbb{E}\exp\{ i \langle f_N, \varphi \rangle \} &= \left( \mathbb{E}\exp\left\{ \frac{i}{\sqrt{N}} \int_{\mathbb{R}^3} \varphi(x) \cos(2\pi k_1 \cdot x+\theta_1) \, dx \right\} \right)^N \\ &= \bigg( 1 + \frac{i}{\sqrt{N}} \mathbb{E}\int_{\mathbb{R}^3} \varphi(x) \cos(2\pi k_1 \cdot x+\theta_1) \, dx \\ &\hspace{4em} - \frac{1}{2N} \mathbb{E}\left( \int_{\mathbb{R}^3} \varphi(x) \cos(2\pi k_1 \cdot x+\theta_1) \, dx \right)^2 + \mathcal{O}(N^{-3/2}) \bigg)^N \\ &= \left( 1 - \frac{1}{4N} \mathbb{E} \left[ \left| \hat{\varphi}(k_1) \right|^2 \right] + \mathcal{O}(N^{-3/2}) \right)^N \\ &\xrightarrow[N\to\infty]{} \exp\left\{ -\frac{1}{4}\mathbb{E} \left[ \left| \hat{\varphi}(k_1) \right|^2 \right] \right\} \\ &\hspace{3em} = \exp\left\{ -\frac{1}{4} \int_{\mathbb{R}^3\times\mathbb{R}^3} \varphi(x)\varphi(y) \frac{\sin(2\pi|x-y|)}{2\pi|x-y|} \, dxdy \right\}. \end{align*}
This strongly suggests that $f_N$ will converge in distribution to a Gaussian field with the covariance kernel
$$ C(x, y) = \frac{\operatorname{sinc}(2\pi|x-y|)}{2} = \frac{\sin(2\pi|x-y|)}{4\pi|x-y|}. $$
(Since Lévy's continuity theorem need not hold in general Banach space and I am not an expert of this topic, I will leave this bold claim open.)
Since we may choose $\varphi$ to be a linear combination of finite masses, we at least know that any finite-dimensional distribution of $f_N$ converges in distribution to the normal distribution with covariance matrix $C$, i.e.,
$$ (f_N(x_1), \cdots, f_N(x_d)) \quad \Rightarrow \quad \mathcal{N}(\mathbf{0}, [C(x_i, x_j)]_{i,j=1}^{d}) $$
Here is a density polt of a simulation of the limit using points $(\frac{i}{10},\frac{j}{10},0)$ for $0 \leq i, j \leq 30$.
$\hspace{8em}$![simulation using mesh size $1/10$ on $[0, 3]^2$](https://i.stack.imgur.com/yQjWD.png)