If you only look inside $[0,2\pi]$, then it is not true that the distribution is Gaussian.

EDIT

Set
$$
S(x;t) = \sum_{p\le x} p^{it},
$$
so that the sum you are looking at is the real part of $S(x;t)$. Then the Prime Number Theorem implies that for $t\in[0,2\pi]$ we have that
$$
S(x;t) = \frac{x^{1+it}}{(1+it)\log x} + O\left(\frac{x}{\log^2x}\right).
$$
Therefore
$$
\mu:=\frac{1}{2\pi} \int_0^{2\pi} S(x;t)dt \ll \frac{x}{\log^2x},
$$
by integration by parts, and
$$
\sigma^2:= \frac{1}{2\pi} \int_0^{2\pi}|S(x;t)|^2 dt \sim \frac{cx^2}{\log^2x},
$$
for some $c>0$. Consequently,
$$
\begin{align}
M_k&:=\frac{1}{2\pi} \int_0^{2\pi}\left| \frac{S(x;t)-\mu}{\sigma}\right|^{2k} dt \\
& = \frac{1}{2\pi} \int_0^{2\pi}\left| \frac{S(x;t)(1+O(1/\log x))}{\sigma}\right|^{2k} dt \\
&\sim \frac{1}{2\pi c^{k}}\int_0^{2\pi} \frac{dt}{(1+t^2)^k}
\in[c_1^k,c_2^k]
\end{align}
$$
for all fixed $k\in\mathbb{N}$, where $c_1$ and $c_2$ are certain positive constants. If the distribution of $S(x;t)$ were Gaussian, then $M_k$ would grow like $c_0^k k!$ for some $c_0>0$, much faster than its current growth.

END OF EDIT

You do get a Gaussian distribution, but you have to look at appropriate ranges of $t$ with respect to $x$ (and we have to look at long intervals of $t$, because we don't understand well what happens inside very short intervals). This is part of the context of the Selberg Central Limit Theorem for the Riemann zeta function. In the proof of his theorem, Selberg shows that the statistical behavior of $\log\zeta(1/2+it)$ for $t\in[T,2T]$ can be modeled by the statistical behavior of $\sum_{p\le T^\epsilon} 1/p^{1/2+it}$ for some $\epsilon=\epsilon(T)$ that tends to 0 slowly. Then he proceeds to estimate moments of the latter sum (which is essentially equivalent to studying the distribution of the sum you're interested in).

The rough idea is that
$$
\begin{align}
\int_T^{2T} S(x;t)^k \overline{S(x;t)}^{\ell} dt
&= \sum_{p_1,\dots,p_{k+\ell}\le x} \int_T^{2T} \left(\frac{p_1\cdots p_k}{p_{k+1}\cdots p_{k+\ell}}\right)^{it} dt \\
&= T\sum_{\substack{ p_1,\dots,p_{k+\ell}\le x \\ p_1\cdots p_k=p_{k+1}\cdots p_{k+\ell}}} 1
+ O\left( \sum_{\substack{ p_1,\dots,p_{k+\ell}\le x \\ p_1\cdots p_k\neq p_{k+1}\cdots p_{k+\ell}}} \frac{1}{|\log(p_1\cdots p_k/(p_{k+1}\cdots p_{k+\ell}))|}\right).
\end{align}
$$
The main term above is 0 unless $k=\ell$, in which case it equals asymptotically $T$ times the $2k$-th moment of a Gaussian. The error term can be shown to be small if $x^{\max\{k,\ell\}}\le T^{1/4}$. Indeed, in this case $p_1\cdots p_k,p_{k+1}\cdots p_{k+\ell}\le T^{1/4}$ and thus $|\log(p_1\cdots p_k/(p_{k+1}\cdots p_{k+\ell}))|\ge 1/T^{1/4}$. So the total error term is $T^{1/4} \pi(x)^{k+\ell} \le T^{3/4}$, which is small enough for our purposes. Hence, if $x=T^{o(1)}$, then an increasingly (as $T\to\infty$) large number of moments is shown to match the moments of the Gaussian, and standard probability results then imply that the distribution of $S(x;t)$ for $t\in[T,2T]$ is Gaussian.

assumethey are iid's, so this is not a complete answer. My question is how to actually go about proving it. Nevertheless your comment is useful, since I was not aware of the paper, but only of later papers by some of the same authors. $\endgroup$ – André LeClair Nov 7 '14 at 0:05