Eigenvalue density of a symmetric tridiagonal matrix Let $A_n\in\mathbb{R}^{n\times n}$ be defined as
$$
A_n=\begin{bmatrix} a & b & 0 & \cdots & \cdots & 0 & 0\\ b & a & b & \cdots & \cdots & 0 & 0\\ 0 & b & a & \cdots & \cdots & 0 & 0\\ \vdots & \vdots & \vdots &\ddots & \ddots & \vdots & \vdots \\ \vdots & \vdots & \vdots &\ddots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 &\cdots & \cdots & a & b \\ 0 & 0 & 0 &\cdots & \cdots & b & a\end{bmatrix},
$$
where $a,b\in\mathbb{R}$. It is well-known that the eigenvalues of $A_n$ are
$$
\text{eig}(A_n) =\left\{ a+2b\cos\left(\frac{\pi}{n+1}k\right), \ k=1,2,\dots,n \right\}.
$$

My question. Does there exist a closed-form expression for the eigenvalue density of the sequence $\{A_n\}$ as $n\to \infty$?

 A: For large $n$ we may treat $x\equiv k/n+1$ as a continuous variable with a uniform density in the interval $0 <x< 1$. The corresponding eigenvalue $\lambda(x)=a+2b\cos\pi x$ ranges from $a-2|b|$ to $a+2|b|$. Since
$$|d\lambda/dx|=\pi\sqrt{4b^2-(\lambda-a)^2}.$$
The eigenvalue density follows from
$$\rho(\lambda)d\lambda=ndx\Rightarrow \rho(\lambda)=\frac{n}{\pi}\frac{1}{\sqrt{4b^2-(\lambda-a)^2}},\;\;a-2|b|<x<a+2|b|.$$
As a check $\int \rho(\lambda)d\lambda=n$.
A: To put this into context, the limit
$$
\lim_{L\to\infty} \frac{\# \textrm{ eigenvalues in }I \textrm{ of the problem on } \{0,\ldots, L\} }{L}
$$
(assuming it exists) is one way of defining the density of states measure $\int_I dN(\lambda)$.
For an ergodic (with respect to the shift) system $\mathcal A$ of operators with probability measure $dP$ this will equal the average
$$
dN(\lambda)=\int_{\mathcal A} d\mu(\lambda;A)\, dP(A)
$$
of the spectral measures $d\mu$.
In your case, for constant coefficients, the system consisting of this single operator is trivially ergodic, so the density of states is the spectral measure
$$
d\mu(\lambda) = \chi_{(a-|b|,a+|b|)}(\lambda) \frac{d\lambda}{\sqrt{4b^2-(\lambda -a)^2}} ,
$$
which Carlo finds by direct computation.
