sed 4.2.2 _{^{-r, 1303 1252 1201 1053 838 764 710 681 629 615}} 536 bytes

s/^/0/                                               #set initial array to test against
:;h;s/\W+//;y/10/ab/;G                               #the start of the main loop
:b;s/(.*)\b([ab]+)( *)/\2\n\3\1/;h;s/\n.*/\n.... ./  #the start of the shuffle loop
:c;s%(.*)\n((.*) (.*))%sed -r 's/^(\3|.)$/\4\4\&/;T;s/\4/\&\\n/2g;q;:;a\3\3\3\2\4'<<<\1%e;/\n$/!bc
s%.*%sed '1h;1d;x;2!s/\\.//;T;/^\\n/{s/$/\\n/p;z};:;$!x;3~2{H;d}'<<<'&'%e
s/(.)\n/\1/g;y/\n/ /;G;s/\n\w*\n/ /;/[ab]{2}/bb      #the end of the shuffle loop
h;s/\n.*//;y/ /\n/                                   #formatting for the calc loop
:d;s!.*!sed -n 'h;:;x;p;x;s/.//;//b'<<<'&'!e         #start of the calc loop; square everything
:e;s/\b\n\b//;s/ba|ab/a/;te                          #add those squares
s/\n(\n*)/\1/g;/\n/bd                                #remove a layer, end calc loop
y/a/#/;G;s/\n.*\n//                                  #grab input to compare against
T;s/^b*(#*)(\w*)\1$/\2/;tz                           #if they match, output the matrix
s/b*#*//;T                                           #delete the calc result
s/10(0*)(1*)#/01\2\1#/;t                             # v v v v v v v v v v v v v v v v v
s/(0*)0(1*)/\21\1/;t                                 #get the next matrix, end main loop
s/1/0000/g;t                                         # ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^
:z                                                   #label for escaping the infinite loop

Try it online!

did you know sed is turing complete? I happen to know that, and it makes me do things I shouldn't...

this does basically the same thing that the other answers are doing. we come up with a way to generate an infinite list of matrices (check the order here), calculate the value of each matrix (check the math here), and then see if it's equal to the input. These two individual steps are fast and easy to check edge cases with; it's only the combined loop that takes a long time.

the tio link has cases 0,1,2,3,4,5,9 loaded in, which finishes in ~30 seconds, and 10 will also finish in just under 60 seconds; but if you want to try any other numbers you'll have to run it on your own computer.

in the explanation next to the code, I reference a 'shuffle loop' and a 'calc loop'. both of these are part of the above 'check the math' link, and together they find the value of an arbitrary matrix. the shuffle loop works by first inserting the correct number of newlines to make a 4n*n matrix, and then re-ordering the lines so that we have 4 separate quadrants. this gets recursively applied until all arrays are one digit long. you can think of this as like making a list of lists of lists of etc until each item is only one bit. the calc loop then takes this list, squares each number, replaces the innermost lists with their sums, and repeats.

Python, 615 bytes

from numpy import*
Z=zeros
r=range
R=lambda n:r(n,-1,-1)
N=lambda M:M.shape[0]
z=lambda M,n:M if N(M)==n else(lambda w:(a:=Z((n,n)),[a.__setitem__((i*w,j*w),M[i,j])for i in r(N(M))for j in r(N(M))],a)[2])([2**i for i in r(n.bit_length())if N(M)*2**i==n][0])
def g(n):
    if n<2:return array([[n]])
    if n<5:r=Z(4);r[:n]=1;return r.reshape(2,2)
    else:a,b,c,d=(lambda n:next((i,j,k,l)for i in R(n)for j in R(i)for k in R(j)for l in R(k)if i*i+j*j+k*k+l*l==n))(n);A,B,C,D=map(g,(a,b,c,d));m=max(map(N,[A,B,C,D]));return block([[z(A,m),z(B,m)],[z(C,m),z(D,m)]])
lambda n:";".join("".join(map(str,map(int,l)))for l in g(n))

Test

Assign the function to G, and

>>> for i in (0, 1, 3, 4, 5, 9, 15, 67, 78):
>>>     print(i, G(i), sep=" -> ")
0 -> 0
1 -> 1
3 -> 11;10
4 -> 11;11
5 -> 1110;0000;0000;0000
9 -> 1100;1000;0000;0000
15 -> 1111;1000;1010;0000
67 -> 11111000;00000000;00000000;00000000;10001000;00000000;00000000;00000000
78 -> 11111010;00000000;00001000;00000000;10101000;00000000;00000000;00000000

05AB1E, 50 49 bytes

∞Íonîü2íε`U0LsãεXô}Σ˜O]í€`.Δ[D˜g5‹#2δô2ô€øOOn}OOQ

Brute-force. Generates all \$2^k\times 2^k\$ sized matrices consisting of 0s and 1s, and will output the first valid one for the given input-integer \$n\$.

-1 byte by changing the ∞<on1š (-1 and prepend \$1\$ manually) to ∞Íonî (-2 and ceil the \$0.25\$ to \$1\$)

Try it online or verify the smaller test cases up to n=7.

Explanation:

Step 1: Generate a list of all possible outputs:

∞          # Push an infinite positive list: [1,2,3,4,5,...]
 Í         # Decrease each by 2 to start at -1: [-1,0,1,2,3,...]
  o        # Take 2 to the power each: [0.5,1,2,4,8,...]
   n       # Square each: [0.25,1,4,16,64,...]
    î      # Ceil each: [1,1,4,16,64,...]
ü2         # Create overlapping pairs: [[1,1],[1,4],[4,16],[16,64],...]
  í        # Reverse each pair: [[1,1],[4,1],[16,4],[64,16],...]
ε          # Map over each pair:
 `         #  Pop and push both to the stack
  U        #  Pop the top smaller value, and save it in variable `X`
  0L       #  Push pair [1,0]
    s      #  Swap so the larger value is at the top of the stack
     ã     #  Cartesian product to create all combinations using the 1s/0s
      ε    #  Map over each list of 1s/0s
       Xô  #   Split it into blocks of size `X`
      }    #  Close the inner map
 Σ         #  Sort by:
  ˜        #   The flattened
   O       #   sum (of 1s)
]          # Close the sort and outer map
 í         # Reverse each inner list
  €`       # Flatten the list of list of lists one level down
           # (which has as side-effect that it reverses the inner lists,
           # hence the need for `í`)

Try just this first step online.

Step 2: Find the first matrix that's a valid reduction for the given input-integer:

.Δ         # Find the first matrix in this infinite list that's truthy for:
  [        #  Start an infinite list:
   D       #   Duplicate the current matrix
    ˜      #   Flatten it to a list
     g     #   Pop and push its length
      5‹   #   If this length is smaller than 5 (e.g. either 4 or 1):
        #  #    Stop the infinite inner loop
   2δô2ä€ø #  Split the matrix into four corner-blocks:
    δ      #   Map over each inner row:
   2 ô     #    Split it into pairs
      2ä   #   Then split the entire list into 2 equal parts as well
        €  #   Map over each inner matrix:
         ø #    Zip/transpose; swapping rows/columns
   OO      #  Sum each inner-most matrix
     n     #  Then square those sums
  }        # After the infinite loop to reduce it to a 2x2 matrix:
   OO      # Sum the four values in this matrix as well
     Q     #  Check whether it now equals the (implicit) input-integer
           # (after which the found matrix is output implicitly)

Verify the resulting \$n\$ for a given matrix to test just this second step (minus leading .Δ and trailing Q). I've added = after both OO, so you can also see the intermediate steps of the reduction.

Charcoal, 101 bytes

Ｎθ≔⦃⁰⟦0⟧¹⟦1⟧⦄ηＷ∧¬§ηθη«≔⦃⦄η≔ＥιλζＦζＦζＦζＦζ«≔ΣＸ⟦νμλκ⟧²ε≔⁺Ｅ§ιν⁺ιμπＥ§ιλ⁺ικπδ¿∨¬§ηε‹ΣΣδΣΣ§ηε§≔ηεδ»»§ηθ

Attempt This Online! Link is to verbose version of code. Explanation: Lots of brute force, so can only solve inputs with 4×4 outputs before the heat death of the universe.

Ｎθ

Input n.

≔⦃⁰⟦0⟧¹⟦1⟧⦄η

Start with solutions for 0 and 1.

Ｗ∧¬§ηθη«

Repeat until a solution for n is found. (Also makes a copy of the solutions so far.)

≔⦃⦄η

Start building solutions of twice (4×) the size.

≔ＥιλζＦζＦζＦζＦζ«

For all possible sets of subsquares:

≔ΣＸ⟦νμλκ⟧²ε

Calculate the value.

≔⁺Ｅ§ιν⁺ιμπＥ§ιλ⁺ικπδ

Calculate the resulting square.

¿∨¬§ηε‹ΣΣδΣΣ§ηε§

If this is a better square, then...

≔ηεδ

... add or replace the previous solution for that value.

»»§ηθ

Output the square for n.

108 bytes for a slightly less forceful version:

Ｎθ≔⦃⁰⟦0⟧¹⟦1⟧⦄ηＷ∧¬§ηθη«≔⦃⦄η≔ΦＥι묛×κκθζＦζＦζＦζＦζ«≔ΣＸ⟦νμλκ⟧²ε≔⁺Ｅ§ιν⁺ιμπＥ§ιλ⁺ικπδ¿∨¬§ηε‹ΣΣδΣΣ§ηε§≔ηεδ»»§ηθ

Attempt This Online! Link is to verbose version of code. Explanation: As above but when generating the list of possible solutions, excludes values whose square exceeds the input.

Pyth, 48 bytes

M?Gho_SsNmsCMc2sgLtGdfqHs*VTT^UhH4]]Hgf>2eSsgTQZ

