Given a (0,1)-matrix $A$, I'll denote by $\mu(A)$ the number of maximal monochromatic polyominoes in $A$ (i.e., the number of connected polyominoes contained in $A$ each of which is either all 0 or all 1 such that each polyomino is as large as possible – so they thus tile $A$). I do not know of any literature on $\mu(A)$, so I apologize if there is already some notation for the quantity that I am not using.
So for example if $A$ is the $n \times n$ identity matrix, $\mu(A) = n + 2$. For an $n \times n$ matrix $A$, we have $\mu(A) \leq n^2$, with equality if and only if $A$ is one of the "checkerboard matrices".
I would like to understand what $\mu(A)$ is on average, as we allow the size of the matrix $A$ to tend toward infinity. In particular, is it true that
$$
\lim_{n \to \infty} \frac{\text{Average value of $\mu$ over all $n\times n$ (0,1)-matrices}}{n^2} = 1
$$
Best Answer
The limit is positive but well below $1$; I show that the limit is in $(0.09448, 0.2646)$, and outline how to approximate it much more closely than this.
The average of $\mu$, call it ${\rm E}(\mu)$, is twice the average number of maximal all-0 polyominos; so the limit of ${\rm E}(\mu) / n^2$ is twice the expected number of maximal all-0 polyominos per unit area of a large square. That limit is the special case $p=q=1/2$ of the expected number, call it $\epsilon(p)$, of maximal all-0 polyominos per unit area when each entry is independently 0 or 1 with probability $p,q$ respectively ($p+q = 1$). We can approximate $\epsilon(p)$ by writing it as a sum $\sum_P \epsilon_P(p)$ over all possible polyomino shapes $P$, and calculating partial sums over all $P$ of size at most $k$.
The contribution $\epsilon_P(p)$ is $p^{|P|} q^{|\delta(P)|}$, where $|P|$ is the size of $P$, and $|\delta(P)|$ is the number of cells at distance $1$ from $P$. (We do not identify different orientations; for example, for $k=4$ there is one square $P$, two straight ones, four S/Z shaped tetrominos, four T-shaped, and eight that are L-shaped.) So we seek $\sum_P p^{|P|} q^{|\delta(P)|}$. Note that $\sum_P |P| p^{|P|} q^{|\delta(P)|} = p$ because this is just the expected number of 0 entries per unit area.
The analysis so far works in any dimension; for example, in dimension $1$ there is a unique $P$ of each size $k \geq 1$, and $|\delta(P)| = 2$ always, so $\epsilon(p) = \sum_{k=1}^\infty p^k q^2 = p q^2/(1-p) = p q$ (and indeed $\sum_{k=1}^\infty k p^k q^2 = p q^2/(1-p)^2 = p$). For example, $\epsilon(1/2) = 1/4$, and the expected number of maximal all-0 or all-1 strings in a random bitstring of length $n$ is asymptotically $2\epsilon(1/2) n = n/2$. (Exercise [noted in my comment]: in fact that expected number is exactly $(n+1)/2$. What happens for arbitrary $p$?)
In dimension $2$, the sum of $\epsilon_P(p)$ over $|P| \leq 4$ is $$ p q^4 + 2 p^2 q^6 + p^3 (2 q^8 + 4 q^7) + p^4 (q^8 + 2 q^{10} + 4 q^8 + 4 q^8 + 8 q^9), $$ and the sum of $|P| \epsilon_p(P)$ up to $|P| \leq 4$ is obtained from by multiplying the $k$-th term by $k$ ($k=1,2,3,4$), and comes to $p - 315 p^5 + 1824 p^6 - 4848 p^7 + \cdots$. (The coefficient $315$ seems to come from the number of rooted 5-ominos, see https://oeis.org/A048664 .) For $p=1/2$ these sums are only $387/2^{13}$ and $612/2^{13}$ respectively. It follows that the $|P| \geq 5$ terms must contribute $1/2 - 612/2^{13}$ to the sum of $|P| \epsilon_p(P)$, and thus can contribute at most $1/5$ of that to the sum of $\epsilon_p(P)$. We conclude that $$ \frac{387}{2^{13}} < \epsilon(1/2) < \frac{387}{2^{13}} + \frac15 \left( \frac12 - \frac{612}{2^{13}} \right) = \frac{5419}{40960}; $$ numerically these bounds are $0.04724+$ and just under $0.1323$. Hence the answer to the OP's question is between $0.09448$ and $0.2646$, as claimed.
Polyominos have been tabulated well beyond $k=4$. It must be feasible to compute $\sum_{|P| \leq k} \epsilon_P(1/2)$ and $\sum_{|P| \leq k} |P| \epsilon_P(1/2)$ for $k$ large enough to obtain a reasonably close approximation to $\epsilon(1/2)$, and thus to the value $2\epsilon(1/2)$ of the limit that the OP seeks.