There's likely no explicit diagonalization of $A_\alpha$ except when
$n$ is very small or in special cases like $\alpha = 0$ and $\alpha
= 1$.  Nevertheless each "limit moment" $\gamma_k$ can be computed as
a rational function of $\alpha$, and this can be used to describe for
each $\alpha$ the distribution of eigenvalues of $A_\alpha$ as $n \rightarrow \infty$.

To diagonalize $A_\alpha$ explicitly we'd need to know the eigenvalues;
these are roots of the degree-$n$ characteristic polynomial $\chi_{A_\alpha}$,
and it's often too much to expect that a family of such polynomials
can be factored for each $n$.  Here $\chi_{A_\alpha}$ does split into
two factors $\chi^\pm_{A_\alpha}$ of equal or nearly equal degree, but
usually that's as far as we can go.  The factorization arises because $A_\alpha$ commutes
with the involution, call it $\iota$, that takes each coordinate $a_k$
to $a_{n+1-k}$, so the $\pm1$ eigenspaces of $\iota$ are invariant subspaces
of $A_\alpha$.  The factor $\chi^\pm_{A_\alpha}$ is the characteristic
polynomial of the restriction of $A_\alpha$ to the $\pm1$ subspace.
 But once $n$ is at all large it seems there's nothing to be done with
these factors  $\chi^\pm_{A_\alpha}$.  For example, trying "random"
rational values for $\alpha$ yields polynomials whose Galois group is
the full symmetric group.  Thus if you ask **gp**

    f(a,n) = factor(charpoly(matrix(n,n,i,j,a^abs(i-j))))
    F = f(1/2,21)
    vector(#F[,1], n, polgalois(F[n,1]))

you'll see that for $n=21$ the factors of $A_{1/2}$ have Galois groups
$S_{10}$ and $S_{11}$.

There are some special values of $\alpha$ for which one can find the
roots of $\chi_{A_\alpha}$ explicitly.  Most obviously, $A_0$ is the
identity matrix.  Also $A_1$ is the all-ones matrix, with one eigenvalue
of $N$ and all other eigenvalues zero.  The OP required $\alpha \in
[0,1]$, but $A_{-1}$ has rank 2 so its eigenvalues are easy too.  In
each of these cases there's no unique diagonalization because there's
an eigenvalue with high multiplicity.

As for the limit moments $\gamma_k$: if $\alpha=1$ then clearly
${\rm Tr}(A_\alpha^k) = n^k$ so $\gamma_k=\infty$ once $k>1$.
So we assume $\alpha < 1$, and then we may as well take
$\alpha \in {\bf C}$ with $|\alpha| < 1$.  Then
$\gamma_1$, $\gamma_2$, $\gamma_3$, $\gamma_4$, $\gamma_5$, etc. are
$$
1, \
\frac{1+\alpha}{1-\alpha},\
\frac{1+4\alpha+\alpha^2}{(1-\alpha)^2},\
\frac{1+9\alpha+9\alpha^2+\alpha^3}{(1-\alpha)^3},\
\frac{1+16\alpha+36\alpha^2+16\alpha^3+\alpha^4}{(1-\alpha)^4}, \ldots
$$
and in general $\gamma_k = P_{k-1}(\alpha) / (1-\alpha)^{k-1}$ where
$$
P_m(X) := \sum_{j=0}^m \left({m \atop j}\right)^2 X^j
$$
is the polynomial obtained from the binomial expansion of $(1+X)^m$
by squaring each coefficient.  These $P_m$ don't have an entirely elementary
formula, but they can be written as hypergeometric polynomials, or
(if memory serves) expressed in terms of Legendre polynomials,
or manipulated using the generating function
$$
\sum_{m=0}^\infty P_m(X) t^m = \left((\alpha-1)^2 t^2 - 2(\alpha+1)t
+ 1\right)^{-1/2}
$$
if I did this right (I guessed the formula using the technique I described
here a few weeks ago:
https://mathoverflow.net/questions/66863).

To get that formula for $\gamma_k$, we first find an integral representation,
which I gather is a special case of the "Szegő-Tyrtyshnikov-Zamarashkin-Tilli
theorem" that F. Poloni mentioned in his comment.  While general Toeplitz
matrices cannot be diagonalized explicitly, circulant ones can.
So we compare $A_\alpha$ with the circulant matrix
$A'_\alpha$ whose $(i,j)$ entry is $\alpha^{\min(|i-j|,n-|i-j|)}$.
For each $\alpha$ and $k$, the $k$-th powers of $A_\alpha$ and
$A'_\alpha$ differ by $O(1)$ as $n \rightarrow \infty$.
*[This was somewhat annoying to check; maybe there's a nice way to do it.]*
Thus $A_\alpha$ and $A'_\alpha$ have the same limit moments --
and the moments of $A'_\alpha$ can be computed by finding its eigenvalues.
Being circulant, $A'_\alpha$ is explicitly diagonalized by
the discrete Fourier transform on ${\bf Z} / n {\bf Z}$, with an eigenvalue
$\lambda_z = \sum_{j=0}^{n-1} \alpha^{\min(j,n-j)} z^j$
for each $n$th root of unity $z = \exp(2\pi i r/n)$.
For large $n$ we can approximate $\lambda_z$ by
$$
f_\alpha(z) = \sum_{j=-\infty}^\infty \alpha^{|j|} z^j
= \frac1{1-\alpha z} + \frac{\alpha z^{-1}} {1 - \alpha z^{-1}}
= \frac{1-\alpha^2}{ (1-\alpha z)(1 - \alpha z^{-1}) }
$$
and deduce that
$$
\gamma_k = \frac{1}{2\pi} \int_{-\pi}^{\pi} f_\alpha(e^{i\theta})^k d\theta.
$$
This also means that as $n \rightarrow \infty$ the eigenvalues of $A_\alpha$
tend to the same distribution as the image of the uniform distribution on the
unit circle $|z|=1$ under $f_\alpha$, in the sense that for any continuous
function $\phi$ on ${\bf C}$ the average of $\phi(\lambda)$ over the
eigenvalues approaches
$(2\pi)^{-1} \int_{-\pi}^{\pi} \phi(f_\alpha(e^{i\theta}))d\theta$.
For real $\alpha$ in $(0,1)$, this distribution is supported on the
interval $((1+\alpha)/(1-\alpha), (1-\alpha)/(1+\alpha))$.

In our case we can evalute the integral for $\gamma_k$ by writing it as
the contour integral

$$
\frac1{2\pi i}\oint_{|z|=1} f_\alpha(z)^k \frac{dz}{z}.
$$

For each $k \geq 1$ the integrand has a pole of order $k$ at $z = \alpha$
and no other poles in $|z| \leq 1$; evaluating the residue at this pole yields
the formula $\gamma_k = P_{k-1}(\alpha) / (1-\alpha)^{k-1}$ given above.