I made a simulation that studies the thermalization of massless particles, assuming isotropic and homogeneous spatial distribution, in 3 dimensions. Namely, 

1. I started with $N$ particles having random distribution in energy $E$ and isotropic directions.
 
2. Next, I assumed that they meet at one point (so I drop any spatial coordinates from the exercise) and randomly chose their pairs, assuming the probability of selecting each pair is pair-independent and equal to $2/N(N-1)$. 

3. Then, I generated collision kinematics: I boosted to the CM frame, and generated random polar (assuming that $\cos(\theta)$ is uniformly distributed from -1 to 1, which is true in the case of a constant matrix element of the process) and azimuthal scattering angles (simply from $-\pi$ to $\pi$), and boosted back. 

4. 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: the density of states in a D-dimensional space is $d^{n}\mathbf{p} \propto E^{n-1}dE$. 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. What may be the reason for this?

**Simulation example**

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]$:

[![enter image description here][1]][1]

It is clearly visible that the interactions turn the distribution from $E^{2}\exp[-E/T]$ to $E\exp[-E/\tilde{T}]$.

P.S. The results do not depend on $N$ as well as on the number of simulated steps. The system falls down to this "equilibrium" state very fast. It looks like an intrinsic property of the pairing and interactions.

**Edit: what if assuming that the pairing probability is not constant?**

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 1,2. 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:

 [![enter image description here][2]][2]

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. I checked, however, many random weights, with the only property that they are $\propto (p_{1}^{\mu}p_{2\mu})^{-1}$, and they all lead to the same Boltzmann shape. 


  [1]: https://i.sstatic.net/0wPiy.png
  [2]: https://i.sstatic.net/xkGZO.png