Uiua, 45 bytes

°√/↥♭∧(⍜⊡◌⊙::×+1/↧⊙◡⊡⊸𝄐⌟⬚0⊡⊸≡⌟-↘1♭₂⇡2_2)♭₂⇡△.

Pad

Uses a dynamic programming algorithm to achieve true \$O(mn)\$ time: let \$A_{ij}\$ be the input matrix, and let \$S_{ij}\$ be the largest square side length including \$(i, j)\$ between \$(0, 0)\$ and \$(i, j)\$.

If cell \$A_{ij}=1\$, then \$S_{ij}=1+\min(S_{(i-1)j}, S_{i(j-1)}, S_{(i-1)(j-1)})\$, because \$S_{ij}=n\ \iff\$ there are squares above, left, and diagonally (above-left) with \$S=n-1\$ (the union of these squares plus cell \$(i, j)\$ creates the new square). Thus, we set \$S_{ij}=A_{ij}*(1+\min(S_{(i-1)j}, S_{i(j-1)}, S_{(i-1)(j-1)}))\$ (if \$A_{ij}=0\$ then \$S_{ij}=0\$.

The code fills a table of \$S\$ lexicographically on input coordinates (as each cell depends only on previous coordinates in this ordering). It stores \$S\$ in the same array as \$A\$, since once we set \$S_{ij}\$ we no longer care about \$A_{ij}\$, and this simplifies stack manipulation significantly.

F ← (
  ♭₂⇡△.      # create array of coords
  ∧(         # for each coord (in ascending order):
    ↘1♭₂⇡2_2 # the array [0_1 1_0 1_1]
    ⊸≡⌟-     # subtract from current coord
    ⊸𝄐⌟⬚0⊡   # lookup squares from table
    ⊙◡⊡      # lookup current coord 1/0
    ×+1/↧    # (min adj + 1)*current
    ⍜⊡◌⊙::   # set in table
  )
  °√/↥♭      # square the last elem in table
)

JavaScript (Node.js), 136 bytes

U=>U.map(L=(r,y)=>r.map((c,x)=>k=k<(t=L[[,y,x]]=Math.min(-~L[[,y-1,x-1]],U[[y,x]]=c*-~U[[y-1,x]],L[[y,x]]=c*-~L[[y,x-1]]))?t:k),k=0)|k*k

Try it online!

Assumes "2D array" L[[y,x]] is O(1)

C, C++ : 231 227 196 bytes

-4 bytes thanks to EzioMercer

-31 bytes thanks to ceilingcat

int f(int t[],int x,int y){int z=1,j=y,i,v,w,m,k,b;for(;j--;)for(i=x;i--;)if(t[i*y+j]){v=k=1;w=0;for(m=x-i<y-j?x-i:y-j;k++<m;w=v?k:w)for(b=k*k;b--;)v=t[y*(i+b%k)+j+b/k]?v:0;z=w>z?w:z;}return z*z;}

TIO Link thanks The Thonnu for pointing it out

Test code :

int t[5][4] = {
    {1,1,1,1},
    {0,0,1,0},
    {1,1,1,0},
    {1,1,1,1},
    {1,1,1,0}
};

int main() {
    printf("r=%d\n", f((int*)t,5,4));
    return 0;
}

Charcoal, 88 80 46 bytes

≔Ｅθ⁰θＷＳ«≔⮌θη≔⟦⁰⟧θＦι⊞θ×Ｉκ⊕⌊⟦↨υ⁰⊟η↨η⁰⟧⊞υ⌈θ»ＩＸ⌈υ²

Try it online! Link to verbose version of code. Takes input as a list of newline-terminated bit strings. Explanation: Now a port of @Jonah's J solution.

≔Ｅθ⁰θ

Start with a row of zeros.

ＷＳ«

Loop over each input row.

≔⮌θη

Reverse the previous row totals to make them easier to process.

≔⟦⁰⟧θ

Start collecting a new row with a zero column.

Ｆι

For each cell...

⊞θ×Ｉκ⊕⌊⟦↨υ⁰⊟η↨η⁰⟧

... multiply its value with the incremented minimum of the cells above and to the left. The cell above is obtained by popping the reversed previous row, the cell to the above left is obtained as the base 0 of the popped reversed previous row and the cell to the left is obtained as the base 0 of the new row.

⊞υ⌈θ

Collect the maximum square found in this row.

»ＩＸ⌈υ²

Output the largest found square from all rows.

Wolfram Language (Mathematica), 109 bytes

Catch[Do[If[AnyTrue[Erosion[ArrayPad[#,1],BoxMatrix[All,n]],#==1&,2],Throw[n]],{n,Min@Dimensions@#,1,-1}]]^2&

Try it online!

Vyxal, 20 bytes

ɖλZɖλ÷nfṪg›*";vt;fG²

Try it Online!

Based on @Jonah's J answer.

For each 2x2 square

a b
c d

it replaces d with d*(min(a,b,c)+1) in sequence going from left to right, top to bottom. Then it takes maximum and squares it.

ɖλ                   # scan by:
  Z                  #   zip
   ɖλ                #   scan by:  receives [a,c], [b,d]
     ÷               #     push each to stack
      n              #     push the function argument [[a,c],[b,d]]
       f             #     flatten
        Ṫ            #     remove tail  [a,c,b]
         g           #     minimum
          ›          #     increment
           *         #     multiply
            "        #     pair [b, new d]
             ;       #   end scan
              vt     #   get the tail of each
                ;    # end scan
                 f   # flatten
                  G  # maximum
                   ² # square

JavaScript, 186 bytes

I know that the complexity should be O(mn) but here the complexity is O((mn)^2) :)

m=>eval('for(x=b=0;x<(r=m.length);++x)for(y=0;y<(c=m[0].length);++y)if(m[x][y]){for(o=s=1;o&&s<(c>r?r:c);++s)for(t=s+1;t--;)if(!m[x+s]?.[y+t]||!m[x+t][y+s]){o=0;--s;break}b=s>b?s:b}b*b')

Try it:

f=m=>eval('for(x=b=0;x<(r=m.length);++x)for(y=0;y<(c=m[0].length);++y)if(m[x][y]){for(o=s=1;o&&s<(c>r?r:c);++s)for(t=s+1;t--;)if(!m[x+s]?.[y+t]||!m[x+t][y+s]){o=0;--s;break}b=s>b?s:b}b*b')

;[
    [
        [1, 0, 1, 1, 1],
        [1, 0, 1, 1, 1],
        [1, 1, 1, 1, 1],
        [1, 0, 0, 1, 0],
    ],
    [
        [1, 0, 1, 0, 1],
        [1, 0, 1, 1, 1],
        [1, 1, 1, 1, 1],
        [1, 0, 0, 1, 0],
    ],
    [
        [1, 0, 1, 1, 1],
        [1, 0, 1, 0, 1],
        [1, 1, 1, 1, 1],
        [1, 0, 0, 1, 0],
    ],
    [
        [1]
    ],
    [
        [0]
    ]
].map(m => console.log(f(m)))

JavaScript, 200 bytes

Here the complexity is much less than O((mn)^2) because if the algorithm find the biggest possible square it will stop

m=>eval('a=(c=m[0].length)>(r=m.length)?r:c;for(x=b=0;b!=a&&x<r;++x)for(y=0;b!=a&&y<c;++y)if(m[x][y]){for(o=s=1;o&&s<a;++s)for(t=s+1;t--;)if(!m[x+s]?.[y+t]||!m[x+t][y+s]){o=0;--s;break}b=s>b?s:b}b*b')

Try it:

f=m=>eval('a=(c=m[0].length)>(r=m.length)?r:c;for(x=b=0;b!=a&&x<r;++x)for(y=0;b!=a&&y<c;++y)if(m[x][y]){for(o=s=1;o&&s<a;++s)for(t=s+1;t--;)if(!m[x+s]?.[y+t]||!m[x+t][y+s]){o=0;--s;break}b=s>b?s:b}b*b')

;[
    [
        [1, 0, 1, 1, 1],
        [1, 0, 1, 1, 1],
        [1, 1, 1, 1, 1],
        [1, 0, 0, 1, 0],
    ],
    [
        [1, 0, 1, 0, 1],
        [1, 0, 1, 1, 1],
        [1, 1, 1, 1, 1],
        [1, 0, 0, 1, 0],
    ],
    [
        [1, 0, 1, 1, 1],
        [1, 0, 1, 0, 1],
        [1, 1, 1, 1, 1],
        [1, 0, 0, 1, 0],
    ],
    [
        [1]
    ],
    [
        [0]
    ]
].map(m => console.log(f(m)))

Explanation

I will not explain each line of code because it is a lot and boring

Here is illustration of algorithm:

When I meet the 1 in matrix (blue) I start checking the right and bottom edges of this cell (red cells) if all red cells are 1 then I start check the right and bottom edges of red cells (green cells) and so on until where I won't meet 0 or meet the matrix edge

J, ^{78 76} 75 bytes

[:>./@,@;($<@,:&2@#:i.@#@,)({:*1+<./@}:)@,;.0`(0<@{1+[)`]}&.>/@,&|.0<@,0&,.

Try it online!

\$\mathcal{O}(mn)\$ using a standard dynamic programming table.

The only challenge was how to translate this essentially procedural algorithm into an array paradigm:

• First add a top-left border of zeros:

   0 0 0 0 0 0
   0 1 0 1 1 1
   0 1 0 1 1 1
   0 1 1 1 1 1
   0 1 0 0 1 0
  

• Next generate the coordinates for what would be the nested loop in a procedural language -- that is, we traverse from left to right, top to bottom. We also add 2 2 to each, since we'll be cutting out 2x2 squares later. For example, the first element represents the 2x2 square whose top-left coordinate is 0 0:

  ┌───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┐
  │0 0│0 1│0 2│0 3│0 4│1 0│1 1│1 2│1 3│1 4│2 0│2 1│2 2│2 3│2 4│3 0│3 1│3 2│3 3│3 4│
  │2 2│2 2│2 2│2 2│2 2│2 2│2 2│2 2│2 2│2 2│2 2│2 2│2 2│2 2│2 2│2 2│2 2│2 2│2 2│2 2│
  └───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┘
  

• Now we reverse the above list and use it in a single reduction of the input from the previous step. In each step of the reduction, we slice out a two by two square. For example in the first step of the reduction we take the top-left square:

  0 0
  0 1
  

• The reduction logic is simple. Call the value of the lower-right cell lr. We update lr as follows: lr * (1 + min(all values except lr)).

• Return the max value of the completed dynamic programming table.

Excel, 142 bytes

=LET(b,SEQUENCE(ROWS(a)),c,TRANSPOSE(b),MAX(MMULT(N(COUNTIF(OFFSET(INDIRECT(TOCOL(CELL("address",OFFSET(a,b-1,c-1)))),,,c,c),"<>1")=0),b^0))^2)

Input a is a worksheet range comprising the matrix. The above should preferably be placed in a worksheet different to that containing the matrix so as to avoid potential circular references.

Complexity is unfortunately \$\mathcal{O}(m^3)\$.

Python NumPy, 148 bytes

lambda X,o=1:[B:=c(c(pad(X,1),0),1),[o:=o+(B[i,j]+B[i-o,j-o]-B[i,j-o]-B[i-o,j]==o*o)for i,j in argwhere(X)+1]]and~-o*~-o
from numpy import*
c=cumsum

Attempt This Online!

Should be linear in the size N (total number of elements) of the input.

How?

Computes first the 2d partial sums B of the input. After this one off O(N) (space and time) investment we can count the number of ones in any grid rectangle in constant time using the four corners of the rectangle: B(top left)+B(bottom right)-B(top right)-B(bottom left).

We can now find the maximum in linear time by traversing the input a single time top left to bottom right. We only check for squares that have their bottom right corner at the current position. If the current point is bottom right of a new maximum square we know that the current maximum must be one less because we have already visited the top left neighbour of the current point. Therefore we need at each point only check a single square, the loop body is constant time and the loop is linear.

MATL, 21 19 bytes

lYa`9I&ZItz}x@UGaa*

2 bytes saved thanks to @lmendo

The complexity of the erosion operation is O(mn), but unfortunately with this approach, we perform erosion as many as max(m, n) times.

Try it at MATL Online

Explanation

        % Implicitly grab input as a multi-dimensional array
lYa     % Pad with a row of zeros on the top and bottom to ensure erosion works
`       % do...while loop
  9I&ZI % Perform image erosion with a 3 x 3 neighborhood
  t     % Duplicate the output
  z     % Check if there are any 1's left in the eroded result, if so, repeat
}       % End loop
x       % Delete the last element on the stack (the eroded matrix of all 0's)
@       % Get the index of the last loop iteration and
U       % Square it to get the size of the largest contiguous square of 1's
Gaa     % Check if there are any 1's in the input
*       % Multiply with the result (to account for a matrix of all 0's)
        % Implicitly display the result