Let me start with a nice experimental observation. Fix a large prime, say $p = 5003$. It turns out that $g = 2$ is a primitive root mod $p$. If we plot the powers of $g \in \Bbb F_p^{\times}$ (considering representatives between $0$ and $p-1$), we obtain this$^{[1]}$ :
$\hspace{2cm}$
We see that it looks pretty random! So we might wonder how powers of a given primitive root modulo $p$ distribute in $\{1, ..., p-1\}$, as $p$ goes to infinity.
To make it more precise, for every prime $p$, fix a primitive root $g_p$ (i.e. it generates $\Bbb F_p^{\times}$), and define $$\phi_p : \Bbb Z / (p-1) \Bbb Z \to \Bbb F_p^{\times}, k \mapsto g_p^k.$$
let $e_p(x) = \exp(2 \pi i x / p)$ and consider the finite subset $$E_p = \{ ( \,e_p(k) \; , \; e_p(\phi_p(k)) \,) \mid 1 \leq k \leq p-1 \} \subset S^1 \times S^1$$ of the complex torus (we might think simply of the fundamental domain $[0,1[^2$ as well, as in the above picture).
The question is :
When $p \to +\infty$, do the sets $E_p$ equidistribute in $S^1 \times S^1$ with respect to Lebesgue (Haar) measure ?
By Weyl theorem, it can be shown that it boils down to prove that for every $(m, n) \in \Bbb Z^2 \setminus \{ (0,0) \},$ one has $$\Sigma_{n, m}(f ; p) := \sum_{1 \leq x \leq p-1} e_p(mx + n g_p^x) = o(p), \quad p \to \infty$$
I did not find any reference dealing with such exponential sums. Any help would be appreciated.
$^{[1]}$ This is done using the following code in SAGE (https://cocalc.com/ or https://sagecell.sagemath.org/)
powers = [[0, 1]];
N = 2;
p = next_prime(5000);
for j in range(1, p):
powers = powers + [[j, mod(N, p)^j]];
Plot1 = list_plot(powers, size=5, color='blue')
show(Plot1, figsize=[8,8])
save(Plot1, "powers_of_2_modulo_5003.pdf", figsize=[8,8])
Taking a look at a little bit different exponential sums that I find more natural in this context. The exponent $x$ more naturally takes values in the additive group $\Bbb{Z}/(p-1)\Bbb{Z}$. So I want to look at the picture as a subset of a $(p-1)\times p$ grid as opposed to a $p\times p$ grid. All this because that is a point of view, where a well-known answer can be given. My guess is that this change would not disturb the equidistribution property. If I'm wrong then we can scratch this.
Let $\xi=e^{2\pi i/(p-1)}$ and $\zeta=e^{2\pi i/p}$. For all intagers $m$ the function $$\psi_m:\Bbb{F}_p^*\to\Bbb{C}^*, g^x\mapsto \xi^{mx}$$ is a multiplicative character of $\Bbb{F}_p$. In other words, $\psi_m(ab)=\psi_m(a)\psi_m(b)$ for all $a,b\in\Bbb{F}_p^*$. This is because $\psi_m(g)=\xi^m$ is a root of unity of order $d\mid p-1$ and in $\Bbb{F}_p^*$ we have $g^x\cdot g^y=g^{x+y}$ with the addition of the exponents done modulo $p-1$.
Similarly, for all integers $n$ the function $$ \chi_n:\Bbb{F}_p\to\Bbb{C}^*, a\mapsto \zeta^{na} $$ is an additive character of $\Bbb{F}_p.$ In other words, $\chi_n(a+b)=\chi_n(a)\chi_n(b)$ for all $a,b\in\Bbb{F}_p$.
Then the substitution $a=g^x$ converts the sums $$ \Sigma_{m,n}(p)=\sum_{x\in[0,p-2]}\xi^{mx}\zeta^{ng^x} $$ to the character sums $$ \Sigma_{m,n}(p)=\sum_{a\in\Bbb{F}_p^*}\psi_m(a)\chi_n(a). $$ These sums are one of the most studied hybrid character sums also known as Gauss sums. It is easy to show that $$ \left|\Sigma_{m,n}(p)\right|= \begin{cases} p-1,&\ \text{if $p-1\mid m$ and $p\mid n$},\\ 0,&\ \text{if $p-1\nmid m$ and $p\mid n$},\\ 1,&\ \text{if $p-1\mid m$ and $p\nmid n$, and}\\ \sqrt{p},&\ \text{when neither character is trivial.} \end{cases} $$ Only the last case is somewhat non-trivial. The divisibility constraints come from the facts that when $p-1\mid m$ we have $\psi_m(a)=1$ for all $a$, and similarly when $p\mid n$ we have $\chi_n(a)=1$ for all $a$.