Distribution of primitive roots, as p varies

Let me adjust your quantity slightly as $$\tilde D_p(f) := \frac{1}{\phi(p-1)}\sum_{x \in \Phi(p)} f \left( \frac{x}{p} \right),$$ and let me initially impose the natural condition $f(0)=f(1)$. Then I claim that $$\lim_{p\to\infty}\tilde D_p(f)=\int_0^1 f(t)\,dt.$$

It suffices to show the claim when $f(t)=e^{2\pi iat}$, where $a$ is a nonzero integer. In this case, by Lemma 2 in Carlitz's paper, $$\tilde D_p(f)=\frac{1}{\phi(p-1)}\sum_{x \in \Phi(p)}e_p(ax)=\frac{1}{p-1}\sum_{d\mid p-1}\frac{\mu(d)}{\phi(d)}\sum_{\substack{\chi\bmod{p}\\\mathrm{ord}(\chi)=d}}\sum_{x=1}^{p-1}\chi(x)e_p(ax).$$ For $p$ exceeding $a$, the inner Gauss sum is always of absolute value at most $\sqrt{p}$, hence $$\left|\tilde D_p(f)\right|\leq\frac{\tau(p-1)}{p-1}\sqrt{p}=p^{-1/2+o(1)},$$ which tends to zero quite optimally. The proof is complete.

Now we can remove the condition $f(0)=f(1)$ by a simple formal argument. Indeed, altering $f$ in a small neighborhood of zero has little effect on $\tilde D_p(f)$ or $\int_0^1 f(t)\,dt$ by what we have proved already (examine bump functions around zero). Thanks to Will Sawin for pointing this out!


A more general result (for general finite fields) is proved in.

Perel’muter, G.I.; Shparlinskij, I.E., The distribution of primitive roots in finite fields, Russ. Math. Surv. 45, No.1, 223-224 (1990); translation from Usp. Mat. Nauk 45, No.1(271), 185-186 (1990). ZBL0705.11078.