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$$?

For large $$n$$ we may treat $$x\equiv k/n+1$$ as a continuous variable with a uniform density in the interval $$0 . 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| As a check $$\int \rho(\lambda)d\lambda=n$$.
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.