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
$$



Average number of tiles of a (0,1)-matrix?
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.


