R, 301 287 277 274 222 217 195 186 178 174 bytes

Nothing particularly creative, including the zero buffering of the peripheral elements of the entry matrix, an earlier version later improved by Robin:

function(x)w=which.max
if(any(s<-!x^.5%%1))
y=cbind(NA,rbind(NA,x,NA),NA)
z=y[i]=y[i<-w(y*y%in%x[s])]^.5
m=i+c(r<--c(1,nrow(y)),-r)
y[j]=y[j<-m[w(-y[m])]]*z
x=p(y[r,r])
x

Try it on-line

Using a sequence of numbers as its entry, and hence removing the call to a function, Nick Kennedy earlier managed a 186 bytes version of the algorithm as follows (with -10 bytes by Robin):

w=which.max;`~`=cbind;x=scan();while(any(s<-!x^.5%%1))y=NA~t(NA~matrix(x,n<-length(x)^.5)~NA)~NA;i=w(y*y%in%x[s]);=i+c(r<--c(1,n+2),-r);y[j]=y[j<-m[w(-y[m])]]*(y[i]=y[i]^.5);x=y[r,r];x

avoiding the definition of a (recursive) function, plus other nice gains.

Try it on-line

Ruby, 140 135 bytes

Takes a flat list as input, outputs a flat list.

->mi=1;(i=m.index m.rejecte.max
m[i+[1,-1,l=m.size**0.5,-l].min_byj]*=m[i]**=0.5if i)while i;m

Try it online!

Explanation:

->m # Anonymous lambda
 i=1; # Initialize i for the while loop
 ( # Start while loop

i=m.index # Get index at...
 m.reject # Get all elements of m, except the ones with...
 e**0.5%1>0 # a square root with a fractional component
 .max # Get the largest of these

m[i+ # Get item at...
 [1,-1,l=m.size**0.5,-l] # Get possible neighbors (up, down, left, right)
 .min_by # If out of range, return matrix max so
 # neighbor isn't chosen
 ]
 *=m[i]**=0.5 # Max square becomes its square root, then multiply
 # min neighbor by it

)while i # End while loop. Terminate when index is nil.
m # Return matrix.

Python 2, 188 bytes

M=input()
l=int(len(M)**.5)
try:
 while 1:m=M.index(max(i**.5%1or i for i in M));_,n=min((M[m+i],m+i)for i in m/l*[-l]+-~m%l*[1]+[l][:m/l<l-1]+m%l*[-1]);M[m]**=.5;M[n]*=M[m]
except:print M

Try it online!

Full program. Takes input and prints as a flat list.

Perl 6, 236 bytes

my@k=.flat;my n=$_;loop my (i,j)=@k>>.sqrt.grep(0==$_,:kv).rotor(2).max(*[1]);last if 0>i;$/=((0,1),(0,-1),(1,0),(-1,0)).map($!=i+n*.[0]+.[1];+$!,n>.[0]+i/n&.[1]+i%n>=0??@k[$!]!!Inf).min(*[1]);@k[i,$0]=j,j*$1;@k.rotor(+n)

Try it online!

MATL, 49 48 bytes

