A paper by Bourgain will help. Here is the [link](https://link.springer.com/article/10.1007/s00039-008-0691-6). 

Bourgain's result is:

Let $p$ be a prime and $H$ be a subgroup of the multiplicative group $\mathbb{Z}_p^*$. There exists positive constants $C_1>1$ and $C$ such that if $|H|> p^{C_1/\log\log p}$ then 
$$
\max_{(k,p)=1} \left| \sum_{a\in H} \exp\left(2\pi i k a /p \right)\right|=O\left( \exp(-(\log p)^C) |H| \right).$$

I also wrote a paper applying Bourgain' result to obtain an equidistribution result. This is Corollary 2.6 of [my paper](http://www.sciencedirect.com/science/article/pii/S0022314X16301354). 

The result is:

Let $p$ be a prime and $y\geq 1$. Let $d|p-1$ and $d> p^{C_1/\log\log p}$ for some $C_1>1$. Then there exists positive constant $C$ such that 
$$
|a<y: a^d\equiv 1 (p)|= \frac {yd}p +O(d\exp(-(\log p)^C)).
$$

The method is a standard method that is frequently used, which links equidistribution result with exponential sum bound by using Erdos Turan inequality. 

However, this result is equidistribution of the elements generated by an element of high order, not the elements of high orders.