I made a simulation that studies the thermalization of massless particles, assuming isotropic and homogeneous spatial distribution, in 3 dimensions. Namely, I started with $N$ particles having random energy distribution and isotropic directions assuming that they meet at one point (so I drop any spatial coordinates from the exercise), then randomly chose their pairs, assuming the probability of selecting each pair is $1/N$, and then generated collision kinematics: I boosted to the CM frame, generated random polar (assuming that $\cos(\theta)$ is uniformly distributed from -1 to 1) and azimuthal scattering angles, and boosted back. I repeat these steps (now remembering the directions) many times.

I expected the final energy distribution of the particles to follow the scaling $E^{2}\exp[-E/T]$, where $T$ is the effective temperature of the system determined from $N$ and the total energy ($E^{2}$ comes from the Liouville theorem). However, I got $E\exp[-E/T]$. When repeating the same exercise for the 2D case, I got $\exp[-E/T]$ instead of $E\exp[-E/T]$. So it looks like I missed some important point, which leads to the reduction of the dimensionality by one.

Below, there is an example of the simulation for the 3D case - the histogram distribution of the initial state (yellow bars), the final state (blue bars), as well as fits of the distributions $E^{n}\exp[-E/T]$ with $n = 0$ (the blue curve), $n = 1$ (the red one), $n = 2$ (the green one). For the initial distribution, I simply took $E^{2}\exp[-E/T]$:

What may be the reason for this?

Edit

In another simulation, I assumed that each pair has the interaction weight $\sim 1/p_{1}^{\mu}p_{2\mu} \sim (1-\cos(\alpha))^{-1}$, where $\alpha$ is the angle between the colliding particles. It comes form the definition of the Lorentz-covariant cross-section. Then I, surprisingly, get the pure Boltzmann exponent independently of the dimensionality of the system:

It sounds weird: in reality, the weight also includes the squared matrix element of the scattering, which may take any form, and in particular be proportional to $p_{1}^{\mu}p_{2\mu}$, but the final distribution should not care about the explicit form of the element.