WolframLanguage (Mathematica), 545245 4445 bytes

#&@@Solve[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.

-1 byte from att.

Try it online!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 arbitrarilyselects chooses the firstlast 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. This gives a different solution from mine because it selects the last solution rather than the first:

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

Last@Solve[#&@@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!-1 byte from att.

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, Lastthis arbitrarily chooses the lastfirst one

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. This gives a different solution from mine because it selects the last solution rather than the first:

(q=s~Array~{#,#})/.Last@Solve[q.q==#2,Integers]& WolframLanguage (Mathematica), 5452 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, Last arbitrarily chooses the last one

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 solutioncode offered by att in the comments:

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