Make all the squares explode You are given a square matrix of width \$\ge2\$, containing square numbers \$\ge1\$.
Your task is to make all square numbers 'explode' until all of them have disappeared. You must print or return the final matrix.
More specifically:

Look for the highest square \$x^2\$ in the matrix.
Look for its smallest adjacent neighbor \$n\$ (either horizontally or vertically and without wrapping around).
Replace \$x^2\$ with \$x\$ and replace \$n\$ with \$n\times x\$.

Repeat the process from step 1 until there's no square anymore in the matrix.
Example
Input matrix:
$$\begin{pmatrix}
625 &amp; 36\\
196 &amp; 324
\end{pmatrix}$$
The highest square \$625\$ explodes into two parts of \$\sqrt{625}=25\$ and merges with its smallest neighbor \$36\$, which becomes \$36\times 25=900\$:
$$\begin{pmatrix}
25 &amp; 900\\
196 &amp; 324
\end{pmatrix}$$
The highest square \$900\$ explodes and merges with its smallest neighbor \$25\$:
$$\begin{pmatrix}
750 &amp; 30\\
196 &amp; 324
\end{pmatrix}$$
The highest square \$324\$ explodes and merges with its smallest neighbor \$30\$:
$$\begin{pmatrix}
750 &amp; 540\\
196 &amp; 18
\end{pmatrix}$$
The only remaining square \$196\$ explodes and merges with its smallest neighbor \$18\$:
$$\begin{pmatrix}
750 &amp; 540\\
14 &amp; 252
\end{pmatrix}$$
There's no square anymore, so we're done.
Rules

The input matrix is guaranteed to have the following properties:

at each step, the highest square will always be unique
at each step, the smallest neighbor of the highest square will always be unique
the sequence will not repeat forever


The initial matrix may contain \$1\$'s, but you do not have to worry about making \$1\$ explode, as it will never be the highest or the only remaining square.
I/O can be processed in any reasonable format
This is code-golf

Test cases
Input : [[16,9],[4,25]]
Output: [[24,6],[20,5]]

Input : [[9,4],[1,25]]
Output: [[3,12],[5,5]]

Input : [[625,36],[196,324]]
Output: [[750,540],[14,252]]

Input : [[1,9,49],[1,4,1],[36,25,1]]
Output: [[3,6,7],[6,2,7],[6,5,5]]

Input : [[81,4,64],[16,361,64],[169,289,400]]
Output: [[3,5472,8],[624,323,1280],[13,17,20]]

Input : [[36,100,1],[49,144,256],[25,49,81]]
Output: [[6,80,2],[42,120,192],[175,21,189]]

Input : [[256,169,9,225],[36,121,144,81],[9,121,9,36],[400,361,100,9]]
Output: [[384,13,135,15],[24,1573,108,54],[180,11,108,6],[380,209,10,90]]

Input : [[9,361,784,144,484],[121,441,625,49,25],[256,100,36,81,529],[49,4,64,324,16],[25,1,841,196,9]]
Output: [[171,19,700,4032,22],[11,210,525,7,550],[176,60,6,63,23],[140,112,1152,162,368],[5,29,29,14,126]]
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&lt;-!x^.5%%1)){
y=cbind(NA,rbind(NA,x,NA),NA)
z=y[i]=y[i&lt;-w(y*y%in%x[s])]^.5
m=i+c(r&lt;--c(1,nrow(y)),-r)
y[j]=y[j&lt;-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&lt;-!x^.5%%1)){y=NA~t(NA~matrix(x,n&lt;-length(x)^.5)~NA)~NA;i=w(y*y%in%x[s]);=i+c(r&lt;--c(1,n+2),-r);y[j]=y[j&lt;-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