`ttX^tt1~*X>X>XJt?wy=(tt5M1Y6Z+*tXzX<=*Jq*+w}**

Try it online! Or verify all test cases.

How it works

` % Do...while
 tt % Duplicate twice. Takes a matrix as input (implicit) the first time
 X^ % Square root of each matrix entry
 tt % Duplicate twice
 1~ % Modulo 1, negate. Gives true for integer numbers, false otherwise
 * % Multiply, element-wise. This changes non-integers into zero
 X>X> % Maximum of matrix. Gives maximum integer square root, or zero
 XJ % Copy into clipboard J
 t % Duplicate
 ? % If non-zero
 wy % Swap, duplicate from below. Moves the true-false matrix to top
 = % Equals, element-wise. This gives a matrix which is true at the
 % position of the maximum that was previously identified, and
 % false otherwise
 ( % Write the largest integer square root into that position
 tt % Duplicate twice
 5M % Push again the matrix which is true for the position of maximum
 1Y6 % Push matrix [0 1 0; 1 0 1; 0 1 0] (von Neumann neighbourhood)
 Z+ % 2D convolution, keeping size. Gives a matrix which is 1 for the
 % neighbours of the value that was replaced by its square root
 * % Multiply. This replaces the value 1 by the actual values of
 % the neighbours
 t % Duplicate
 XzX< % Minimum of non-zero entries
 = % Equals, element-wise. This gives a matrix which is true at the
 % position of the maximum neighbour, and zero otherwise
 * % Multiply, element-wise. This gives a matrix which contains the
 % maximum neighbour, and has all other entries equal to zero
 J % Push the maximum integer root, which was previously stored
 q % Subtract 1
 * % Multiply element-wise. This gives a matrix which contains the
 % maximum neighbour times (maximum integer root minus 1)
 + % Add. This replaces the maximum neighbour by the desired value,
 % that is, the previously found maximum integer square root
 % times the neighbour value
 w % Swap
 } % Else. This means there was no integer square root, so no more
 % iterations are neeeded
 ** % Multiply element-wise twice. Right before this the top of the
 % stack contains a zero. Below there are the latest matrix with
 % square roots and two copies of the latest matrix of integers,
 % one of which needs to be displayed as final result. The two
 % multiplications leave the stack containing a matrix of zeros
 % and the final result below
 % End (implicit). The top of the stack is consumed. It may be a
 % positive number, which is truthy, or a matrix of zeros, which is
 % falsy. If truthy a new iteration is run. If falsy the loop exits
 % Display (implicit)

JavaScript (ES6), 271 259 250 245 bytes

m=>for(l=m.length;I=J=Q=-1;)0)[y]

Thanks to Luis felipe De jesus Munoz for −14 bytes!

Explanation:

m => // m = input matrix
 // l = side length of square matrix
 // I, J = i, j of largest square in matrix (initialized to -1 every iteration)
 // Q = square root of largest square in matrix
 for (l = m.length; (I = J = Q = -1); ) 1 / 0);
 // x = i, y = j of smallest adjacent neighbor of largest square
 [x, y] = d[D.indexOf(Math.min(...D))];
 // multiply smallest adjacent neighbor by square root of largest square
 m[x][y] *= Q;
 // set largest square to its square root
 m[I][J] = Q;
 // repeat until no remaining squares in matrix
 // no return necessary; input matrix is modified.
;

C# (Visual C# Interactive Compiler), 220 bytes

n=>l=>for(int g;n.Any(x=>Math.Sqrt(x)%1==0);n[n.Select((a,b)=>(x:Math.Abs(b/l-g/l)+Math.Abs(b%l-g%l)==1?a:1<<30,y:b)).OrderBy(x=>x).First().y]*=n[g])n[g=n.IndexOf(n.Max(x=>Math.Sqrt(x)%1==0?x:0))]=(int)Math.Sqrt(n[g]);

Try it online!

Wolfram Language (Mathematica), 224 bytes

(l=#;While[(c=Length)[m=Select[Join@@l,IntegerQ[Sqrt@#]&]]>0,t=##&@@#&@@SortBy[Select[(g=#&@@Position[l,f=Max@m])+#&/@1,0,0,1,-1,0,0,-1,Min@#>0&&Max@#<=c@l&],l[[##]]&@@#&];l[[##&@@g]]=(n=Sqrt@f);l[[t]]=l[[t]]n];l)&

Try it online!

JavaScript (Node.js), 157 bytes

a=>g=(l,m=n=i=j=0)=>a.map((o,k)=>m>o||o**.5%1||[m=o,i=k])|m&&a.map((o,k)=>n*n<o*n|((i/l|0)-(k/l|0))**2+(i%l-k%l)**2-1||[n=o,j=k])|[a[i]=m**=.5,a[j]=m*n]|g(l)

Try it online!

-14 bytes thanks the @Arnauld who also wrote a nice test harness :)

Anonymous function that takes a 1-dimensional array as input and a length parameter specifying number if columns/rows.

Curried input is specified as f(array)(length).

