Does it really satisfy detailed balance? If I am understanding the algorithm correctly, it seems like the probability of an L-tetromino transforming to square-tetromino is twice as large as a square-tetromino transforming to L-tetromino.

I think the only problem with the quoted construction is that the algorithm sometimes gets stuck. But as someone else pointed out, this does not mean that it is not sampling correctly––you just have to start over whenever you get stuck. All one needs to check is that every $n$-omino gets counted exactly $n$ times, which seems straightforward unless I am missing something. Still, this does not answer your questions about the expected number of $2 \times 2$ squares as $ n \to \infty$.

Brendan McKay, yes this my intention, so this would be useful. Thank you.

Igor Pak. Thanks for the reference. I have been looking at this book, but still haven't found the answer to my question. Maybe it is right in front of me, and I am missing it. As far as I understand it, self-avoiding random walks can be "nearly uniformly" sampled via relatively easy to implement Markov chain Monte Carlo methods such as the pivot algorithm, for example: arxiv.org/abs/cond-mat/0109308 . I guess I am asking if there are similar methods for simulating uniform polyominoes. I am guessing that proving things about mixing time is hard, and I am happy to just do experiments.

Aaron Meyerowitz, it seems that this keeps you from getting stuck, but then not every polyomino is equally likely. For example, the square tetromino is more likely than the straight tetromino.