2 Added analysis of triple and higher coincidences, including explanation of 0.672 for 5th powers

[Edited to describe triple and higher-order coincidences for prime $k$, recovering the observed $0.672$ proportion for $k=5$]

To see thisExtending the analysis to triple and higher-order coincidences also yields the asymptotic proportion of $k$-th powers that arise in $\lbrace a^k \phantom. \bmod p : a < p/k \rbrace$. For example, note when $k$ is an odd prime, the proportion of $k$-th powers that do not have a $k$-th root in $(0,p/k)$ is asymptotic to $((k-1)^k+1)/k^k$; for $k=5$ that's $41/125$, so the proportion with such a $k$-th root is $84/125$, which matches A.Caicedo's observed $0.672$ exactly. It also gives $1 - \frac{8+1}{27} = 2/3$ for $k=3$, matching the proportion of cubes reported by Greg Martin in comments below; as $k \rightarrow \infty$ the proportion of $k$-th powers with small $k$-th roots approaches $1 - (1/e)$.

Here's how to estimate the number of pairs. Begin with the observation that $a^k = b^k$ iff $b \equiv ma \bmod p$ where $m$ is one of the $k-1$ solutions of $m^k \equiv 1 \bmod p$ other than $m=1$. If $k$ is even, we exclude also $m=-1$, which is impossible with $0<a,b<p/k$. Then $b \equiv ma \bmod p$ defines a lattice of index $p$ in ${\bf Z}^2$ all of whose nonzero vectors have length $\gg p^{\epsilon(k)}$, because for such a vector $p$ divides the nonzero number $a^k-b^k$, which factors into homogeneous polynomials in $a,b$ each of degree at most $\phi(k)$. [This is where we use $m \neq -1$: if $a=-b$ then $a^k-b^k=0$.] Thus the solutions of $b \equiv ma \bmod p$ with $a,b \in (0,p/k)$ are the lattice points in a square of area $(p/k)^2$, and their number is estimated by $p^{-1} (p/k)^2 = p/k^2$, with an error bound proportional to (perimeter)/(shortest perimeter)/(length of shortest nonzero vector), i.e. proportional to $p^{1-\epsilon(k)}$. The total of $C_k \phantom. p + O_k(p^{1-\epsilon(k)})$ then follows by summing over all $k-1$ or $k-2$ solutioms solutions of $m^k=1 \bmod p$ other than $m = \pm 1$, and dividing by 2 because we've counted each coincidence twice, as $(a,b)$ and $(b,a)$.

Likewise one can estimate the counts of triples etc. One must be careful with subsets of the $k$-th roots of unity that have integer dependencies, but at least when $k$ is prime there are no dependencies except that all $k$ of them sum to zero. If I did this right, the result for $j<k$ is that the number of $j$-element subsets of $\lbrace 1, 2, \ldots, (p-1)/k \rbrace$ with the same $k$-th power is asymptotic to ${k \choose j} p / k^{j+1}$, while there are no such subsets with $j=k$ because the sum of all $k$ solutions of $a^k \equiv c \bmod p$ vanishes. An exercise in generatingfunctionological inclusion-exclusion then produces the formula $((k-1)^k+1)/k^k$ for the asymptotic proportion of $k$-th powers that have no $k$-th roots at all in $(0,p/k)$.

1

Darij's pretty argument, extended by GH, nicely answers the question for $k$-th powers modulo a large prime $p \equiv 1 \bmod k$ for each fixed $k>2$. Yet more can be said: that approach yields the existence of one coincidence $a^k \equiv b^k$ with $0 < a < b < p/k\phantom.$; but in fact the number of coincidences is asymptotically proportional to $p$: the count is $C_k \phantom. p + O_k(p^{1-\epsilon(k)})$, where $C_k = (k-1)/(2k^2)$ or $(k-2)/(2k^2)$ according as $k$ is odd or even, and $\epsilon(k) = 1/\varphi(k) \geq 1/(k-1)$.

To see this, note that $a^k = b^k$ iff $b \equiv ma \bmod p$ where $m$ is one of the $k-1$ solutions of $m^k \equiv 1 \bmod p$ other than $m=1$. If $k$ is even, we exclude also $m=-1$, which is impossible with $0<a,b<p/k$. Then $b \equiv ma \bmod p$ defines a lattice of index $p$ in ${\bf Z}^2$ all of whose nonzero vectors have length $\gg p^{\epsilon(k)}$, because for such a vector $p$ divides the nonzero number $a^k-b^k$, which factors into homogeneous polynomials in $a,b$ each of degree at most $\phi(k)$. [This is where we use $m \neq -1$: if $a=-b$ then $a^k-b^k=0$.] Thus the solutions of $b \equiv ma \bmod p$ with $a,b \in (0,p/k)$ are the lattice points in a square of area $(p/k)^2$, and their number is estimated by $p^{-1} (p/k)^2 = p/k^2$, with an error bound proportional to (perimeter)/(shortest nonzero vector), i.e. proportional to $p^{1-\epsilon(k)}$. The total of $C_k \phantom. p + O_k(p^{1-\epsilon(k)})$ then follows by summing over all $k-1$ or $k-2$ solutioms of $m^k=1 \bmod p$ other than $m = \pm 1$, and dividing by 2 because we've counted each coincidence twice, as $(a,b)$ and $(b,a)$.

The same technique also works for $0 < a < b < M$ with $M$ considerably smaller than $p/k$; and the resulting coincidences, when they exist, can be calculated efficiently using lattice basis reduction (which as it happens I mentioned on this forum a few days ago).