Why does energy exchange mechanism appear to affect energy distribution in a gas?

Question

I'm trying to reproduce Boltzmann distribution by simulating an abstract gas, i.e. a system of particle-like objects that can randomly exchange energy between each other.

The simulation is arranged as follows. The state of the gas is a collection $\mathcal E$ of energies of $N$ particles:

$$\mathcal E=\{E_n\}_{n=1}^N.\tag1$$

Initially, $E_n=E\delta_{n,1}$, where $\delta_{n,m}$ is Kronecker delta.

At each iteration a pair of particles is chosen by generating two random uniformly distributed in $[1,N]$ integers $n,m$, such that $n\ne m$ (they are generated independently and regenerated in case of $n=m$). The total energy $E_\Sigma=E_n+E_m$ is computed and redistributed using a random real $\alpha\in[0,1]$ (uniformly distributed) as

$$ \begin{align} E'_n&=E_\Sigma\alpha,\\ E'_m&=E_\Sigma(1-\alpha). \end{align}\tag2 $$

These new energies are written in place of old ones into $\mathcal E$. After many iterations ($1.5\times10^7$ iterations for $10^5$ particles), $\mathcal E$ is assumed to have reached equilibrium, and its histogram is plotted.

Here's one such result, with 40 histogram bins plotted (joined blue points) against the Boltzmann distribution $A\exp\left(-\frac{\varepsilon}{\langle \mathcal E\rangle}\right)$ (orange curve):

Here the result nicely matches the Boltzmann distribution.

Now, instead of $(2)$, I choose a different mechanism of energy exchange:

$$ \begin{align} E'_n&=E_\Sigma\alpha^k,\\ E'_m&=E_\Sigma(1-\alpha^k). \end{align}\tag3 $$

Setting $k=2$, I get the following result:

As you can see, now there is a discrepancy between the histogram and the ideal Boltzmann curve. I've tried increasing number of iterations tenfold, but the result doesn't change, so apparently, the system has come to its equilibrium.

Increasing $k$ further, e.g. to $k=3$, increases the discrepancy even more:

This makes me wonder: I thought that regardless of the mechanism of the energy exchange, the equilibrium energy distribution should be the same. But for some reason the mechanism does affect the result. How can then Boltzmann distribution work for actual gases, which all have different atomic repulsive potentials (different hardness at least)? We don't see a large deviation from this distribution there, do we?

If the energy exchange mechanism shouldn't affect the outcome, then what could be wrong in my simulation?

Answer 1

The Boltzmann distribution is derived from the canonical ensemble, which in turn is derived by assuming your system has ergodic Hamiltonian dynamics on phase space. That in turn is true for a wide variety of real-world systems.

If you discard all of this structure and simply postulate an arbitrary time evolution, you shouldn't expect the Boltzmann distribution to hold. It would simply be asking too much for such a simple result to hold for all possible dynamical systems!

For example, in the limit $k \to \infty$ we have $\alpha^k \to 0$, which means that in any interaction, all of the energy goes to the latter particle. After a long time, all of the particles end up with zero energy except for one, which has all the energy. This is absolutely not a Hamiltonian time evolution (the phase space volume shrinks to zero), and the final result is absolutely not Boltzmann distributed.