A toy model of heat death Motivation:
This is a toy model of how a closed system will always evolve towards the distribution of maximal entropy, where no further transfer of heat/energy is possible.
Problem set up:
Fix a positive integer $N$, and denote by $[N]$ the set $\{1, \dots, N\}$.
Let $\mathcal L := [N] \times [N]$ be a 2D lattice.
We model a flow of heat as follows - Initially at time $0$, all heat is concentrated at $(1, 1)$. At each time step $t$, for $t \in \mathbb Z_+$, an element of $\mathcal L$ is picked uniformly at random. At that stage, it and all its immediate neighbours to its east, west, north and south average their heat.
Formally, we have a sequence of iid uniformly distributed $\mathcal L$-valued random variables $\varepsilon_n$, for $n \in \mathbb Z$, and a sequence $X_n$ of functions $\Omega \times \mathcal L \to [0, 1]$, defined as follows:
$X_0 = \mathbf 1_{(1, 1)}$, almost surely.
Assume now $X_0, \dots, X_n$ have already been defined.
Write $X_n = \sum_{z \in \mathcal L} \lambda_z \mathbf 1_{z}$ for (random) $\lambda_z \in [0, 1]$ with $\sum_{x \in \mathcal L} \lambda_z = 1$, and $\mathcal N(\varepsilon_n)$ for the set consisting of $\varepsilon_n$ and its two to four, depending on the location of $\varepsilon_n$, neighbours.
Then define $X_{n+1} := \left[\underset{x \in \mathcal N(\varepsilon_n)}{\sum}  \underset{y \in \mathcal N(\varepsilon_n)}{\sum}  \frac{\lambda_y}{|\mathcal N(\varepsilon_n)|} \mathbf 1_x \right]+ \underset{z \in \mathcal L \setminus \mathcal N(\varepsilon_n)}{\sum} \lambda_z \mathbf 1_z$.
where $|S|$ denotes the cardinality of a finite set $S$.

Question: Let $Y$ be the “uniform distribution of heat” given by $Y := \underset{z \in \mathcal L}{\sum} \frac{1}{|\mathcal L|} \mathbf 1_z$. Is it true that almost surely, $X_n \to Y$ uniformly?

Thus the system evolves almost surely toward a distribution where no further transfer of heat is possible.
 A: The Ehrenfest model (in discrete time, for simplicity) is just a Markov chain with the finite state space $\{0,1,\dots, N\}$ and the transition probabilities
