Limit of a sequence of determinants.
We answer this question using the information in the addendum. We proceed in three steps. (1) We compute the eigenvalues of each $M_n$ and notice that the problem reduces to evaluating a limit of a sequence of products. (2) We show that the product formula proven in another math.SE question can be used to evaluate each term in this sequence of products. (3) We take the $n\to\infty$ limit of the resulting expression.
Step 1.
For each $n\geq 2$, let $U_n = ((U_n)_{ij}) $ denote the $n\times n$ matrix defined by \begin{align} (U_n)_{ij} = \left\{\begin{array}{cc} 1, & j=i+1 \\ 1, & (i,j) = (n,1) \\ 0, & \mathrm{otherwise} \end{array}\right. \end{align} then we have \begin{align} M_n = (2+\epsilon^2)I_n-U_n-U_n^\dagger \end{align} where $I_n$ is the $n\times n$ identity, and $\dagger$ denotes complex conjugate + transpose (which in this case is of course just the transpose since the matrices are real). Let's suppress the subscript $n$'s for notational simplicity for a moment.
It is straightforward to show that each $U_n$ has the property $(U_n)^n = I$; it's $n^\mathrm{th}$ matrix power is the identity. This tells us that each eigenvalue $\lambda$ satisfies $\lambda^n = \lambda$, namely its eigenvalues are the $n^\mathrm{th}$ roots of unity $e^{2\pi ik/n}$ with $k=0,1,\dots, n-1$ with corresponding eigenvectors $v_k$. It is also not hard to see that $U_n^\dagger$ shares the same eigenvectors $v_k$ with the complex conjugated corresponding eigenvalues $e^{-2\pi ik/n}$. It follows that the eigenvalues of $M_n$ are \begin{align} \lambda_k = 2+\epsilon^2 - e^{2\pi ik/n} - e^{-2\pi ik/n} = 4\sin^2\left(\frac{\pi k}{n}\right)+\epsilon^2 \end{align} The determinant of $M_n$ is the product of its eigenvalues, so we obtain \begin{align} \Delta_n = \prod_{k=0}^{n-1}\left[4\sin^2\left(\frac{\pi k}{n}\right)+\left(\frac{\beta}{n}\right)^2\right] \end{align} as claimed in the addendum. So problem reduces to evaluating the limit of a sequence of finite products; \begin{align} \Delta = \lim_{n\to \infty}\left(\prod_{k=0}^{n-1}\left[4\sin^2\left(\frac{\pi k}{n}\right)+\left(\frac{\beta}{n}\right)^2\right]\right). \end{align}
Step 2.
The following product identity holds: \begin{align} \prod_{k=0}^{n-1}\left[\sinh^2y+\sin^2\left(x+\frac{k\pi}{n}\right)\right]=2^{1-2n}(\cosh(2ny) -\cos(2nx)) \end{align} see, for example, the following math.SE question (which, incidentally, was inspired by the question as hand):
How would one discover this finite product identity?
If we set $x=0$, $y=\sinh^{-1}(\beta/(2n))$ and multiply this product by $4^n$, then we obtain the following identity as a special case: \begin{align} \prod_{k=0}^{n-1}\left[4\sin^2\left(\frac{\pi k}{n}\right)+\left(\frac{\beta}{n}\right)^2\right] = 2\left\{\cosh\left[2n\sinh^{-1}\left(\frac{\beta}{2n}\right)\right]-1\right\} \end{align} The desired expression is thus the $n\to\infty$ limit of the expression on the right hand side.
Step 3.
For each $\beta >0$, the large $n$ limit allows us to approximate $\sinh^{-1}(\beta/(2n)) \sim \beta/(2n)$ so the expression on the right hand side goes to $2(\cosh(\beta)-1)$ as $n\to\infty$. Now simply recall the half-angle formula $\sinh(\theta/2) = ((\cosh\theta-1)/2)^{1/2}$ which gives the conjectured result; \begin{align} \Delta = 4\sinh^2\left(\frac{\beta}{2}\right). \end{align}