Local phase statistics of the nontrivial Riemann zeros
This is an extended version of my comments above. The phase is the imaginary part of $\log(\zeta^\prime(\rho)))$. For the real part, we have the following theorem due to Hejhal (1):
Assuming the Riemann Hypothesis and a technical condition on the spacing of zeros which is a weak consequence of the Montgomery Pair Correlation Conjecture, then $\log |\zeta^\prime(1/2+i t)|$, suitably normalized, converges in distribution over fixed ranges to a standard normal variable. More precisely, for $\alpha<\beta$ we have
$ \lim_{N\to\infty}\left|\frac1N\left\{n:N\le n\le 2N,\alpha<\frac{\log\left|2\pi\zeta^\prime(1/2+i\gamma_n)/\log(\gamma_n/(2\pi))\right|}{\sqrt{\log\log(N)/2}}<\beta\right \}\right|= \frac{1}{\sqrt{2\pi}}\int_\alpha^\beta\exp(-x^2/2)\, dx$
This is an analog for $\zeta^\prime(s)$ of an (unconditional) result of Selberg for both the real and imaginary parts of $\log \zeta(s)$. Hiary and Odlyzko have investigated this numerically, and find the convergence rather slow. They also observe a surplus of large values and a deficit of small values. See the screenshot of Figure 1 from their paper:
I tried doing my own calculation in Mathematica, using the first 10,000 [update: 100,000] zeros above height $10^6$. The convergence of the mean to $0$ is so slow I first thought I had omitted a logarithm:
I would expect a result similar to Hejhal's to be true for the imaginary part of $\log \zeta^\prime(\rho)$, suitably defined. (I'm interested in knowing if if this has already been done.) I hope also to do a numerical calculation in Mathematica soon; if so I will provide an update. But again the convergence is expected to be slow.
Update: Below is a histogram of $\arg(\zeta'(1/2+i\gamma))-\pi+\gamma\log(2)$, where the argument is computed via continuous variation up the line $\sigma=4$ and then horizontally, for the first 100,000 zeros above height $10^6$. If the asymptotic behavior is indeed Gaussian, the lower order terms are contributing something interesting.
Update July 2020: This paper extends the computation described above to $5\cdot 10^6$ zeros.
(1) D.A. Hejhal, "On the distribution of $\log |\zeta^\prime(1/2+i t)|$", in Number Theory, Trace Formulas, and Discrete Groups, K.E. Aubert, E. Bombieri, D.M. Goldfeld, eds., Proc. 1987 Selberg Symposium, Academic Press, 1989, 343-370.
Inspired by the OP question 'Why are they distributed as such?' I took another look at the paper of Hejhal on the real part of $\log\zeta^\prime$ cited in the previous answer, and (after 30 years) I finally understand a little. Here's my exposition of Hejhal's exposition (p. 346) of the basic idea behind the proof. First some notation:
With $\chi(s)=\pi^{s-1/2}\frac{\Gamma((1-s)/2)}{\Gamma(s/2)}$, define $\phi(s)$ by $\chi(s)^{-1/2}=\exp(i\phi(s))$, so $\phi$ is real on the critical line.
Let $M$ be a large constant, $t$ an auxiliary random variable with $T\le t\le 2T$, and $W_t$ the 'window' $[t-M/\log T,t+M/\log T]$.
Let $A(t)=\frac{t}{2\pi}\log\left(\frac{t}{2\pi}\right)-\frac{t}{2\pi}$.
Let $x=A(u)$, and $\theta(u)=\phi(1/2+iu)$. Let $P_t(x)$ be the polynomial approximation $\Pi_{\gamma\in W_t}(x-A(\gamma))$, and define $\Omega_t(u)$ to be the correction to the approximation so that
$ \zeta(1/2+iu)=\exp(\Omega_t(u)-i\theta(u))P_t(x). $
Computing logarithmic derivative (in $u$, being careful with the chain rule) we see
$ i\zeta^\prime(1/2+iu)/\zeta(1/2+iu)=\Omega_t^\prime(u)-i\theta^\prime(u)+P_t^\prime(x)/P_t(x)\cdot A^\prime(u) $
Rearranging gives
$ i\zeta^\prime(1/2+iu)/A^\prime(u)=\zeta(1/2+iu)\left(\Omega_t^\prime(u)/A^\prime(u)-i\theta^\prime(u)/A^\prime(u)+P_t^\prime(x)/P_t(x)\right). $
Hejhal will make an estimate (see below) of the term in parenthesis on the right to argue that
$\log\left|\zeta(1/2+iu)\right|/\sqrt{\log\log T}$ and $\log\left|\zeta^\prime(1/2+iu)/A^\prime(u)\right|/\sqrt{\log\log T}$ are (in effect) the same random variable, and so Selberg's theorem (mentioned in the previous answer) applies.
For this estimate, Hejhal claims he and Bombieri showed previously that the total variation of $\Omega_t(u)$ on $W_t$ is $O_M(1)$ for 'most' $t$. From this,
$ \int_{W_t}\left|\Omega^\prime_t(u)\right|du=O_M(1) $ for 'most' $t$, and so on 'most' windows $W_t$, $\left|\Omega^\prime_t(u)\right|=O_M(\log T)$.
The above was the hard part; $\theta^\prime(u)/A^\prime(u)$ is elementary. And $P_t^\prime(x)/P_t(x)=\sum_{\gamma\in W_t}1/(x-A(\gamma)),$ with average spacing between $A(\gamma)$ being 1 and the number of terms in the sum $O_M(1)$. Hejhal argues heuristically that
$ \log\left|\Omega_t^\prime(u)/A^\prime(u)-i\theta^\prime(u)/A^\prime(u)+P_t^\prime(x)/P_t(x)\right|=O_M(1) $ except for a subset of small measure.
From this one see the difficulty in extending this to the imaginary part of $\log\zeta^\prime(1/2+iu)$, defined (say) by continuous variation up the vertical line from $4$ to $4+iu$ and along the horizontal line from $4+iu$ to $1/2+iu$. Namely, the estimates depend on $u$ being inside the window $W_t$.