Is there a relativistic generalization of the Maxwell-Boltzmann velocity distribution? The Maxwell-Boltzmann velocity distribution in 3D space is
$$
f(v)dv = 4\pi \left(\frac{m}{2\pi k_B T}\right)^{3/2} v^2 \exp\left(-\frac{m v^2}{2k_B T}\right)dv$$
It gives the probability for a single particle to have a speed in the intervall $[v,v+dv]$. But this probability is not zero for speeds $v > c $ in conflict with special relativity.
Is there a generalization of the Maxwell-Boltzmann velocity distribution which is valid also in the relativistic regime so that $f(v) = 0$ for $v>c$ ? And how can it be derived? Or can a single particle distribution simply not exist for relativistic speeds, because for high energies, we always have pair-production meaning the particle number is not conserved and a single particle distribution can not be defined in a consistent way? 
 A: The point of a Boltzmann distribution is that it maximizes entropy given a fixed energy.  The concept applies to systems with other degrees of freedom besides translational kinetic energy.  The general distribution, from Wikipedia is

Thus, the simple adjustment to the Maxwell-Boltzmann distribution you cited is to replace the Newtonian kinetic energy $\frac{mv^2}{2}$ with the relativistic kinetic energy $(\gamma - 1)mc^2$ everywhere it appears in the distribution.
Pair creation is a separate issue that I'll leave to someone else.
A: The usual, or text book generalisation (Juttner), which was first derived in 1911 is not covariant. Presumably a fully covariant distribution would cover all these details assuming that it exists in the first place. The most recent attempt I am aware of is by Ewald Lehmann who goes back to basics, "Covariant equilibrium statistical mechanics", Journal of Mathematical Physics 47, 023303,2006.
David Sher
A: Here's the relativistic generalization:
$$f(\beta)=\left(\frac{\pi\alpha\beta^2}{2(1-\beta^2)}\right)^{(d-1)/2}\frac{e^{-\alpha/\sqrt{1-\beta^2}}}{2\sqrt{1-\beta^2}\,K_{(d+1)/2}(\alpha)}$$
where $\alpha=\frac{mc^2}{kT}$, $d$ is the spatial dimension, and $K$ is a Bessel function. I'll leave the derivation to the reader, but this integral is handy:
$$\int_0^\infty e^{-\alpha\cosh\theta}\sinh^{2\nu}\theta\,d\theta=\frac{1}{\sqrt{\pi}}\left(\frac{2}{\alpha}\right)^\nu\Gamma(\nu+1/2)\,K_\nu(\alpha)$$
For $d=3$ I've included a plot with several values of $\alpha$.

A: Per Mark's request, I'll provide an answer.
First, I think it's not possible to obtain such a formula. One can of course naively try to extrapolate various classical formulae but all of these attempts are bound to fail. Here's why: even when constructing relativistic quantum mechanics one finds that the theory is not consistent. For example, in non-relativistic quantum mechanics one has a position operator that can be used to get information on precise position of the particle (up to the uncertainty given by Heisenberg's uncertainty principle). But when one includes relativity into the picture the theory stops being consistent. This reflects the fact that to localize the particle with great precision one has to make experiments with increasingly high energy and at certain point the energy is enough for the appearance of new particles. In fact, creation and annihilation of particles is inevitable in relativistic regime (and in some systems it's not even clear what particles should be and one has to talk about fields instead). The situation is even more pronounced in statistical physics where there are huge number of particles present.
More importantly, there's no need to get such a formula. Consider systems for which it would be useful. Such systems would have to be extremely non-classical (like neutron stars, black holes interior, quark-gluon plasma, etc.) and the concept of velocity would have no meaning as there's no way to observe individual particles of these systems; which is in contrast to the classical case where you can test Maxwell distribution by letting the particles out of the box one by one and seeing how fast they are (the actual experiment is lot more sophisticated of course, but that's unimportant here).
A: There is no particle count invariance in general relativity. One can substitute energy density and from a classical point of view an equilibrium energy density defines a specific temperature. Lehmann's paper, which I referred to in a comment posted, above clearly attempts to do this. From a practical point of view such extreme conditions preclude any attempt at experimental verification. Numerical simulations usually refer to impenetrable classical particles and disregard loss of count invariance.
Richard Feynman is reported to have said something to the effect that if he can derive a result in several different ways he is probably right. Lehmann's results can be duplicated, at least qualitatively, using Maxwell's original (teenage) argument based on symmetry alone. Both give distributions which become flatter and both have a critical temperature of logT 11-12, roughly what one may expect at the heart of a core collapse supernova.
An upper limit should not be a surprise. For example if we go to more and more extreme conditions we will reach those prevalent in the early universe at the time of the big bang. Any hope of even a theoretical equilibrium is likely to be lost well before this. 
