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: Determining a generating function (of a restricted form)).
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.
