I started with $N$ particles having random distribution in energy $E$ and isotropic directions.
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)$. Alternatively, I assumed that the pairing probability is not a constant and is determined by some weight, and selected pairs according to this weight.
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.
I repeat these steps (now remembering the directions) many times.
Key update
It looks like the key point is in the fact that I must use the weights $\propto |\mathbf{v}_{\text{rel}}|$, where $\mathbf{v}_{\text{rel}}$ is the relativistic relative velocity. I do not understand how this fact might affect the shape of the energy distribution (see Edit).
Simulation example
Also, the results do not depend on $N$ or the number of simulated steps. The system quickly falls to this "equilibrium" state, which seems like an intrinsic property of the pairing and interactions.
How I implement the weigths
This is how I implement the weights.
First, I make all possible combinations of pairs of 4-momenta. Then, I calculate the weights for these combinations. Then, I randomly sample the pairs according to the weights. Then, I select $N/2$ pairs out of the sampled $N(N-1)/2$ pairs so that each of the 4-momenta occurs only once. This is not a Hungarian algorithm -- the total weight of the selected pairs would not be maximal. The algorithm's timing grows as $N^{2}$.
Edit: some highlights from the literature
It turns out that my simulation is a simplified version of the direct simulation Monte Carlo approach (DSMC), where I restrict myself to one "particle cell". Some relevant discussions are provided in this paper; in particular, see Fig. 1 and the text around, which discusses typical errors. They have found in their simulation the shape $E^{2}\exp[-E/T]$. What they claim around Fig. 1 is that using the standard definition of the relative velocity $p_{1}\cdot p_{2}$ would lead to a wrong equilibrium distribution that looks like the distribution I obtained; they consider instead the so-called Moller velocity, which is a Lorentz-invariant generalization of the relative velocity assuming that the velocities $v_{1},v_{2}$ are not collinear:
$$ v_{\text{m}} = \sqrt{(\mathbf{v}_{1}-\mathbf{v}_{2})^{2} - (\mathbf v_{1}\times \mathbf v_{2})^{2}} $$
If I try to use this velocity as a weight of pairing, I indeed get the desired result:
It is interesting that if I include any particular Lorentz-invariant pre-factor to the velocity (say, the invariant mass of the colliding pair), the final result does not change. If I, however, include the Lorentz-violating factor (e.g., simply the product $E_{1}E_{2}$), then I again get a wrong shape.
However, I do not understand why the relative velocity matters at all. Namely, the probability of pairing would be $\propto \sigma v_{\text{m}}$. Since, by definition, $\sigma \propto v_{\text{m}}^{-1}$ (as the flux is proportional to $v_{m}$), the velocity seems to cancel out. Also, it is independent of the energy and depends only on the directions, so it is not clear to me how it would affect the shape of the energy distribution.