Define $n\times n$ matrices $L , M$ as follows: $\,$ L$_{ij}$= $1$ if $\,$i + j $\geq$n + 1 and is zero otherwise; $\,$ $\qquad$ $\qquad$ $\quad$M$_{ij}$= $\min(i,j)$ for $1 \leq i , j \leq n$ . It is straightforward to check that $M = L^2$ . For example, when $n=3$ we have
$\qquad$ $\qquad$ $\qquad$$\begin{bmatrix}1 & 1 & 1\\1 & 2 & 2\\1 & 2 & 3\\ \end{bmatrix}$ = $\begin{bmatrix}0 & 0 & 1\\0 & 1 & 1\\1 & 1 & 1\\ \end{bmatrix}^2$ .
[Remarks: $\det(M) = 1$ ; $M$ is positive definite but $L$ is not ; the eigenvalues of $M$ are distinct ( to see this, it's easier to work with $M^{-1}$ which is tridiagonal ). ]
General theory then predicts that $L$ can be written as a polynomial in $M$. A priori, this polynomial might be expected to have algebraic numbers as coefficients, but the examples below suggest that the situation may be nicer than that. First, when $n=2$ we find by inspection that $L = M - I.$ Next, for $n=3$, a calculation shows that $L = M^2 - 5M + 2I$. And yet again for $n=4$, we find $L = 2M^3 - 19M^2 + 21M - 5I.$
In each case (so far) the polynomial turned out to have integral coefficients.
Questions: (1) Is it always true, for arbitrary $n$ , that $L$ can be expressed as a polynomial in $M$ with integral coefficients?
(2) Presumably, matrices like $M$ above which have an integral square root are the exception rather than the rule. Is it easy to show that almost all ( in any suitable sense) matrices in $\text{SL}(n,Z)$ do not possess a square root in $\text{GL}(n,Z)$ whenever $n \geq 2$ ?
If $L$ and $M$ are $n \times n$ matrices, $g(M) = L$ for polynomial $g(z) = c_0 + \ldots + c_{n-1} z^{n-1}$ translates to the set of linear equations $\sum_{k=0}^{n-1} c_k (M^k)_{ij} = L_{ij}$ for $i,j = 1\ldots n$. If $L$ and $M$ have rational entries and solutions exist, then rational solutions exist. The only surprise is that the solutions are in integers.
If $L$ has minimal polynomial $f(\lambda)$ (in your case these appear to always be the characteristic polynomial), then $g(L^2) = L$ if and only if $f(\lambda)$ divides $g(\lambda^2) - \lambda$. If $f$ has degree $n$ and $g(\lambda) = \sum_{j=0}^{n-1} c_j \lambda^j$, we take $a_{kj}$ so that $\lambda^{2k} \equiv \sum_{j=0}^{n-1} a_{kj} \lambda^j \mod f(\lambda)$. Then $c_j$ are obtained by solving $A c = (0,1,0,\ldots,0)^T$ where $A$ is the matrix of $a_{kj}$, $k,j = 0 \ldots n-1$. In your case it appears that $A$ always has determinant $\pm 1$, which will cause $c$ to have integer entries. I have verified that the determinant is $\pm 1$ for $n$ up to $200$. That seems not to be the case for "most" matrices in $SL(n,\mathbb Z)$, and the polynomials $g$ in those cases will tend to have non-integer entries.
The characteristic polynomial $f_n(\lambda)$ of $M$ appears to satisfy a recurrence: $f_{n+4}(\lambda) + (1 - 2 \lambda^2) f_{n+2}(\lambda) + \lambda^4 f_n(\lambda) = 0$. I'm not sure if that will help.