# [R][1], <s>301</s> <s>287</s> <s>277</s> <s>274</s> <s>222</s> <s>217</s> 195 bytes

Nothing particularly creative, including the zero buffering of the peripheral elements of the entry matrix:
 
 p=function(x){
 w=which.max;d=cbind
 n=sum(x|1)^.5
 s=which(!x^.5%%1)
 if (sum(s|1)){
 y=d(0,t(d(0,x,0)),0)
 i=w(y*y%in%x[s])
 y[!y]=max(x)
 n=-c(1,n+2)
 m=i+c(n,-n)
 j=m[w(-y[m])]
 z=y[i]=y[i]^.5;y[j]=y[j]*z
 x=p(y[n,][,n])}
 x}

[Try it on-line][2]

Nick Kennedy managed [a better version][3] of the algorithm as follows:

 w=which.max;`~`=cbind;x=scan();while(any(s<-!x^.5%%1)){y=(k=Inf)~t(k~matrix(x,n<-sum(x|1)^.5)~k)~k;i=w(y*y%in%x[s]);m=i+c(-1,-(n=n+2),1,n);y[j]=y[j<-m[w(-y[m])]]*(y[i]=y[i]^.5);x=y[q<-3:n-1,q]};x

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

[Try it on-line][3]


 [1]: https://www.r-project.org/
 [2]: https://tio.run/##JY3RasMwDEXf/RXZQ0BulFK3adno9CXGg86h1AFroUmJna3fninbwxWcy@HqvvS0XB/sx/DFkPS3mmi6BX/bxks6t@Q/A7eKaXhESD9Gf2yPavg34CUJlaXRKlxhFQYRZKHI1ML8XsuCLFYGRxDGhLPWElUEmiBvchm4THZw0jDVHgxytReIFCoPjDULdBTtBHW20WmnipmyDe7vyPNztt0KndvMqkjUQ7aMziI7/VTpufQQL@M9JPBwOKHZ7dBg84amaXB/PElWejUaD1ovvw
 [3]: https://tio.run/##FU5basMwELyK@mHYTawQJXZIKu8BegahEtd5WHG0kNhgiVJf3ZVhZ2AeDPue55HG1jXtxtdBn6czNT@OLzpQ39QMqFP4vELNEfpKfoTvTZllCvE3EnT0xTecBugmXw9vFyDkXMnnle9DCwFTF6cunXY0QlzFzHEWTG9Re3LrBqTKJTDxeoe5yhl1NA9LiSrpzQgyGm/R2hVE4xbf2WUy/RbNq5L7T04DL/unw7w/CLXdCiWKk1BFIXblIWFRRzX/Aw