Simulation of Boltzmann distribution I'm tackling with the problem in Blundell_Concepts in Thermal Physics textbook. Let me introduce the context first. Below figure is just taken from Blundell textbook section 4.6.

As you can see in fig(a), energy quanta is distributed on a 20*20 lattice initially. Then I start to choose two sites randomly and make a quanta exchange between each site(fig (b)). Then, after many repetitions of such process, I get the final result in fig(c). Such distribution really resembles the Boltzmann distribution case.
So the question is following:

For the case of $N\gg1$ sites with initially one quantum per site, show that after many repetitions one could expect there to be $N(n)$sites with n quanta, where $N(n)\approx2^{-n}N$.

I know that this is just a form of Boltzmann distribution but how could I derive such formula? I guess that if I show that $N(n+1)=aN(n)$ holds(where a is constant), it's straightforward to show that $N(n)\approx2^{-n}N$ holds. But I don't know how to do it.(This problem is listed on exercise 4.9)
Um... and one more question. If I start the above situation not initially with one quantum per site but $Q\gg1$quantum per site, is it okay to say $N(n)\approx2^{-n}NQ$?
 A: Simple approach (does not account for temperature)
Let's focus on one particular site, say site $m$, where $0 \leq m < N$, and look at a single iteration of the quanta exchange process. The probability that $m$ loses a quantum in that iteration is $1/N$, because you could have chosen any of the $N$ sites. Likewise, the probability of $m$ gaining a quantum is $1/N$, because whichever site is chosen to lose a quantum, one of the $N$ sites has to be given it.
For $m$ to get $n$ quanta, it has to be chosen to gain a quantum $n$ times in a row without ever losing one. Since the probability of gain and loss is the same, the probability for this whole chain of events is $2^{-n}$. By multiplying the total number of sites with the probability to find a site with $n$ quanta, you get $N(n) = 2^{-n} N$ sites with $n$ quanta.
Remark 1: Actually, the probability of $m$ losing a quantum is higher or lower, because some sites will already have 0 quanta and cannot lose any more and $m$ might or might not be among them. I guess for sufficiently large $N$ this difference is simply negligible. This probably also is the reason for writing $N(n) \approx 2^{-n} N$ instead of "$=$".
Remark 2: In the process of $m$ getting $n$ quanta, there could also be "loops", in which $m$ first loses quanta and then gains them again or just does not participate in any quantum exchange processes, but since these loops may be arbitrarily long, for each step it does not matter if $m$ loses or gains or nothing happens, so their probability is just 1.
More sophisticated approach
If you have $N$ sites and want to distribute $NQ$ quanta on them, the number of possibilities to do this is
$$
\Omega(NQ) = \begin{pmatrix} NQ + N - 1 \\ N - 1 \end{pmatrix} \approx \begin{pmatrix} NQ + N \\ N \end{pmatrix} = \frac{(NQ + N)!}{(NQ)! N!}~,
$$
where the parts in parenthesis are binomial coefficients, not vectors. To see why this is true, encode the system's state as a sequence of dots and bars (e.g. "....|..|.|..."), where dots represent quanta and bars separate sites (in the example, site 1 has 4 quanta, site 2 has 2, site 3 has 1 and so on). Out of the $NQ + N - 1$ characters of the code, you have to choose $N-1$ to be bars, hence the binomial coefficient.
Now let's assume that site $m$ has $n$ quanta. The number of possibilities to distribute the remaining $NQ-n$ quanta on the remaining $N-1$ sites is
$$
\Omega_s(n) = \begin{pmatrix} NQ - n + N - 2 \\ N-2 \end{pmatrix} \approx \begin{pmatrix} NQ - n + N \\ N \end{pmatrix} = \frac{(NQ - n + N)!}{(NQ - n)! N!}~.
$$
Next we calculate
$$
\log \Omega_s(n) = \log(NQ - n + N)! - \log(NQ-n)! - \log N! \\
\overset{\text{Stirling's approximation}}\approx (NQ - n + N) \log(NQ - n + N) - NQ + n - N - (NQ - n) \log(NQ - n) + NQ - n - N \log N + N \\
=(NQ - n + N) \log(NQ - n + N) - (NQ - n) \log(NQ - n) - N \log N \\
=NQ \left( 1 - \underbrace{\frac{n}{NQ}}_{=: x}  + \frac 1 Q \right) \log \left( NQ \left( 1 - \frac{n}{NQ} + \frac 1Q \right) \right) - NQ \left( 1 - \frac{n}{NQ} \right) \log \left( NQ \left( 1 - \frac{n}{NQ} \right) \right) - N \log N \\
= NQ ( 1 - x + 1/Q) (\log(NQ) + \log(1-x+1/Q)) - NQ(1-x) (\log(NQ + \log(1-x)) - N \log N \\
\overset{\text{Taylor expansion in $x$ around 0}}= \log \Omega(NQ) + (-NQ \log(NQ + N) + NQ \log NQ)x + \mathcal O(x^2) \\
=\log \Omega(NQ) + \log \left( \frac{NQ}{NQ + N} \right) n + \mathcal O(x^2)~.
$$
Thus, there follows
$$
\log \frac{\Omega_s(n)}{\Omega(NQ)} = \log \Omega_s(n) - \log \Omega(NQ) = \log \left(\frac{NQ}{NQ + N} \right) n~,
$$
and exponentiation yields
$$
\frac{\Omega_s(n)}{\Omega(NQ)} = \exp \left( \log \left( \frac{NQ}{NQ + N} \right) n \right) = \left( \frac{NQ}{NQ + N} \right)^n~.
$$
Because all of the different possibilities in $\Omega_s(n)$ and $\Omega(NQ)$ are equally probable, this is the probability of finding site $m$ with $n$ quanta. If you set $Q = 1$, this becomes $(N / (2N))^n = 2^{-n}$, like before.