$$
p(k,k-1)=k/N, \quad p(k,k+1)=1-k/N
$$
described by a single transition (averaging) operator $P$. Its stationary distribution $m$ is the binomial one with the parameters $\frac12,N$, and $\frac12 \theta (P^n+P^{n+1}) \to m$ for any initial distribution $\theta$. [Since the operator $P$ has period 2, one has to take the average of $\theta P^n$ and $\theta P^{n+1}$.]
In your situation there is a family of averaging operators $P_x$ indexed by the points from the state space $X$ which have a unique common invariant measure $m$ (the uniform distribution on $X$). You take a sequence $\boldsymbol x=(x_1,x_2,\dots)$ of iid $X$-valued uniformly distributed random variables, and ask whether, given an initial distribution $\theta$ on $X$, the sequence
$$
\theta P_{x_1} P_{x_2} \dots P_{x_n} 
$$
converges to $m$ almost surely. Note that since we are talking about measures on a finite state space, all reasonable kinds of convergence are equivalent (in particular, the $\ell^1$ convergence in the total variation norm $\|\cdot\|$ and the $\ell^\infty$ "uniform" convergence).
Let
$$
f_n(\boldsymbol x) = \| \theta P_{x_1} P_{x_2} \dots P_{x_n} - m \| \;.
$$
The sequence $f_n$ is non-increasing, and therefore convergent. By Kolmogorov's 0-1 law its limit $f_\infty$ is almost surely constant. Let $k$ be the minimal number such that for any $x\in X$ there is a sequence $x_1,x_2,\dots, x_k\in X$ with
$$
\text{supp}\,\delta_x P_{x_1} P_{x_2} \dots P_{x_k} = X \;.
$$
Then there is $\varepsilon > 0$ such that
$$
\mathbf E [ f_{n+k} | f_n ] \le (1-\varepsilon) f_n \qquad \forall\,n\ge 0 \;,
$$
whence $f_\infty=0$.
EDIT. The fact that the sequence $f_n$ is non-increasing is a consequence of the following inequality:
$$
\frac1d \sum_{i=1}^d |\theta_i - C| \ge \left| \frac1d \sum_i \theta_i - C \right| \;.
$$
Here $d$ is the cardinality of the averaging set (i.e., between 3 and 5 in your example), and $C=1/N^2$ is the common value of the weights of the uniform distribution. After removing $C$ and the division by $d$ the above inequality amounts to the well-known
$$
\sum_i |\theta_i| \ge \left| \sum_i\theta_i \right| \;.
$$
The expectation bound is just a constructive version of this inequality: if $f_n(\boldsymbol x)=F>0$, then there are two points $z_1,z_2\in X$ such that
$$
\theta P_{x_1} P_{x_2} \dots P_{x_n}(z_i) - m(z_i) \qquad i=1,2\;,
$$
have absolute values comparable with $F$ and opposite signs. Therefore by the definition of $k$ there is at least one choice of $x_{n+1},\dots, x_{n+k}$ with
$$
\begin{aligned}
&\|\theta P_{x_1} P_{x_2} \dots P_{x_n+k} - m \| \\ \\ &< (1-\epsilon) \cdot \|\theta P_{x_1} P_{x_2} \dots P_{x_n} - m \| \;,
\end{aligned}
$$
where $\epsilon$ is an appropriate constant (which only depends on $N$).
A: $\newcommand\R{\mathbb R}\newcommand{\ep}{\varepsilon}\newcommand{\de}{\delta}\newcommand{\ze}{\zeta}\newcommand{\om}{\omega}\newcommand{\Om}{\Omega}$
The answer is yes.
Indeed, let $\ze_n:=\ep_n$. Let $L:=\mathcal L$ and let $H$ be the Hilbert space of all functions $f\colon L\to\R$ with the norm $\|f\|:=\sqrt{\sum_{z\in L}f(z)^2}$, so that $\dim H=|L|=M:=N^2$. Then the $X_n$'s may be viewed as random vectors in $H$. Moreover, for all $n=0,1,\dots$,
\begin{equation}
    X_n=P_{\ze_n}\cdots P_{\ze_1}X_0,\tag{1}
\end{equation}
where, for each $z\in L$, $P_z$ is the orthogonal projector (o.p.) of $H$ onto the space of all functions $f\in H$ that are harmonic at $z$. We say that a function $f\in H$ is harmonic at $z$ if $f(z)$ coincides with the average of the values of $f$ at all neighbors of $z$.
Let us also say that a function $f\in H$ is harmonic if its harmonic at every $z\in L$.
Remark 1: Every harmonic function $f\in H$ is constant. Indeed, such a function $f$ attains its maximum $M$ at some point $z\in L$. Then $f(y)=M$ for all neighbors $y$ of $z$. Since any point of $L$ can be connected to any other point of $L$ by a chain linking neighbors, we conclude that $f$ is the constant $M$ on $L$. $\quad\Box$
It follows from (1) that $\|X_0\|\ge\|X_1\|\ge\cdots$ and hence there is a (random) limit
\begin{equation}
    l:=\lim_n\|X_n\|^2. \tag{2}
\end{equation}
Next, for each real $\ep>0$ there is some real $\de(\ep)>0$ such that the following implication holds:
\begin{equation}
    \|f\|\le1\ \&\ \max_{z\in L}\|Q_z f\|\le\de(\ep)\implies \|Qf\|<\ep, \tag{3}
\end{equation}
where $Q_z:=I-P_z$, $Q:=I-P$, $I$ is the identity operator on $H$, and $P$ is the o.p. of $H$ onto the space of all constant functions in $H$.
Indeed, otherwise there would exist some real $\ep>0$ and a sequence $(f_k)$ in $H$ such that   $\|f_k\|\le1$, $\max_{z\in L}\|Q_z f_k\|\le1/k$, and $\|Qf_k\|\ge\ep$ for all $k$. By compactness, without loss of generality $f_k\to f$ for some $f\in H$, so that $\max_{z\in L}\|Q_z f\|=0$ and $\|Qf\|\ge\ep$, so that $f$ is a non-constant harmonic function, which contradicts Remark 1.
Note that for all $n=0,1,\dots$
\begin{equation}
    1\ge\|X_n\|^2\ge \|PX_n\|^2=1/M, \tag{4}
