As we have known, the Sato-Tate measure for GL(2) turned out to be the half circle measure
$\frac{1}{2\pi} \sqrt{4-x^2}dx$ on [-2,2],
which appears in various versions of equi-distribution problems in GL(2).
My question is what the corresponding measure for GL(3) should be.
Firstly we shall note that Hecke eigenvalues in GL(3) can be imaginary. And we have to replace [-2,2] with {$z\in \mathbb{C}:|z|<=3$} the ball of radius 3 in complex plane. I don't know anything further yet.
(2017-11-26 edit by j.c.: earlier versions of this answer consisted of David Hansen's screenshot of the following, with the text "Here is a screenshot of a semi-answer which froze my computer when I hit 'post'":)
Things get more complicated. $\newcommand{\SU}{\mathop{\rm SU}\nolimits} \newcommand{\C}{\bf C} \newcommand{\Irr}{\mathop{\rm Irr}\nolimits} $
Let $\pi=\pi_\infty\otimes\bigotimes'_p\pi_p$ be a unitarily normalized cuspidal automorphic representation of $GL_n(\bf{A}_{\bf Q})$; suppose for simplicity that $\pi$ has trivial central character. The local factor $\pi_p$ (for $p$ outside a finite set $S$) is determined by a matrix $s(\pi_p)\in GL_n(\C)$, which is well-defined up to conjugacy. The generalized Ramanujan conjecture predicts that $s(\pi_p)$ may be chosen as an element of the unitary group $\SU_n$. Suppose this is true, and let $f:\SU_n\rightarrow\C$ be a continuous function which is invariant under conjugation. Here is a generalized Sato-Tate conjecture:
If $\pi$ is not a functorial lift from a "smaller" group, then $\frac{1}{\pi(X)}\sum_{p\leq X}f(s(\pi_p))\rightarrow \int_{\SU_n}f(g)dg$ where $dg$ is the Haar probability measure on $\SU_n$.
If $n=2$ and $\pi$ is associated with an elliptic curve $E$, this simplifies drastically: the conjugacy class $s(\pi_p)\in\SU_2$ is determined uniquely by its trace, which is $\frac{a_p(E)}{\sqrt{p}}$. Likewise, any $f$ determines a unique continuous $f':[-2,2]\rightarrow\C$, with $f=f'\circ\mathrm{tr}$. The Haar probability measure $dg$ pushes forward to the measure $\frac{1}{2\pi}\sqrt{4-x^2}$ on $[-2,2]$, and the conjecture now takes the more familiar form.
What's the meaning of the vague condition "$\pi$ is not a functorial lift from a smaller group"? Well, let's try to give a plausibility argument for the above conjecture. Let $\Irr_n$ denote the set of irreducible algebraic representations of $\SU_n$; for $\sigma\in\Irr_n$, write $\chi_\sigma$ for its character. In $L^2(\SU_n)$ we have
$$f(g)=\sum_{\sigma\in\Irr_n}\chi_\sigma(g)\langle f,\chi_\sigma\rangle.$$
Changing $g$ to $s(\pi_p)$, summing over $p\leq X$ and switching the order of summation gives
$$\frac{1}{\pi(X)}\sum_{p\leq X}f(s(\pi_p)) \approx \int_{\SU_n}f(g)dg + \sum_{\sigma\text{ nontrivial}}\frac{1}{\pi(X)}\sum_{p\leq X}\chi_\sigma(s(\pi_p)).$$
The behavior of the sum $\sum_{p\leq X}\chi_\sigma(s(\pi_p))$ is ontrolled by the Langlands L-function $L(s,\pi,\sigma)$, and in particular this sum is $o(X)$ if the L-function is nonvanishing and without poles in $\mathrm{Re}(s)\geq 1$. The truth of this latter property for all $\sigma$ is expected to be equivalent to $\pi$ not arising via functoriality from a smaller group.
Based on the comments I'm going to assume that the measure you're interested in is the pushforward of Haar measure to the conjugacy classes of $\text{SU}(3)$ (I know nothing about why you should expect this). There are a lot of things you can say about this measure without writing down an explicit formula for it. In particular, you can compute the moments of various random variables with respect to this measure.
The basic observation is that the map sending a (complex, unitary, finite-dimensional) representation of $\text{SU}(3)$ to its character induces a homomorphism $\chi$ from the representation ring of $\text{SU}(3)$ to the algebra of class functions on $\text{SU}(3)$, and by Peter-Weyl the span of the image is dense. Moreover, we know the integral over $\text{SU}(3)$ of the character of every representation: it's the dimension of the invariant subspace. So to compute moments of class functions on $\text{SU}(3)$ (random variables on the conjugacy classes) it suffices in principle to understand invariant subspaces of tensor products of representations of $\text{SU}(3)$ (given that you know how to express your class functions in terms of characters).
For example, suppose you're interested in the moments of twice the real part of the trace of a random element of $\text{SU}(3)$. This is the character of the representation $W = V \oplus V^{\ast}$ where $V$ is the defining representation. The $n^{th}$ moment is then precisely the dimension of the invariant subspace of $W^{\otimes n}$, and to compute this it suffices to count the number of walks of length $n$ in a Weyl chamber from the origin to itself by highest weight theory. The first few moments are $$1, 0, 2, 2, 12, 30, 130, ....$$
This is A151366 in the OEIS. There is a generating function given there but it involves hypergeometric functions, so like Marty I don't expect a particularly nice explicit formula for the corresponding measure.
(Doing the above for $\text{SU}(2)$ you reduce the problem to computing the moments of a single random variable, the trace. The corresponding representation is just the defining representation, the Weyl chamber is $\mathbb{Z}_{\ge 0}$, and walks of length $n$ in the Weyl chamber from the origin to itself are counted by Catalan numbers, which are the moments of the Sato-Tate measure.)
Edit: Emerton says in the comments that the trace of an element of $\text{SU}(3)$ takes values in $[-3, 3]$ and doesn't determine the conjugacy class. In fact, it takes values in the ball of radius $3$ in $\mathbb{C}$ and it does determine the conjugacy class!
Recall that $g \in \text{SU}(3)$ has eigenvalues $\lambda_1, \lambda_2, \lambda_3$ which are unit complex numbers satisfying $\lambda_1 \lambda_2 \lambda_3 = 1$. It follows that $$\lambda_1 \lambda_2 + \lambda_2 \lambda_3 + \lambda_3 \lambda_1 = \overline{\lambda_1 + \lambda_2 + \lambda_3}$$
hence that $\text{tr}(g) = \lambda_1 + \lambda_2 + \lambda_3$ determines the characteristic polynomial of $g$, hence its conjugacy class. So to understand the pushforward of Haar measure to the conjugacy classes of $\text{SU}(3)$ it suffices to understand the joint moments of the real and imaginary parts of the trace. For example, the moments of twice the imaginary part are (up to sign) the dimensions of the invariant subspaces of tensor powers of the virtual representation $V \ominus V^{\ast}$ and so one must count walks in a Weyl chamber with certain weights $\pm 1$. The first few moments are $$1, 0, 2, 0, 12, 0, 98, ...$$
and this sequence does not appear to be in the OEIS.
I can give a closed formula for the pushforward of $SU(3)$ Haar measure under trace. I don't understand the number theoretic issues well enough to know whether I want $SU(3)$ or $U(3)$, though.
I'll write $(e^{i \alpha}, e^{i \beta}, e^{i \gamma})$ for the eigenvalues of the unitary matrix, with $\alpha+\beta+\gamma=0$, and write $z=x+iy$ for the trace. Set $$\Delta = \left(e^{i \alpha} - e^{i \beta}\right) \left( e^{i \alpha} - e^{i \gamma}\right) \left( e^{i \beta} - e^{i \gamma}\right)$$ so $$|\Delta| = 8 \sin \left( \frac{\alpha-\beta}{2} \right) \sin \left( \frac{\alpha-\gamma}{2} \right) \sin \left( \frac{\alpha-\gamma}{2} \right) .$$ By the Weyl integration formula, conjugacy class $(\alpha, \beta, \gamma)$ has volume proportional to $|\Delta|^2$.
We have $$\begin{array}{rcl} z &=& e^{i \alpha} + e^{i \beta} + e^{i \gamma} \\ x &=& \frac{z+\bar{z}}{2} \\ y &=& \frac{z - \bar{z}}{2i} \\ \end{array}$$ so $dx dy$ is $\frac{1}{2i} dz \bar{dz}$. We compute
$$ \begin{pmatrix} \frac{\partial z}{\partial \alpha} & \frac{\partial z}{\partial \beta} \\ \frac{\partial \bar{z}}{\partial \alpha} & \frac{\partial \bar{z}}{\partial \beta} \\ \end{pmatrix} = \begin{pmatrix} e^{i \alpha} - e^{i \gamma} & e^{i \beta} - e^{i \gamma} \\ e^{-i \alpha} - e^{-i \gamma} & e^{- i \beta} - e^{-i \gamma} \\ \end{pmatrix}.$$ So the Jacobian of $(\alpha, \beta) \mapsto (z,\bar{z})$ is $$ \begin{array}{rcl} (e^{i \alpha} - e^{i \gamma})(e^{-i \beta} - e^{-i \gamma}) - (e^{-i \alpha} - e^{-i \gamma})(e^{i \beta} - e^{i \gamma}) &=& (e^{i \alpha} - e^{i \gamma})(e^{i \beta} - e^{i \gamma})(-e^{-i (\beta+\gamma)} + e^{-i(\alpha+\gamma)})\\ &=& (e^{i \alpha} - e^{i \gamma})(e^{i \beta} - e^{i \gamma})(- e^{i \alpha} + e^{i \beta}) \\ &=& - \Delta \\ \end{array}$$
I think I might have lost a sign somewhere, but the point is that $\Delta dx dy$ is proportional to $d \alpha d \beta$. So the Weyl measure, $|\Delta^2| d \alpha d \beta$, is proportional to $|\Delta| dx dy$. (If you really want to practice your high school trig, try doing that computation with $x$ and $y$ rather than $z$ and $\bar{z}$.)
Now, $\Delta^2$ is the discriminant of the cubic $$t^3 - (x+iy) t^2 + (x-iy) t - 1.$$ Mathematica tells me that is $$-27 + 18 x^2 - 8 x^3 + x^4 + 18 y^2 + 24 x y^2 + 2 x^2 y^2 + y^4$$
So the desired measure on $\mathbb{C}$ is proportional to $$\sqrt{27 - 18 x^2 + 8 x^3 - x^4 - 18 y^2 - 24 x y^2 - 2 x^2 y^2 - y^4} dx dy.$$ The region where the quantity in the square root is nonnegative is exactly the deltoid John Baez describes.
In principal, one should be able to carry out this computation carefully enough to get the constant, but I'm not going to do it.
Conceptually, this computation is similar to the classical computation that the map sending $(\alpha_1, \ldots, \alpha_r)$ to the coefficients of the polynomial $\prod (z-\alpha_i)$ has Jacobian $\prod_{i<j} (\alpha_i - \alpha_j)$. However, because of the constraint that the eigenvalues have product $1$, the details don't quite match.
