Limit of series of sine terms with random arguments 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.
 A: 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}$ 
