This animation shows the evolution of the energy distributions of a 2D gas that initially
has a mono-energetic distribution (MED where all atoms have initial energy `E_0`).
Let the indices `n` refer to the interaction generation number and 1 and 2 refer to the particle numbers.
That means that the energies of all generations of the atoms will be random values
in the range between 0 and `E_[n,1]+E[n,2]` or

and where r=random(0,1)
and random(0,x) provides a randomly selected floating point number
between 0 and x. Note that the generation n+1 energy values are chosen so that neither is negative.
This non-negativity is the most important property driving the final distribution since it results in
a compaction at lower energies and a spreading at higher energies which are unrestricted in
value except that they must obey energy conservation on each interaction. The distribution after many generations always fits an exponential
`N(E)=N_0*exp(-E/E_(avg))`
where `N_0` is a normalizing constant and `E_(avg)` means the average value of `E`.
This would seem to explain energy distributions for one dimension (1D) and three dimensions (3D) as well.
Unfortunately this is not true for the other dimensions.
To explain 3D energy distributions, one has to use the full conservation of momentum equations
as well as the energy conservation equations.
The result for 3D is a distribution of form `sqrt(E/E_(avg))exp[E/(2/3E_(avg))]`.
The individual momentum probabilities are all like `f(px)dpx=exp(-(p_x*p_x)/[p_x*p_x]_(avg))dpx`
and the 3D momentum probability density is `dF(p_x,p_y,p_z)=f(p_x)f(p_y)f(p_z)dp_xdp_ydp_z`.
Because `p^2dp=m^(3/2)sqrt(E)dE` where m is the mass, we can change the momentum
differential to energy differential and rewrite this as `G(E)=exp(-p^2/[p^2]_(avg))p*p*dp=exp(-E/E_(avg))sqrt(E/E_(avg))dE`
thus getting the expression used to fit the 3D energy here.
Using the same collision equations elevated to higher dimensions, we could model particles with access
to 4 or more dimensions. Of course, these higher dimensions usually refer to rotation and vibration modes.

For 1D the situation is even more complicated for 2 reasons:
First, the particles can interact only with their nearest initial neighbors.
Second, there is no orthogonal dimension (like x,y) for the momentum to be transferred to during a collision.
As a result, the momentum distribution does not approximate a gaussian `exp(-p^2/([p*p]_(avg)))` which is a
correct distribution for, say, the y component of the momentum in 2D. So I have used a function involving
the square root of the energy to describe the 1D energy distribution.
Converting the momentum differential, dp, to energy we would have dp=dE/sqrt(E) which is what
I have used for 1 dimension energy distributions.
It should be noted that, while the energy distributions for 2D and 3D are mono-energetic,
the initial distribution for 1D is flat. with the same average energy as for 2D and 3D.

You will motice that I have not included particle motion graphics here. The reason is that the graphics take a lot of computer time.
The scattering partner particle is randomly selected from the group of all particles of the same dimensional group.