Try it online!

Can run up to about 30 before timing out online. Try a more performant, 5 byte longer version here (can easily do 1000 without timing out). The faster version is the same except instead of looking at all permutations of 4 numbers in the entire range of the input, we only look at all combinations (with repetition) of 4 numbers in the range up to the square root of the input.

This answer uses a recursive approach.

Explanation

M?Gho_SsNmsCMc2sgLtGdfqHs*VTT^UhH4]]Hgf>2eSsgTQZQ    # implicitly add Q
                                                     # implicitly assign Q = eval(input())
M                                                    # define g(G,h)
 ?G                               ]]H                #   recursive base case, if G==0: return [[H]]
                              UhH                    #   list(range(H+1))
                             ^   4                   #   4x repeated cartesian product
                     f                               #   filter this list on lambda T (this keeps only quadruplets whose squares sum to the desired value)
                        s                            #     sum of
                         *VTT                        #     vectorized product of T and T
                      qH                             #     is equal to H
         m                                           #   map over lambda d (this constructs possible outputs from the recursive calls)
                gLtGd                                #     map g over d with a first argument of G-1
          sCMc2s                                     #     take four square matrices and put them in a bigger square
    o                                                #   sort by lambda N
       sN                                            #     flattened N
      S                                              #     sorted
     _                                               #     reversed
   h                                                 #   take the first element
                                      f        Z     # find the first integer which satisfies lambda T, counting up from 0
                                            gTQ      #   g(T, Q)
                                           s         #   flatten list
                                         eS          #   largest element
                                       >2            #   is less than 2
                                     g          Q    # apply g on the result and the input