How to combine thermal velocity concept and special relativity? We know that in classical thermodynamics $$v_{rms} = \sqrt{\frac{3k_B T}{m}}$$ However we immediately see that this is wrong for high temperatures as there is no upper bound on velocity. How do I get the exact equation?
My approach- 
We have, $E = \sqrt{m_o^2c^4 + p^2c^2}$
Now from thermal energy we have total energy to be(sum of rest energy and thermal energy) $E = m_o c^2 + \frac{3}{2}k_B T$
Thus, $$m_oc^2 + \frac{3}{2}k_B T = \sqrt{m_o^2c^4 + p^2c^2}$$
Here,  $p = mv$ & $m = \frac{m_o}{\sqrt{1-\frac{v^2}{c^2}}} $  
Then we can solve for $v$ $ ( \sim v_{rms})$
I am not sure if this is right. Can someone correct me? Can you give me atleast the final result if not the entire drivation?
 A: What you are looking for is the mean of the Maxwell–Jüttner distribution, $$f(\gamma)=\frac{\gamma^2\beta}{\theta K_2(1/\theta)}e^{-\gamma/\theta}$$ where $\beta=v/c=\sqrt{1-1/\gamma^2}$, $\theta=k_BT/mc^2$ and $K_2$ is the Bessel function of the second kind. 
So the expected $\gamma$ would be $$E[\gamma]=\frac{1}{\theta K_2(1/\theta)}\int_1^\infty \gamma^3 \left (\sqrt{1-1/\gamma^2} \right )e^{-\gamma/\theta}d\gamma$$
Unfortunately there doesn't seem to be any closed form formula for it. Here is a plot of my numeric results in Matlab:

Here is a derivation of anisotropic Maxwell–Jüttner distributions, which may or may not be useful. 
A: The assumption that the thermal energy is $\frac{3}{2}k_bT$ is actually only valid at non-relativistic temperatures. In general we have to use the equipartition theorem to find the relation between temperature and energy:
\begin{equation}
\left< x_m \frac{\partial E}{\partial x_n} \right> = \delta_{mn} k_BT,
\end{equation}
where $x$ can be a coordinate or conjugate momentum.
Just taking the one-dimensional case for simplicity, in the Newtonian regime, $E = \frac{mv^2}{2}$, so that $v=\sqrt{\frac{k_BT}{m}}$. But in the relativistic case, 
$E = \sqrt{p^2c^2 + m_0^2c^4}$. This means that
\begin{equation}
\frac{c^2p^2}{\sqrt{p^2c^2 + m_0^2c^4}} = k_BT,
\end{equation}
so 
\begin{equation}
p^2 = \frac{k_B^2T^2c^2 \pm \sqrt{k_B^4T^4c^4 + 4k_B^2T^2m_0^2c^8}}{2c^4}.
\end{equation}
As $T \rightarrow \infty$ we get $p = k_BT/c$, or $E=k_BT$ (as the mass-term in the energy becomes negligible compared to the momentum). This is the well-known energy-equation for the ultra-relativistic gas. 
As $T \rightarrow 0$ we get $v = \sqrt{\frac{k_BT}{m_0}}$, the Newtonian result.
A: If I am correct, the Maxwell–Jüttner distribution has
$$E[\gamma]=\frac{1}{\theta K_2(1/\theta)}\int_1^\infty \gamma^3 \left (\sqrt{1-1/\gamma^2} \right )e^{-\gamma/\theta}d\gamma=\frac{K_1\left(\frac{1}{\theta}\right)}{K_2\left(\frac{1}{\theta}\right)}+3\theta,$$
where $K_n$ is a modified Bessel function of the second kind.
