Jelly, 17 bytes

FAṀŒRṗLƲṁ€ðæ*2iị⁸

Try it online!

It's been over a week, so I figured I'd reveal my brute-force Jelly approach. This makes the following assumption, which I haven't managed to prove, but can't find a counter example:

Let \$B\$ be the input \$n \times n\$ matrix, and \$A\$ be an \$n \times n\$ matrix such that \$A^2 = B\$

Let \$b = \max\{|B_{ij}|, 1 \le i,j \le n\}\$ i.e. the largest absolute element of \$B\$. This answer assumes that, assuming \$A\$ exists, there is at least one matrix \$A\$ such that all its elements are in the inclusive range \$[-b, b]\$

For example, for $$B = \left[ \begin{matrix} 9 & -8 \\ 0 & 1 \end{matrix} \right],$$ \$b = 9\$ and there exists at least one \$A\$ such that all elements are between \$-9\$ and \$9\$ (the example in the question, for example)

How it works

FAṀŒRṗLƲṁ€ðæ*2iị⁸ - Main link. Takes B on the left
F                 - Flatten B
 A                - Absolute values of each.
       Ʋ          - Last 4 links as a monad f(abs(flat(B)):
  Ṁ               -   Maximum
   ŒR             -   Bounced range; [-max, +max]
      L           -   Length i.e. number of elements of B
     ṗ            -   Powerset
         €        - Over each list:
        ṁ         -   Mold it like B
                    Call this list of matrices M
          ð       - Begin a new dyadic chain with M on the left and B on the right
           æ*2    - Matrix square of each matrix in M
              i   - Index of B in M
                ⁸ - Yield M
               ị  - Index back into M

One byte longer is a slightly more conventional way of extracting the root:

FAṀŒRṗLƲṁ€æ*2⁼ɗƇ⁸Ḣ

Try it online!

Here, the dyadic filter æ*2⁼ɗƇ means that it would save a byte over the dyadic chaining with the indexing. However, as we should only output one solution, and the Ḣ would chain to the filter without either ⁸ or ¹, the indexing here actually saves a byte.

WolframLanguage (Mathematica), 54 52 45 bytes

Last@Solve[(q=s~Array~{#,#}).q==#2,Integers]&

–2 bytes from ovs.

–7 bytes due to reading the directions and realizing I can take n as an input.

Try it online! [The request exceeded the 60 sec time limit, so it could only solve the first eight of the ten test cases.]

How it works:

q=s~Array~{#,#} //creates a dummy-indexed matrix with dimensions equal to *n* , calls it "q"
Solve[(q=s~Array~{#,#}).q==#2,Integers] //Solves q.q=#2 for q, where #2 is the input matrix, restricting the solutions to matrices with integer elements
Last@... //since there can be multiple solutions, this selects  the last one (the first ones contain conditionals)

Here I leaned heavily on the OP's allowance that output could be in "any convenient format", since the above code specifies the solution matrix in terms of its indexed elements, i.e., s[row, column]. For example:

For 3 additional bytes (48 total), a nicer output can be obtained with a small modification of the code offered by att in the comments:

(q=s~Array~{#,#})/.Last@Solve[q.q==#2,Integers]&

Julia 1.0, 109 100 bytes

I=Iterators
m%n=for i=I.countfrom(),j=I.product(fill(-i:i,n^2)...)
r=+m;r[:].=j;r^2==m&&return r end

Try it online!

9 bytes golfed by MarcMush.

Boring bruteforce solution. Takes input as matrix m and its size n. But here's also a bonus:

Julia 0.7, 214 194 bytes

~=round
m%n=.~filter(j->j≈.~j,1-2digits.(0:2^n-1,2,n).|>h->((t,z)=schur(m+0im);u=√t.*h;for i=n:-1:1,j=i+1:n a&k=a-u[i,k]*u[k,j];u[i,j]=reduce(&,t[i,j],i+1:j-1)/(u[i,i]+u[j,j])end;z)*u*z')[1]

Try it online!

This is an attempt at an algebraic solution making use of Julia's rich linear algebra library. In fact, Julia's sqrt or √ already works on matrices, but obviously only returns the principal root which may not be the one we need. Therefore, we resort to Schur decomposition, and our task boils down to finding the square roots of the upper triangular matrix \$T\$ (notation is different from Wiki, but corresponds to Julia's):

\$M^\frac{1}{2}=ZT^\frac{1}{2}Z^{-1}\$

The square roots \$U\$ of matrix \$T\$ can be found by Björck-Hammarling method, which is what Julia itself uses in the implementation of sqrt:

\$U_{ii}=T_{ii}^\frac{1}{2}\$

\$U_{ij}=\frac{T_{ij} - \sum_{k=i+1}^{j-1} U_{ik}U_{kj}}{U_{ii}+U_{jj}}, j>i\$

In order to find all solutions, we also need to apply all possible permutations of the signs of the roots in the first equation. The upper triangle can then be fixed according to the second equation.

All these computations introduce some roundoff errors, but when checking for integers, we can resolve this by using inexact comparison operator ≈. There is a problem with zero, for which it behaves the same as exact equality and comparisons fail. (It works just fine when applied to entire matrix that has both zeros and non-zeros).

The results are composed of complex numbers with zero imaginary part, so in TIO footer, we extract the real part for nicer display. It would be -2 bytes if it's acceptable to output the values without cleaning the roundoff errors, and further -3, if we can actually output all found integer solutions instead of just one.

Haskell, ⁹⁸ 96 bytes

• -2 bytes thanks to Bubbler.

s a=[x|m<-[1..],x<-q(q[-m..m]a)a,a==[foldr1(z(+))$z(map.(*))r x|r<-x]]!!0
q=mapM.const
z=zipWith

Try it online!

The relevant function is s, which takes as input a matrix a as a list of lists of integers and returns a square root of a in the same form.

The matrix multiplication part is taken almost verbatim from xnor's answer to this challenge.