// a: 1-dimensional array of values
// g: recursive function that explodes once per recursive call
// l: number of columns, user specified
// m: max square value
// n: min neighbor
// i: index of max square
// j: index of min neighbor
a=>g=(l,m=n=i=j=0)=>
 // use .map() to iterate and find largest square
 a.map((o,k)=>
 // check size of element
 m>o||
 // check if element is a square
 o**.5%1||
 // new max square found, update local variables
 [m=o,i=k])|
 // after first .map() is complete, continue iff a square is found
 // run .map() again to find smallest neighbor
 m&&a.map((o,k)=>
 // check size of element
 n*n<o*n|
 // check relative position of element
 ((i/l|0)-(k/l|0))**2+(i%l-k%l)**2-1||
 // a new smallest neighbor found, update local variables
 [n=o,j=k])|
 // update matrix in-place, largest square is reduced,
 // smallest neighbor is increased
 [a[i]=m**=.5,a[j]=m*n]|
 // make recursive call to explode again
 g(l)

Java 8, 299 297 bytes

m->for(int l=m.length,i,j,I,J,d,M,t,x,y;;m[x][y]*=d)for(i=l,I=J=d=0;i-->0;)for(j=l;j-->0;d=t>d*d&Math.sqrt(t)%1==0?(int)Math.sqrt(m[I=i][J=j]):d)t=m[i][j];if(d<1)break;for(M=-1>>>1,m[x=I][y=J]=d,t=4;t-->0;)tryM=m[i=t>2?I-1:t>1?I+1:I][j=t<1?J-1:t<2?J+1:J]<M?m[x=i][y=j]:M;catch(Exception e)

Modifies the input-matrix instead of returning a new one to save bytes.

Try it online.

Explanation:

m-> // Method with integer-matrix input and no return-type
 for(int l=m.length, // Dimension-length `l` of the matrix
 i,j,I,J,d,M,t,x,y; // Temp integers
 ; // Loop indefinitely:
 m[x][y]*=d) // After every iteration: multiply `x,y`'s value with `d`
 for(I=J=d=0, // (Re)set `I`, `J`, and `d` all to 0
 i=l;i-->0;) // Loop `i` in the range (`l`, 0]:
 for(j=l;j-->0; // Inner loop `j` in the range (`l`, 0]:
 d= // After every iteration: set `d` to:
 t>d*d // If `t` is larger than `d` squared
 &Math.sqrt(t)%1==0?
 // And `t` is a perfect square:
 (int)Math.sqrt(m[I=i][J=j])
 // Set `I,J` to the current `i,j`
 // And `d` to the square-root of `t`
 :d) // Else: leave `d` the same
 t=m[i][j]; // Set `t` to the value of `i,j`
 if(d<1) // If `d` is still 0 after the nested loop
 // (which means there are no more square-numbers)
 break; // Stop the infinite loop
 for(M=-1>>>1, // (Re)set `M` to Integer.MAX_VALUE
 m[x=I][y=J]=d, // Replace the value at `I,J` with `d`
 t=4;t-->0;) // Loop `t` in the range (4, 0]:
 tryM= // Set `M` to:
 m[i=t>2? // If `t` is 3:
 I-1 // Go to the row above
 :t>1? // Else-if `t` is 2:
 I+1 // Go to the row below
 : // Else (`t` is 0 or 1):
 I] // Stay in the current row
 // (and save this row in `i`)
 [j=t<1? // If `t` is 0:
 J-1 // Go to the column left
 :t<2? // Else-if `t` is 1:
 J+1 // Go to the column right
 : // Else (`t` is 2 or 3):
 J] // Stay in the current column
 // (and save this column in `j`)
 <M? // And if the value in this cell is smaller than `M`:
 m[x=i][y=j] // Set `x,y` to `i,j`
 // And `M` to the current value in `i,j`
 :M; // Else: leave `M` the same
 catch(Exception e) // Catch and ignore IndexOutOfBoundsExceptions

Jelly, 70 67 bytes

’dL½$}©+Ø.,U$;N$¤%®‘¤<®Ạ$Ƈæ.®‘ịÐṂḢ;ḷ;€ị½¥×÷ƭ⁹Ḣ¤¦Ṫ}¥ƒ
×Ʋ$MḢçɗ⁸Ẹ?ƊÐL

Try it online!

I’m sure this can be done much more briefly, but I found this harder than it first appeared. Explanation to follow once I’ve tried to golf it better.

A full program that takes a list of integers corresponding to the square matrix and returns a list of integers representing the final exploded matrix.l