\end{equation}
\begin{equation}
    \|X_n\|^2-\|X_{n+1}\|^2=\|Q_{\ze_{n+1}}X_n\|^2. \tag{6}
\end{equation}
On the event $\{l>1/M\}$, for all $n$ we have $\|QX_n\|^2=\|X_n\|^2-\|PX_n\|^2=\|X_n\|^2-1/M\ge d:=l-1/M>0$ and hence, by (3),
\begin{equation}
    \max_{z\in L}\|Q_z X_n\|>\de_*:=\de(\sqrt{d})>0. \tag{8}
\end{equation}
So, there is a sequence $(g_n)$ of (deterministic) functions $g_n\colon L^n\to L$ such that on the event $\{l>1/M\}$
\begin{equation}
    \|Q_{Z_n} X_n\|>\de_*>0 \tag{9}
\end{equation}
for all $n$, where $Z_n:=g_n(\ze_1,\dots,\ze_n)$.
For natural numbers $n,m$, let
\begin{equation}
    p_{n,m}:=P(\ze_{n+1}\ne Z_n,\dots,\ze_{n+m+1}\ne Z_{n+m}). \tag{10}
\end{equation}
Then, conditioning on $\ze_1,\dots,\ze_{n+m}$, we get $p_{n,m}=(1-1/M)p_{n,m-1}$ if $m\ge2$ and hence $p_{n,m}\le(1-1/M)^m\to0$ as $m\to\infty$. It follows that the events $\{\ze_{n+1}= Z_n\}$ occur for all natural $n$ in a random set $\nu$ which is infinite almost surely (a.s.).
Thus, in view of (6) and (9), on the event $\{l>1/M\}$ we a.s. have
\begin{equation}
    1=\|X_0\|^2\ge\sum_{n=0}^\infty(\|X_n\|^2-\|X_{n+1}\|^2) \\
    \ge\sum_{n\in\nu}(\|X_n\|^2-\|X_{n+1}\|^2) \\ 
        =\sum_{n\in\nu}\|Q_{\ze_{n+1}}X_n\|^2 \\ 
        =\sum_{n\in\nu}\|Q_{Z_n}X_n\|^2 \\
        \ge \sum_{n\in\nu}\de_*=\infty.
\end{equation}
Thus, the event $\{l>1/M\}$ has probability $0$. That is, $l=1/M$ a.s. So, $\|X_n\|^2\to1/M$ a.s. So, $\|X_n-PX_n\|^2=\|X_n\|^2-\|PX_n\|^2=\|X_n\|^2-1/M\to0$ a.s. Since $PX_n$ is the constant function equal $1/M$, we have proved the desired result: $X_n\to1/M$ a.s.
A: After t steps, Let $$S_t=S(X_t) =\sum_{ z \in L}  \lambda_z(t)^2] - 1/N^2=\sum_{ z \in L} [\lambda_z(t) -1/N^2]^2 \,.$$,   Then $S_t$ is nonincreasing, so it must converge to some limiting random variable S*. Moreover, if $S_t>0$ then with positive probability $S_{t+1}<(1-c)S_t$ for some $c=c(N)$. This implies that that $S_*=0$ almost surely, whence $\lambda_z(t)$ must tend to $1/N^2$ almost surely for each $z$.
Details: If $S_t=a^2$, then there must exist vertices $z,w$ such that $[\lambda_z(t)-\lambda_w(t)] \ge a/N$, so there must exist two adjacent vertices $u,v$ such that $[\lambda_u(t)-\lambda_v(t)] \ge a/(2N^2)$. If $u$ or $v$ is the next vertex selected, then $S_{t+1} \le S_t-a^2/(8N^4)$. In other words, given the history up to time $t$, with probability at least $2/N^2$, we have $S_{t+1} \le S_t [1-1/(8N^4)]$. If $S*>0$ with positive probability, then on this event, find a random $t$ such that $S_t<S*[1+1/(8N^4)]$ and then wait until the first time one of $u,v$ is selected to get a contradiction.
