Diagonalizing a Certain Real and Symmetric Toeplitz Matrix Consider $0\leq \alpha\leq 1$, and let $A_{\alpha}$ be the Toeplitz $n\times n$ matrix given by
$$
A_\alpha := \begin{bmatrix} 
1 & \alpha & \alpha^2 & \ldots &\alpha^{n-1} \\\ 
\alpha & 1 & \alpha & \ddots & \vdots \\\ 
\alpha^2 & \alpha & \ddots & \ddots & \alpha^2 \\\
\vdots & \ddots & \ddots  & 1 & \alpha \\\ 
\alpha^{n-1} & \ldots  & \alpha^2 & \alpha & 1 
\end{bmatrix}. 
$$
We can decompose $A_{\alpha}=U_{\alpha}D_{\alpha}U_{\alpha}^{*}$ where $U_{\alpha}$ is a unitary matrix and $D_{\alpha}$ is a diagonal matrix. 


*

*Are the matrices $U_{\alpha}$ and
$D_{\alpha}$ explicitly known as a
function of $\alpha$ and $n$? Can
anyone point me to the right
reference.

*What is known about the spectral
limit distribution of these matrices
as $n\to\infty$? More specifically,
can we compute the limit moments $$
   \gamma_{k}:=\lim_{n\to\infty}{\frac{1}{n}\mathrm{Tr}\Big(A_{\alpha}^{k}\Big)}
   $$ for $k\geq 1$?
Thanks!
 A: While there seems to be no explicit diagonalization, here is a simple idea to obtain an approximation.
To simplify notation, define $T=A_x$ to be the original Toeplitz matrix (I'm using $x$ instead of $\alpha$ to save typing). 
After some playing around, it can be seen that the inverse enjoys remarkable structure, namely,
$$
  T^{-1} = \frac{1}{1-x^2}
  \begin{bmatrix}
    1  & -x & \cdots & \cdots & 0\\\\
    -x & 1+x^2 & -x & \cdots & 0\\\\
    & \ddots & \ddots & \ddots\\\\
    0 & \cdots & &1+x^2 & -x\\\\
    0 & \cdots & &-x & 1
  \end{bmatrix}
$$
Consider, therefore, the following Toeplitz matrix
\begin{equation*}
M :=
\begin{bmatrix}
  1+x^2  & -x & \cdots & \cdots & 0\\\\
  -x & 1+x^2 & -x & \cdots & 0\\\\
  & \ddots & \ddots & \ddots\\\\
  0 & \cdots & &1+x^2 & -x\\\\
  0 & \cdots & &-x & 1+x^2
  \end{bmatrix},
\end{equation*}
for which one has closed form eigenvalues and eigenvectors, given by
\begin{equation*}
  \lambda_k = (1+x^2)-2x\cos\left(\frac{k\pi}{n+1}\right),\quad 1 \le k \le n,
\end{equation*}
and
\begin{equation*}
v_{ik} = \sin\left(\frac{ik\pi}{n+1}\right),\quad 1 \le i \le n, 1 \le k \le n.  
\end{equation*}
These eigenvalues (after scaling by $1-x^2$) and eigenvectors may be approximately substituted for those of $T^{-1}$.
Simple experimentation reveals that the eigenvalues of the matrix
\begin{equation*}
  M' := 
  \begin{bmatrix}
    1  & -x & \cdots & \cdots & 0\\\\
    -x & 1+x^2 & -x & \cdots & 0\\\\
    & \ddots & \ddots & \ddots\\\\
    0 & \cdots & &1+x^2 & -x\\\\
    0 & \cdots & &-x & 1
  \end{bmatrix}
\end{equation*}
satisfy $$|\lambda(M) - \lambda(M')| \le 4x/n,$$
where the bound can be made tighter by closer analysis (note that $\lambda(M) \ge \lambda(M')$ also holds). Similar results can also be shown for the eigenvectors, but I haven't had the time to prove that. 


EDIT
This journal article shows how to compute eigenvalues and eigenvectors for matrices that look like $T^{-1}$ above. In particular, it implies for example, that the explicit eigenvalues of $T^{-1}$ above are given by 
$$\lambda_k = 1 - 2x\cos\theta_k,$$
where, for $k=1,\ldots,n$, the angle $\theta_k$ is a root of
$$\sin(n+1)\theta + x^2\sin(n-1)\theta - x\sin(n\theta).$$
Formulae for eigenvectors can also be found in terms of the $x$ and $\theta_k$ as stated above.
However, it seems that in our case, we'll have to numerically solve for $\theta_k$. Modulo that, I guess, this is as close as we'll get to explicit eigenvalues and eigenvectors of $T^{-1}$ (and thereby of $A_x$).
A: 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.
A: This class of matrices may have additional nice properties.  You might look at the minors when alpha is less than 1; I am thinking in this case, you may have all positive minors, in other words, a totally positive matrix.   
